A Novel Development in Three-Dimensional Analytical Solutions for Air Pollution Dispersion Modeling

Abstract

A novel analytical solution has been developed for the three-dimensional dispersion of atmospheric pollutants, providing a new approach to understanding and addressing this important environmental issue. The central concept of the study is to divide the planetary boundary layer into multiple vertical sub-layers, each characterized by its own average wind speed and eddy diffusivity. This allows for a more comprehensive and nuanced examination of atmospheric processes within the boundary layer. The validity of the model is thoroughly evaluated through a comparison of its predictions with data collected from the Copenhagen Diffusion and Prairie Grass experiments. This approach ensures that the model accurately reflects the complexities of atmospheric dispersion in real-world scenarios. The results of the study demonstrate a strong correlation between the predicted and measured crosswind-integrated concentrations. Furthermore, the statistical indices computed for the model fall within an acceptable range, indicating a high level of accuracy in the model’s predictions. These findings reinforce the validity of the analytical solution for modeling atmospheric pollutant dispersion.

1. Introduction

The scientific community has been driven to act urgently due to the increase in ecological disasters worldwide. A key approach to addressing this issue is the development of reliable models that enable quantitative prediction of pollution-related phenomena, either through analytical descriptions or simulations using powerful and operational tools. These models, which can be based on simulations or analytical descriptions, have strong quantitative predictive power for understanding pollution-related phenomena. The atmosphere is considered the primary means of dispersing pollutants in the environment, which can come from industrial sources or accidental events, leading to progressive contamination of the ecosystem, fauna, flora, and all populations. Therefore, accurately evaluating how pollutants move through the atmosphere in the boundary layer is vital to preserve, protect, and restore the integrity of ecosystems. To accomplish this, it is essential to create atmospheric dispersion models tailored to the parameters and weather conditions specific to the region being studied, using parameters, weather conditions, and local topographic information. These models should produce realistic outcomes on environmental consequences, helping to minimize the negative impact of potential disasters like forest fires. Furthermore, implementing simulations of models based on specific cases can be beneficial in establishing particular emission limits for industrial sites, limiting the release of pollutants into the air.

The aim of our work is to propose a closed-form analytical solution for three-dimensional advection-dispersion transport problems in finite, multilayered media using rigorous mathematical tools. We ensure that wind velocity profiles and vertical diffusivity coefficients take average values in each sub-layer. We have made innovative contributions by using a governing function to generate the eigenvalues associated with our problem and overcoming the common difficulty of missing some of the eigenvalues when they are calculated by developing a transcendental equation for each layer. Finally, our approach has enabled us to obtain a closed-form analytical solution that differs from previous solutions that required the determination of integral coefficients. This method could have important implications in many scientific and technical fields for solving complex advection-dispersion transport problems in finite, multilayered media.

In recent years, there has been an increasing emphasis on creating new analytical approaches that can be used for a range of wind speeds and turbulent diffusivity coefficients. For the most part, these approaches involve the projection of the solution onto a basis of orthogonal polynomials, such as the GILTT technique [1, 2, 3, 4, 5, 6, 7, 8]. However, an important drawback of this method is that it requires a large number of eigenvalues to ensure convergence, which can reach up to 250 eigenvalues. In contrast, our new solution is able to provide better results by using a much smaller number of eigenvalues, specifically between 10 and 15, to ensure convergence. This is a significant improvement over existing approaches, making our method more efficient and practical for analytical applications.

Our model was developed using the Fourier transform and separation of variables technique, which resulted in the Sturm-Liouville problem. In order to understand the factors influencing pollutant dispersion, we considered the following parameters: (i) the wind speed profile of Deaves and Harris [9], (ii) the vertical turbulent diffusivity coefficient, which is considered as an explicit function of the downwind distance and vertical height under convective conditions, as described by Mooney and Wilson [10] and Degrazia et al. [11], and (iii) the lateral eddy diffusivity coefficient, which also depends on the downwind distance and vertical height, as described in the work of Huang [12] and Brown et al. [13]. We first proceed to the exposition of the general formulation of the equation that governs the dispersion of pollutants in the atmospheric boundary layer, with a comprehensive presentation in Section 2. We then turn to the explicit solution in Section 3, while exploring the model parameterizations in Section 4. We provide an in-depth discussion of our numerical results in Section 5. Finally, the last section proposes our conclusion, a synthesis of our investigation.

2. Presentation of the general problem

The movement of pollutants within the planetary boundary layer (PBL) is characterized by turbulent dispersion, which is governed by a mathematical equation known as the advection-diffusion equation. This equation considers the combined influence of two primary mechanisms, advection and diffusion, which respectively refer to the transport of pollutants by wind and the spreading of pollutants due to turbulence.

By applying the advection-diffusion equation, researchers can gain a more comprehensive understanding of the intricate and dynamic behavior of pollutants within the PBL. This equation serves as a mathematical framework that accounts for various factors that affect pollutant dispersion, such as wind direction and speed, Eddy diffusivities, and atmospheric stability, and pollutant concentration, presented in the following form:

$$\frac{\partial C}{\partial t}+\nabla \mathbf{.}\left({\mathbf{U}}_{\mathbf{w}}C\right)=\nabla \mathbf{.}\left(D\nabla C\right)+S,$$E1

where ${\mathbf{U}}_{\mathbf{w}}={\left(U,V,W\right)}^T$ is the wind speed vector (m/s) representing the components U, V, and W in the east-west, north-south and vertical directions, respectively; D is the molecular diffusion coefficient; S is the source term; and $\nabla$ is the gradient operator.

By use of the time average and fluctuation values, $U=u+u^{\prime }$, $V=v+v^{\prime }$, $W=w+w^{\prime }$ and $C=c+c^{\prime }$, the wind speed vector ${\mathbf{U}}_{\mathbf{w}}$ is expressed as:

$${\mathbf{U}}_{\mathbf{w}}=\overline{{\mathbf{U}}_{\mathbf{w}}}+{{\mathbf{U}}_{\mathbf{w}}}^{\text{'}},\text{with}\overline{{\mathbf{U}}_{\mathbf{w}}}={\left(u,v,w\right)}^T\text{and}{{\mathbf{U}}_{\mathbf{w}}}^{\text{'}}={\left(u^{\prime },v^{\prime },w^{\prime }\right)}^T$$E2

Applying the Reynolds averaging rules [14] to the vertical mass flow term, denoted as ${\mathbf{U}}_{\mathbf{w}}$C, results in the following expression: It turns out that turbulent diffusion can be described with Fick’s laws of diffusion as follows [15].

$$\overline{u^{\prime }{c}^{\prime }}=-K_x\frac{\partial c}{\partial x};\overline{v^{\prime }{c}^{\prime }}=-K_y\frac{\partial c}{\partial y};\overline{w^{\prime }{c}^{\prime }}=-K_z\frac{\partial c}{\partial z}.$$E3

where $K_x$, $K_y$, and $K_z$ are the eddy diffusivities components along x, y, and z directions, respectively.

It should be noted that in a turbulent boundary layer where advection is occurring, K will be larger than D and eddy diffusion will dominate solute transport. In this case, the molecular diffusion coefficient $\nabla .\left(D\nabla C\right)$ is then to be replaced by an eddy or turbulent diffusivity. The source term could be eliminated from Eq. (1) and should be added to the boundary conditions as a delta function: At the point $\left(0,0,H_s\right)$, there is a source rejecting the pollutant with a continuous flow Q:

$$uc\left(0,y,z\right)=Q\delta \left(y\right)\otimes \delta \left(z-H_s\right),$$E4

where $\otimes$ is the tensor product of two distributions and $H_s$ is the source height. By application of the Reynolds averaging and the divergence operator to Eq. (3), Eq. (1) may be written as follows.

$$\frac{\partial c}{\partial t}=-\frac{\partial }{\partial x}\overline{u^{\prime }{c}^{\prime }}-\frac{\partial }{\partial y}\overline{v^{\prime }{c}^{\prime }}-\frac{\partial }{\partial z}\overline{w^{\prime }{c}^{\prime }}-\overline{U_w}.\nabla c.$$E5

In the remainder of this paper, the following assumptions are considered:

1. the steady state condition (i.e., $\frac{\partial c}{\partial t}=0$);

2. the two terms $v\frac{\partial c}{\partial y}$ and $w\frac{\partial c}{\partial z}$ are neglected since the x-axis coincides with the wind flow average, therefore the w and v wind velocity components are less important; and

3. the turbulent diffusion in the direction of the mean wind is neglected compared to the advection transport mechanism (i.e., $u\frac{\partial c}{\partial x}>>\frac{\partial }{\partial x}\left(K_x\frac{\partial c}{\partial x}\right)$).

These assumptions lead to the steady-state advection-diffusion equation defined as:

$$\left\{\begin{array}{c}u\left(z\right)\frac{\partial c}{\partial x}=\frac{\partial }{\partial y}\left(K_y\frac{\partial c}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_z\frac{\partial c}{\partial z}\right)\mathbf{;}\left(x,y,z\right)\in ]0,L_x\left[\times \right]-L_y,L_y\left[\times \right]0,H_{\mathit{mix}}[,\\ \\ \text{subject to the boundary conditions}:\left\{\begin{array}{ll}u\left(z\right)c\left(x,y,z\right)=Q\cdot \delta \left(y\right)\otimes \delta \left(z-H_s\right);& x=0,\\ K_y\left(x,z\right)\frac{\partial c\left(x,y,z\right)}{\partial y}\to 0;& \mid y\mid \to L_y,\\ K_z\left(x,z\right)\frac{\partial c\left(x,y,z\right)}{\partial z}=0;& z\in \left\{0,H_{\mathit{mix}}\right\}.\end{array}\right.\end{array}\right.$$E6

where $z_0$ is the surface roughness length and $H_{\mathit{mix}}$ is the PBL height.

We consider that the eddy diffusivities have the following separable formulations:

$$K_y\left(x,z\right)={\zeta }_y\left(x\right)u\left(z\right),$$E7

$$K_z\left(x,z\right)=\xi \left(x\right){\phi }_z\left(z\right).$$E8

We vertically divide the PBL into H intervals such that for each one the eddy diffusivity and wind speed assume average values. For $h=1,\cdots ,H$,

$$u_h=\frac{1}{z_h-z_{h-1}}{\int }_{z_{h-1}}^{z_h}u\left(s\right)\mathit{ds},$$E9

$${\phi }_{z_h}=\frac{1}{z_h-z_{h-1}}{\int }_{z_{h-1}}^{z_h}{\phi }_z\left(s\right)\mathit{ds}.$$E10

By use of the formulations of $K_y$ and $K_z$ given by Eqs. (7) and (8), Eq. (5) is written as:

$$u_h\frac{\partial c_h}{\partial x}={\zeta }_y\left(x\right)u_h\frac{{\partial }^2{c}_h}{\partial y^2}+\xi \left(x\right){\phi }_{z_h}\frac{{\partial }^2{c}_h}{\partial z^2},$$E11

with $u_h$ and ${\phi }_{z_h}$ (given by Eqs. (9) and (10)) are constants.

This later Eq. (11) is subject to the first boundary conditions Eq. (6) on the one hand, and on the other hand, the continuity of both the concentration and the flux at the interface level is applied.

$$\left\{\begin{array}{l}{\phi }_{z_1}\frac{\partial c_1\left(x,y,z_0\right)}{\partial z}=0,\\ \\ \left\{\begin{array}{ll}{c}_{h-1}\left(x,y,z_{h-1}\right)=c_h\left(x,y,z_{h-1}\right)& \\ & ,\forall h\in \left\{2,\cdots ,H\right\},\\ {\phi }_{z_{h-1}}\frac{\partial c_{h-1}\left(x,y,z_{h-1}\right)}{\partial z}={\phi }_{z_h}\frac{\partial c_h\left(x,y,z_{h-1}\right)}{\partial z}& \end{array}\right.\\ \\ {\phi }_{z_H}\frac{\partial c_H\left(x,y,z_H\right)}{\partial z}=0\end{array}\right.$$

3. A mathematical approach to solving the pollutant dispersion equation: analytical solution

We start this section by applying Fourier transform to Eq. (11). Let ${\widehat{c}}_h^{\omega }\left(x,z\right)$ denote the Fourier transformation of $c_h$ with respect to y.

$${\widehat{c}}_h^{\omega }\left(x,z\right)={\int }_{-\infty }^{+\infty }{c}_h\left(x,y,z\right)e^{-2\mathit{\text{iπωy}}}\mathit{dy},h\in \left\{1,\cdots ,H\right\}$$

which gives

$$u_h\left[\frac{\partial {\widehat{c}}_h^{\omega }}{\partial x}+{\left(2\pi \right)}^2{\omega }^2{\zeta }_y\left(x\right){\widehat{c}}_h^{\omega }\right]=\xi \left(x\right){\phi }_{z_h}\frac{{\partial }^2{\widehat{c}}_h^{\omega }}{\partial z^2},z\in \left[z_{h-1},z_h\right].$$E12

Let

$${\chi }_h^{\omega }\left(x,z\right)={\widehat{c}}_h^{\omega }\left(x,z\right).exp\left({\left(2\pi \right)}^2{\omega }^2{\int }_0^x{\zeta }_y\left(s\right)\mathit{ds}\right)$$E13

then

$$\frac{\partial {\chi }_h^{\omega }}{\partial x}=\left[\frac{\partial {\widehat{c}}_h^{\omega }}{\partial x}+{\left(2\pi \right)}^2{\omega }^2{\zeta }_y\left(x\right){\widehat{c}}_h^{\omega }\right]exp\left({\left(2\pi \right)}^2{\omega }^2{\int }_0^x{\zeta }_y\left(s\right)\mathit{ds}\right).$$

By multiplying both sides of Eq. (12) by $exp\left({\left(2\pi \right)}^2{\omega }^2{\int }_0^x{\zeta }_y\left(s\right)\mathit{ds}\right)$ and since $u_h$ and ${\phi }_{z_h}$ are constants for each interval, we show easily that, for all $h\in \left\{1,\cdots ,H\right\}$:

$$u_h\frac{\partial {\chi }_h^{\omega }}{\partial x}=\xi \left(x\right){\phi }_{z_h}\frac{{\partial }^2{\chi }_h^{\omega }}{\partial z^2},z\in \left[z_{h-1},z_h\right]$$E14

We proceed in the same way with the boundary conditions, and find

$$\left\{\begin{array}{l}{\phi }_{z_1}\frac{\partial {\chi }_1^{\omega }\left(x,z_0\right)}{\partial z}=0\\ \begin{array}{l}\left\{\begin{array}{l}{\chi }_{h-1}^{\omega }\left(x,z_{h-1}\right)={\chi }_h^{\omega }\left(x,z_{h-1}\right)\\ {\phi }_{z_{h-1}}\frac{\partial {\chi }_{h-1}^{\omega }\left(x,z_{h-1}\right)}{\partial z}={\phi }_{z_h}\frac{\partial {\chi }_h^{\omega }\left(x,z_{h-1}\right)}{\partial z},h\in \left\{2,\cdots ,H\right\}\end{array}\right.\\ {\phi }_{z_H}\frac{\partial {\chi }_H^{\omega }\left(x,z_H\right)}{\partial z}=0\end{array}\end{array}\right.$$

The solution of Eq. (14) is assumed to be in the form:

$${\chi }_h^{\omega }\left(x,z\right)=\sum _{n=0}^{\infty }{G}_{h,n}^{\omega }\left(x\right)P_{h,n}\left(z\right),h\in \left\{1,\cdots ,H\right\}$$

This separated form gives two ordinary differential equations to be solved:

$$\frac{d{G}_{h,n}^{\omega }}{\mathit{dx}}+{\gamma }_n^2\xi \left(x\right)G_{h,n}^{\omega }=0$$E15

and

$${\phi }_{z_h}\frac{d^2{P}_{h,n}}{{\mathit{dz}}^2}+{\gamma }_n^2{u}_h{P}_{h,n}=0,$$E16

where ${\gamma }_n$ is a separation constant.

The first-order ordinary differential Eq. (15) has the solution

$$G_{h,n}^{\omega }\left(x\right)={\mu }_n\left(\omega \right).exp\left(-{\gamma }_n^2{\int }_0^x\xi \left(s\right)\mathit{ds}\right),$$

where ${\mu }_n$ is an arbitrary function depending on $\omega$.

Eq. (16) represents a Sturm-Liouville problem. Solutions of such problem form an eigenfunction basis of the form:

$$P_{h,n}\left(z\right)={\alpha }_{h,n}cos\left({\lambda }_{h,n}z\right)+{\beta }_{h,n}sin\left({\lambda }_{h,n}z\right),$$E17

where, ${\lambda }_{h,n}={\gamma }_n\sqrt{\frac{u_h}{{\phi }_{z_h}}}$.

The Eq. (17) satisfies the following boundary conditions:

$$\left\{\begin{array}{l}{\phi }_{z_1}\frac{{\mathit{dP}}_{1,n}\left(z_0\right)}{\mathit{dz}}=0\left(\mathrm{a}\right)\\ \begin{array}{l}\left\{\begin{array}{l}{P}_{h-1,n}\left(z_{h-1}\right)=P_{h,n}\left(z_{h-1}\right)\\ {\phi }_{z_{h-1}}\frac{{\mathit{dP}}_{h-1,n}\left(z_{h-1}\right)}{\mathit{dz}}={\phi }_{z_h}\frac{{\mathit{dP}}_{h,n}\left(z_{h-1}\right)}{\mathit{dz}},h\in \left\{2,\cdots ,H\right\}\left(\mathrm{b}\right)\end{array}\right.\\ {\phi }_{z_H}\frac{{\mathit{dP}}_{H,n}\left(z_H\right)}{\mathit{dz}}=0\left(\mathrm{c}\right)\end{array}\end{array}\right.$$E18

To calculate the expression of $P_{h,n}$, it comes down to calculate the values of ${\alpha }_{h,n}$ and ${\beta }_{h,n}$, on each of the sub-layer $\left[z_{h-1},z_h\right],h\in \left\{1,\cdots ,H\right\}$.

By solving the recursive system resulting from substitution of Eq. (17) in Eq. (18), we obtain respectively the formulations of ${\alpha }_{h,n}$ and ${\beta }_{h,n}$. More specifically,

For the first sub-layer, ${\alpha }_{1,n}$ and ${\beta }_{1,n}$ satisfy the equation:

$${\alpha }_{1,n}sin\left({\lambda }_{1,n}{z}_0\right)-{\beta }_{1,n}cos\left({\lambda }_{1,n}{z}_0\right)=0$$E19

from which, we can take ${\beta }_{1,n}=sin\left({\lambda }_{1,n}{z}_0\right)$, so that ${\alpha }_{1,n}=cos\left({\lambda }_{1,n}{z}_0\right)$, which means:

$$P_{1,n}\left(z\right)=cos\left({\lambda }_{1,n}\left(z-z_0\right)\right).$$

For the last sub-layer ($H^{\mathit{th}}$ sub-layer):

$${\alpha }_{H,n}=cot\left({\lambda }_{H,n}{z}_H\right){\beta }_{H,n}.$$E20

And for the intermediate sub-layers:

$$\begin{array}{ll}{\alpha }_{h,n}=& \frac{{\phi }_{z_h}{\lambda }_{h,n}-{\phi }_{z_{h-1}}{\lambda }_{h-1,n}}{2{\phi }_{z_h}{\lambda }_{h,n}}\left[\left(cos\left(\left({\lambda }_{h,n}+{\lambda }_{h-1,n}\right)z_{h-1}\right)+\frac{{\phi }_{z_h}{\lambda }_{h,n}+{\phi }_{z_{h-1}}{\lambda }_{h-1,n}}{{\phi }_{z_h}{\lambda }_{h,n}-{\phi }_{z_{h-1}}{\lambda }_{h-1,n}}.\right.\right.\\ & cos\left(\left({\lambda }_{h,n}-{\lambda }_{h-1,n}\right)z_{h-1}\right)){\alpha }_{h-1,n}+(sin\left(\left({\lambda }_{h,n}+{\lambda }_{h-1,n}\right)z_{h-1}\right)-\\ & \left.\left.\frac{{\phi }_{z_h}{\lambda }_{h,n}+{\phi }_{z_{h-1}}{\lambda }_{h-1,n}}{{\phi }_{z_h}{\lambda }_{h,n}-{\phi }_{z_{h-1}}{\lambda }_{h-1,n}}sin\left(\left({\lambda }_{h,n}-{\lambda }_{h-1,n}\right)z_{h-1}\right)\right){\beta }_{h-1,n}\right],\end{array}$$E21

and

$$\begin{array}{ll}{\beta }_{h,n}=& \frac{{\phi }_{z_h}{\lambda }_{h,n}-{\phi }_{z_{h-1}}{\lambda }_{h-1,n}}{2{\phi }_{z_h}{\lambda }_{h,n}}\left[\left(sin\left(\left({\lambda }_{h,n}+{\lambda }_{h-1,n}\right)z_{h-1}\right)+\frac{{\phi }_{z_h}{\lambda }_{h,n}+{\phi }_{z_{h-1}}{\lambda }_{h-1,n}}{{\phi }_{z_h}{\lambda }_{h,n}-{\phi }_{z_{h-1}}{\lambda }_{h-1,n}}.\right.\right.\\ & sin\left(\left({\lambda }_{h,n}-{\lambda }_{h-1,n}\right)z_{h-1}\right)){\alpha }_{h-1,n}-(cos\left(\left({\lambda }_{h,n}+{\lambda }_{h-1,n}\right)z_{h-1}\right)-\\ & \left.\left.\frac{{\phi }_{z_h}{\lambda }_{h,n}+{\phi }_{z_{h-1}}{\lambda }_{h-1,n}}{{\phi }_{z_h}{\lambda }_{h,n}-{\phi }_{z_{h-1}}{\lambda }_{h-1,n}}cos\left(\left({\lambda }_{h,n}-{\lambda }_{h-1,n}\right)z_{h-1}\right)\right){\beta }_{h-1,n}\right].\end{array}$$E22

The eigenvalues ${\gamma }_n,n\in {\mathrm{\mathbb{N}}}^{\ast }$ of this problem are real and discrete and the eigenfunctions are mutually orthogonal. The orthogonality relation developed by [16] for this class of (self-adjoint) problems with respect to the density $u_h$ on each interval $\left[z_{h-1},z_h\right],h\in \left\{1,\cdots ,H\right\}$ leads to

$$\sum _{h=1}^H{\int }_{z_{h-1}}^{z_h}{u}_h{P}_{h,m}\left(s\right)P_{h,n}\left(s\right)\mathit{ds}=\parallel P_{h,m}\parallel .\parallel P_{h,n}\parallel .{\delta }_{m,n},$$

where ${\delta }_{m,n}$ is the Kronecker symbol. Then, $\parallel P_{h,n}{\parallel }^2$ is written as follows.

$$\begin{array}{rcl}{‖P_{H,n}‖}^2& =& \sum _{h=1}^H{\int }_{z_{h-1}}^{z_h}{u}_h{\left(P_{h,n}\left(s\right)\right)}^2\mathit{ds}\\ & =& \{\begin{array}{l}\sum _{h=1}^H{u}_h\left(z_h-z_{h-1}\right),\text{when}{\gamma }_n=0\\ \sum _{h=1}^H\frac{u_h}{2{\lambda }_{h,n}}[sin\left({\lambda }_{h,n}\left(z_h-z_{h-1}\right)\right)(\left({\alpha }_{h,n}^2-{\beta }_{h,n}^2\right)cos\left({\lambda }_{h,n}\left(z_h+z_{h-1}\right)\right)+\\ 2{\alpha }_{h,n}{\beta }_{h,n}sin\left({\lambda }_{h,n}\left(z_h+z_{h-1}\right)\right))+{\lambda }_{h,n}\left({\alpha }_{h,n}^2+{\beta }_{h,n}^2\right)\left(z_h-z_{h-1}\right)],\text{when}{\gamma }_n\ne 0.\end{array}\end{array}$$

The eigenvalues of each sub-layer can be obtained by integrating Eq. (17) on each of the intervals $\left[z_{h-1},z_h\right],h\in \left\{1,\cdots ,H\right\}$ taking into account the boundary conditions Eqs. (19)–(21). Whereas the eigenfunctions $P_{H,n}$ form a complete set, so ${\chi }_H^{\omega }$ can be developed as:

$${\chi }_H^{\omega }\left(x,z\right)=\sum _{n=0}^{\infty }{a}_{n}exp\left(-{\gamma }_n^2{\int }_0^x\xi \left(s\right)\mathit{ds}\right)P_{H,n}\left(z\right).$$

Then, the coefficients $a_n$ are given by:

$$\begin{array}{ll}{a}_n& =\frac{1}{\parallel P_{h,n}{\parallel }^2}\sum _{\mathrm{\ell }=1}^H{\int }_{z_{\mathrm{\ell }-1}}^{z_{\mathrm{\ell }}}{u}_{\mathrm{\ell }}{\chi }_{\mathrm{\ell }}^{\omega }\left(0,z\right)P_{\mathrm{\ell },n}\left(z\right)\mathit{dz}\\ & =\left\{\begin{array}{l}\frac{Q}{{\sum }_{h=1}^H{u}_h\left(z_h-z_{h-1}\right)},\text{when}{\gamma }_n=0\\ \frac{Q}{\parallel P_{h,n}{\parallel }^2}{\sum }_{\mathrm{\ell }=1}^H{P}_{\mathrm{\ell },n}\left(H_s\right),\text{when}{\gamma }_n\ne 0\end{array}\right.\end{array}$$

The inverse Fourier transform of $exp\left({\left(2\pi \right)}^2{\omega }^2{\int }_0^x{\zeta }_y\left(s\right)\mathit{ds}\right)$ is:

$$\frac{1}{2\sqrt{\pi {\int }_0^x{\zeta }_y\left(s\right)\mathit{ds}}}exp\left(-\frac{y^2}{4{\int }_0^x{\zeta }_y\left(s\right)\mathit{ds}}\right)$$

Consequently,

$$\begin{array}{l}{c}_H\left(x,y,z\right)=\frac{Q}{2\sqrt{\pi {\int }_0^x{\zeta }_y\left(s\right)\mathit{ds}}}exp\left(-\frac{y^2}{4{\int }_0^x{\zeta }_y\left(s\right)\mathit{ds}}\right).\\ \sum _{n=0}^{\infty }exp\left(-{\gamma }_n^2{\int }_0^x\xi \left(s\right)\mathit{ds}\right)\frac{P_{H,n}\left(z\right)}{{‖P_{H,n}‖}^2}\left(\sum _{h=1}^H{P}_{h,n}\left(H_s\right)\right).\end{array}$$E23

The Crosswind-Integrated Concentration (CIC) model is a valuable tool for analyzing the behavior of pollutants in the atmosphere. It takes into account the horizontal wind component, which is a significant factor in pollutant transport. Mathematically, the crosswind-integrated concentration $c_H^y\left(x,z\right)$ is obtained by integrating $c_H\left(x,y,z\right)$ in Eq. (20) to y from $-\infty$ to $+\infty$:

$$c_H^y\left(x,z\right)={\int }_{-\infty }^{+\infty }{c}_H\left(x,y,z\right)\mathit{dy}$$

The (CIC) model considers the dispersion of pollutants in the crosswind direction due to the horizontal wind component and produces a two-dimensional concentration map that can be used to assess the impact of pollutants in both the downwind (x-axis) and crosswind (z-axis) directions by integrating the concentration of pollutants over the crosswind direction (y), it gives

$$c_H^y\left(x,z\right)=Q\sum _{n=0}^{\infty }exp\left(-{\gamma }_n^2{\int }_0^x\xi \left(s\right)\mathit{ds}\right)\frac{P_{H,n}\left(z\right)}{\parallel P_{H,n}{\parallel }^2}\sum _{h=1}^H{P}_{h,n}\left(H_s\right).$$E24

4. Parameterizing pollutant dispersion: Exploring model inputs

The wind speed profile u(z) is determined using the Deaves and Harris model [9], which is based on experimental measurements of wind speed profiles. This model is advantageous because it provides a more accurate representation of the logarithmic law at low heights, and improves the accuracy of both the logarithmic and power law models at moderate heights. The profile is given by the equation

$$u\left(z\right)=\frac{u_{\ast }}{k}\left[ln\left(\frac{z+z_0}{z_0}\right)+5.75\left(\frac{z}{H_{\mathit{mix}}}\right)-1.88{\left(\frac{z}{H_{\mathit{mix}}}\right)}^2-1.33{\left(\frac{z}{H_{\mathit{mix}}}\right)}^3+0.25{\left(\frac{z}{H_{\mathit{mix}}}\right)}^4\right],$$E25

where $u_{\ast }$ is the friction velocity, k is the von Karman constant, $z_0$ is the roughness length, and $H_{\mathit{mix}}$ is the height of the planetary boundary layer (PBL). The Deaves and Harris model includes terms that account for the effects of turbulence and thermal stratification on the wind profile. By using this model, we can better understand and predict the behavior of the wind at different heights within the PBL.

The modified form of the vertical eddy diffusivity coefficient $K_z$, given in Eq.(8), is adopted from [11], where the functional ${\phi }_z\left(z\right)$ is expressed as:

$${\phi }_z\left(z\right)=0.22H_{\mathit{mix}}{w}_{\ast }{\left(\frac{z}{H_{\mathit{mix}}}\left(1-\frac{z}{H_{\mathit{mix}}}\right)\right)}^{1/3}\left(1-exp\left(-\frac{4z}{H_{\mathit{mix}}}\right)-0.0003exp\left(\frac{8z}{H_{\mathit{mix}}}\right)\right),$$E26

where $w_{\ast }$ is the convective velocity. The integrable correction dimensionless function in Eq. (8) is defined in terms of the along-wind length scale $L_1$ as follows [10]:

$$\xi \left(x\right)=1-exp\left(-\frac{x}{L_1}\right).$$E27

The length $L_1$ is given in terms of u, ${\phi }_z$ and ${\sigma }_w$ as [10]:

$$L_1=\frac{1}{{\sigma }_w^2\left(H_s\right)}u\left(H_s\right){\phi }_z\left(H_s\right),$$

where ${\sigma }_w$ is the vertical turbulent intensity. Note that, there exist many expressions of ${\sigma }_w$, we adopt here the expression given by [17]:

$${\sigma }_w=0.96w_{\ast }{\left(\frac{3z}{H_{\mathit{mix}}}+\frac{\mid L\mid }{H_{\mathit{mix}}}\right)}^{1/3},$$

The Monin-Obukhov length (L) is a commonly parameter used in atmospheric dispersion modeling [18]. It plays a crucial role in characterizing the stability of the atmospheric boundary layer, which is essential for understanding and predicting the dispersion of pollutants or particles in the atmosphere. Specifically, the sign of L indicates the atmospheric stability state. A negative value of L indicates an unstable atmosphere, where turbulence is generated due to buoyancy forces. A positive value of L indicates a stable atmosphere, where turbulence is suppressed due to stratification. A large absolute value of L (i.e., $\mid L\mid \gg 1$) indicates a neutral atmosphere, where turbulence is neither generated nor suppressed due to buoyancy effects [18].

It should be noted that when $L_1\to 0$, the term $exp\left(-\frac{x}{L_1}\right)$ becomes negligible in (27), and $K_z$ will dependent only on $z$ (i.e., $K_z\equiv {\phi }_z$ in Eq. (8)).

The lateral eddy diffusivity is given by [12]:

$$K_y\left(x,z\right)=\frac{1}{2}u\left(z\right)\frac{d{\sigma }_y^2\left(x\right)}{\mathit{dx}},$$E28

where ${\sigma }_y$ is the standard deviation in the crosswind direction (depends only on x).

By identifying Eqs. (7) and (28), we obtain

$${\zeta }_y\left(x\right)=\frac{1}{2}\frac{d{\sigma }_y^2\left(x\right)}{\mathit{dx}}$$

5. Validating analytical solutions for air pollution dispersion modeling: analyzing experimental data

This section offers a comprehensive analysis of the data generated by our model, which is founded on the Deaves and Harris wind velocity profile. The model is constructed using two different formulations of vertical turbulent diffusivity. The first formulation, with $L_1>0$, factors in both the downwind distance and the vertical height, while the second formulation, with $L_1\to 0$, considers only the vertical height. The performance of the model, as presented in Eq. (24), is assessed and confirmed through the use of datasets from the Copenhagen diffusion experiments in Denmark and the Prairie Grass experiments in the USA.

Atmospheric dispersion experiments were conducted in Copenhagen, Denmark, from June 27 to November 9, 1978, to study air pollutant transport in the lower atmosphere. The tracer gas used in these experiments was sulfur hexafluoride (SF6), which was released from a tower located at a height of 115 meters ($H_s$), with a roughness length of 0.6 meters ($z_0$). Additional details of the parameters used in these experiments can be found in the publications of Gryning and Lyck [19] and Gryning et al. [20].

The Prairie Grass experiments, carried out near O’Neil, Nebraska, between March and August 1956, were a groundbreaking series of studies on atmospheric pollution behavior and transport. These experiments measured the dispersion of sulfur dioxide (SO₂) released from a 1.5 meter tower at five downwind distances from the source: 50, 100, 200, 400, and 800 meters. The roughness length of the area was 0.006 meters, which is critical in determining the spread of air pollutants. The experiments were extensively documented by Barad et al. [21] and Nieuwstadt et al. [22], providing essential initial insights into atmospheric dispersion processes.

Overall, the Copenhagen and Prairie Grass experiments provided essential knowledge on the transport and dispersion of air pollutants in the lower atmosphere. The Copenhagen experiments used an inert tracer gas to examine general dispersion patterns, while the Prairie Grass experiments focused on the downwind transport of SO₂, providing dispersion data across a range of distances. Both experiments relied on the roughness length of the terrain, a crucial parameter in determining how pollutants spread near the ground. The findings of these landmark experiments made significant contributions to scientists’ understanding of atmospheric pollution behavior, particularly in the PBL near the Earth’s surface.

It should be noted that the computations were carried out using two distinct programming environments, namely MATLAB and Wolfram Mathematica. Specifically, the function fzero, which utilizes a combination of the bisection, secant, and inverse quadratic interpolation methods, was defined and executed in the MATLAB environment. Conversely, the FindRoot function, which combines Brent’s and Newton’s methods along with a method that approximates the Jacobian, was defined and executed in the Wolfram Mathematica environment.

To implement the analytical solution in a practical sense, we have employed a discretization technique for the atmospheric boundary layer. Specifically, we have discretized the layer into two and four sublayers, respectively. In the case of two sub-layers, the boundary layer height is discretized using the following approach: For the two sub-layers case, the boundary layer height is divided into two sub-layers, and the discretization is taken as follows: ${\mathit{dz}}_1=\frac{7}{13}\left(H_{\mathit{mix}}-z_0\right)$ and ${\mathit{dz}}_2=\frac{6}{13}\left(H_{\mathit{mix}}-z_0\right)$. On the other hand, for the four sub-layers case, the boundary layer height is divided into four sub-layers, and the discretization is taken as follows: ${\mathit{dz}}_1=\frac{3}{15}\left(H_{\mathit{mix}}-z_0\right)$, ${\mathit{dz}}_2=\frac{3}{15}\left(H_{\mathit{mix}}-z_0\right)$, ${\mathit{dz}}_3=\frac{7}{15}\left(H_{\mathit{mix}}-z_0\right)$, and ${\mathit{dz}}_4=\frac{2}{15}\left(H_{\mathit{mix}}-z_0\right)$. These formulas are used to accurately model the behavior of the boundary layer height in numerical simulations.

In this work, the results of the numerical convergence for concentration predictions of atmospheric dispersion models are indicated as a function of the number of eigenvalues used in the model for different downwind distances from the source. Reference is made to experiments conducted in Copenhagen (experiments 1 to 3) and Prairie Grass (experiments 1, 5, and 7).

For the Copenhagen experiments, the models used four sub-layers in the vertical and tested two different formulations for calculating the vertical eddy diffusivity, a parameter that represents turbulence and mixing in the vertical direction. Figure 1 shows the numerical convergence of normalized crosswind-integrated concentration as a function of the number of eigenvalues used in the model. The convergence is shown for four different downwind distances (200, 500, 1000, and 1500 m) from the source. The sentence indicates that for Prairie Grass, only one formulation of vertical diffusivity was used because the mixing height ($H_m\mathit{ix}$) was much greater than the surface layer height ($H_s$), so the two formulations gave the same result.

Figure 1.

Copenhagen experiments number 1, 2, and 3: Numerical convergence of the proposed solution for two vertical turbulence formulations as a function of the number of eigenvalues at four distances using four sub-layers.

Figure 2 shows the numerical convergence of ground-level crosswind-integrated concentration for Prairie Grass experiments as a function of the number of eigenvalues at four downwind distances (50, 100, 200, and 800 m).

Figure 2.

Prairie grass experiment numbers 1, 5, and 7: Numerical convergence of the proposed solution as a function of the number of eigenvalues at four distances using four sub-layers.

Numerical simulations were conducted to investigate the influence of the vertical eddy diffusivity formulation on the normalized crosswind-integrated concentration of air pollution near a point source. Specifically, the simulations were performed for Copenhagen experiments numbered 1, 2, and 3, and the results are presented in Figure 3.

Figure 3.

Normalized concentration Eq. (24) as a function of distance using four sub-layers for Copenhagen experiment numbers 1, 2, and 3.

The figure depicts the normalized crosswind-integrated concentration $\frac{c_H^y\left(x,z_0\right)}{Q}$ at ground level as a function of downwind distance, using two different vertical turbulence formulations. The results indicate that the choice of vertical eddy diffusivity formulation has a significant impact on the concentration near the point source. Notably, the figure shows that the concentration is higher when using the first formulation of the vertical eddy diffusivity, compared to the second one. Furthermore, the figure illustrates the expected decrease in crosswind-integrated concentration with increasing downwind distance due to the dispersion of air pollutants caused by turbulent mixing in the atmosphere. In conclusion, these simulations provide valuable insights into air pollution dispersion behavior and can inform future experiments and modeling efforts.

The quality and performance of models are typically evaluated by comparing their predictions to observed (or measured) values. When observations are available, a common approach to evaluating model performance is to plot the predicted values against the observed values in a scatter diagram.

In this study, the performance of the models for the Copenhagen and Prairie Grass experiments was evaluated using scatter diagrams, which are shown in Figures 4 and 5, respectively. Figure 4 compares the predicted and observed crosswind-integrated concentrations for two different formulations of the vertical diffusivity parameter and two different numbers of sub-layers used to represent the atmospheric boundary layer (with $H=2$ and $H=4$), while Figure 5 compares the predicted and observed ground-level crosswind integrated concentrations for two different numbers of sub-layers in the atmospheric boundary layer (again with $H=2$ and $H=4$).

Figure 4.

Comparison of observed and predicted crosswind-integrated concentrations in Copenhagen (top panel) and prairie grass (bottom panel) experiments using scatter plots. The one-to-one line ($y=x$) and factor-of-two lines ($y=0.5\times x$ and $y=2\times x$) are shown.

Figure 5.

Isolines of non-dimensional predicted crosswind-integrated concentration as a function of non-dimensional distance $\frac{w_{\ast }x}{u_{\mathit{\text{mean}}}{H}_{\mathit{mix}}}$ and non-dimensional depth $\frac{z}{H_{\mathit{mix}}}$ are presented for Copenhagen experiment number 2 (top panel) and prairie grass experiment number 7 (bottom panel) at different source heights: (a) $\frac{H_s}{H_{\mathit{mix}}}=0.1$, (b) $\frac{H_s}{H_{\mathit{mix}}}=0.25$, (c) $\frac{H_s}{H_{\mathit{mix}}}=0.5$, and (d) $\frac{H_s}{H_{\mathit{mix}}}=0.75$, using four sub-layers.

The scatter diagrams show good agreement between the predicted and observed values, indicating that the models are well-parameterized. However, a visual inspection of the figures also reveals that the observed concentrations tended to be slightly higher than the predicted concentrations.

While scatter diagrams provide a useful visual evaluation, a more rigorous evaluation of model performance is achieved through statistical metrics. In this study, the statistical indices used to evaluate the models’ performance are defined in the reference [23]. These statistical metrics provide a more detailed and quantitative assessment of the models’ performance, which can be used to identify any areas of weakness or sources of error in the models. We defined the statistical indices used in our study as:

• Normalized Mean Square Error (NMSE): A measure of the accuracy of a model or prediction, calculated as the mean square error divided by the variance of the observed data:

  $$\mathit{\text{NMSE}}=\frac{\overline{{\left(c_o-c_p\right)}^2}}{\overline{c_o}\times \overline{c_p}},$$

• Mean Relative Square Error (MRSE): A measure of the accuracy of a model or prediction, calculated as the mean square error divided by the mean of the observed data squared:

  $$\mathit{\text{MRSE}}=4\overline{{\left(\left(c_o-c_p\right)/\overline{\left(c_o+c_p\right)}\right)}^2},$$

• Correlation Coefficient (COR): A measure of the linear relationship between two variables, usually the predicted and observed values of a model or prediction. It ranges from −1 to 1, with 1 indicating a perfect positive correlation and − 1 indicating a perfect negative correlation:

  $$\mathit{COR}=\frac{\overline{\left(c_p-\overline{c_p}\right)\left(c_o-\overline{c_o}\right)}}{\left({\sigma }_o{\sigma }_p\right)},$$

• Fractional Bias (FB): A measure of the systematic bias in a model or prediction, calculated as the difference between the mean predicted and observed values divided by the mean observed value:

  $$\mathit{FB}=2\frac{\left(\overline{c_o}-\overline{c_p}\right)}{\left(\overline{c_o}+\overline{c_p}\right)},$$

• Fractional Standard Deviation (FS): A measure of the variability in a model or prediction, calculated as the standard deviation of the predicted values divided by the mean observed value:

  $$\mathit{FS}=2\frac{\left({\sigma }_o-{\sigma }_p\right)}{\left({\sigma }_o+{\sigma }_p\right)},$$

• Mean Geometric Ratio (MG): A measure of the logarithmic bias in a model or prediction, calculated as the exponential of the mean of the differences between the logarithms of the predicted and observed values:

  $$\mathit{MG}=exp\left(\overline{log\left(c_o\right)}-\overline{log\left(c_p\right)}\right),$$

• Variance of Geometric Ratio (VG): A measure of the variability in the logarithmic bias of a model or prediction, calculated as the variance of the differences between the logarithms of the predicted and observed values:

  $$\mathit{VG}=exp\left[\overline{{\left(ln\left(c_o\right)-ln\left(c_p\right)\right)}^2}\right],$$

• Factor of Two (FAC2): A measure of the accuracy of a numerical weather prediction model, defined as the distance between the observed and predicted position of a weather system divided by the radius of the weather system. If the distance is less than or equal to twice the radius, the forecast is considered accurate and has a FAC2 score of 1. If the distance is greater than twice the radius, the forecast is considered inaccurate and has a FAC2 score of less than 1:

$$0.5\le \left(c_p/c_o\right)\le 2.$$

The satisfactory numerical results presented in Figure 4 are supported by the values obtained from the statistical performance measures. Upon examining these values, it is observed that they mostly fall within the range of acceptable performance for the two proposed formulations for vertical turbulent diffusivity. However, it should be noted that the simulated concentrations for experiments conducted with Prairie Grass tend to slightly overestimate the measured quantities. Specifically, when considering the turbulent diffusivity that depends on both x and z, the model produces results that are relatively more promising than those obtained with other methods.

Figure 5 depicts the spatial distribution of pollutants for two different experiments—Copenhagen experiment number 2 and Prairie Grass experiment number 7, using contour lines that represent the adimensional transverse integrated concentration predicted based on the adimensional distance for experiment 2 and adimensional depth for experiment 7, where the adimensional distance for experiment 2 is defined as $w_{\ast }x/\left(u_{\text{mean}}{H}_{\mathit{mix}}\right)$.

6. Conclusions

A solution has been developed for the 3-dimensional stationary state atmospheric diffusion equation, which incorporates more realistic formulations of wind velocity profiles and two formulations of vertical turbulent diffusivity. This solution was achieved through a combination of an appropriate auxiliary eigenvalue problem with mathematical induction, while a transcendental equation was developed to determine the eigenvalues for any number of subdomains. The advection-diffusion equation was solved by dividing the planetary boundary layer into discrete layers and assuming average values of turbulent diffusivity and wind velocity in each subdomain. The solution was evaluated using the Copenhagen and Prairie Grass experiments for two and four layers, and the numerical convergence of the solution was verified based on the number of eigenvalues. The results indicated good agreement between predicted and observed values, and most of the calculated statistical indices were within an acceptable range for model performance. This current model could be a promising approach for accurately predicting atmospheric pollutant dispersion and may also be applicable to other continuous flows (Table 1).

	Copenhagen experiments
Models	NMSE	MRSE	COR	FB	FS	MG	VG	FAC2
$K_z=\xi \left(x\right){\phi }_z\left(z\right)$$\left(H=2\right)$	0.0533	0.0533	0.8657	−0.0135	0.0018	0.9660	1.0885	1.096
$K_z={\phi }_z\left(z\right)$$\left(H=2\right)$	0.0660	0.0660	0.8498	0.0432	0.03897	1.0235	1.1007	1.055
$K_z=\xi \left(x\right){\phi }_z\left(z\right)$$\left(H=4\right)$	0.0768	0.0768	0.8308	0.0543	0.1686	1.0273	1.1022	0.9235
$K_z={\phi }_z\left(z\right)$$\left(H=4\right)$	0.0892	0.0890	0.8274	0.0973	0.2065	1.0726	1.1124	0.8882
[24] ($H=2$)	0.08	—	0.86	−0.02	0.05	—	—	1.0
[24] ($H=4$)	0.1	—	0.82	0.1	0.04	—	—	0.92
[25]	0.069	—	—	−0.009	0.051	0.996	1.055	1.009
	Prairie Grass experiments
$K_z=\xi \left(x\right){\phi }_z\left(z\right)$$\left(H=2\right)$	0.0210	0.0209	0.9809	−0.0503	−0.0108	0.9018	1.0645	0.9807
$K_z=\xi \left(x\right){\phi }_z\left(z\right)$$\left(H=4\right)$	0.0613	0.0608	0.9763	−0.1809	– 0.0514	0.8616	1.3228	0.9373
[26]	0.25	—	0.92	0.03	0.20	—	—	0.68
[27]	0.04	—	0.96	−0.09	0.13	—	—	0.79

Table 1.

Statistical measures for the Copenhagen and Prairie grass experiments.