Generalized 3D Model of Crosswind Concentrations and Deposition in the Atmospheric Boundary Layer

1. Introduction

Understanding the dispersion, transport, and fate of atmospheric pollutants necessitates a thorough analysis of deposition processes. These mechanisms are essential in determining plume behavior and concentration patterns within the atmospheric boundary layer (ABL) [29]. However, Gaussian plume model estimates often fall short due to particle settling, deposition, or vertical gradients in diffusivity and wind speed [27]. The Earth’s surface acts as the primary reservoir for trace gases and particles, which transfer through wet and dry deposition pathways. Wet deposition is driven by precipitation, whereas dry deposition involves the adhesion or chemical reaction of trace gases and particles on surfaces like vegetation, soil, water bodies, and structures. Factors such as atmospheric turbulence intensity, the properties of gases and particles, and surface characteristics influence the efficiency of dry deposition [24]. In the absence of turbulence, even surfaces with high adsorptive capacity are ineffective in removing suspended materials, as the process is limited by molecular diffusion. Gravitational settling also plays a role in the dry deposition of particles, further complicating these dynamics [10]. The role of deposition processes in the dispersion of passive contaminants within the ABL is conventionally modeled using the advection-diffusion equation, as dictated by gradient transport theory, with deposition fluxes incorporated as boundary conditions at the surface [9]. To address this, a modified Gaussian model was formulated, initially presuming constant wind speed and homogeneous eddy diffusivity [13]. For enhanced applicability, eddy diffusivities were subsequently adjusted based on dispersion parameters relative to the downwind distance, and wind speed was modeled using a power-law relationship with vertical height [17]. However, this method introduces mathematical inconsistencies due to the employment of unrealistic parameterizations [26].

In this paper, we present an innovative approach aimed at improving the understanding of pollutant behavior in the atmospheric boundary layer (ABL). Our research introduces a generalized, one-of-a-kind model capable of adapting to all existing parameterization schemes. This model holistically integrates various wind speed profiles, turbulent diffusion coefficients, and dry deposition criteria, thereby offering unmatched flexibility and accuracy in the study of pollutant dispersion dynamics. Our model subdivides the atmospheric boundary layer (ABL) into distinct sub-layers, facilitating a granular examination of pollutant dispersion dynamics. Each sub-layer is parameterized with specific wind speed profiles and turbulent diffusivity coefficients, adopting their average values. This segmented approach allows for a detailed analysis of vertical and horizontal variations in atmospheric properties, thereby providing a more precise understanding of the mechanisms of pollutant transport and deposition. The key contributions of this work include the development of a generalized analytical model, the incorporation of closed-form solutions, and the meticulous validation of these models through empirical data. By addressing the intricacies of turbulent dispersion and dry deposition, our research offers a comprehensive framework that enhances the predictability and understanding of atmospheric pollutant behavior. In the course of this contribution, we have demonstrated the equivalence between two formulations of the pollutant dispersion problem in the atmospheric boundary layer. The problems ($\mathcal{PS}$) and ($\mathcal{P}$) differ in their treatment of the source term: in ($\mathcal{PS}$), the source term $Q·\delta \left(x\right)\otimes \delta \left(y\right)\otimes \delta (z-H_s)$ is directly integrated into the main differential equation, representing a point source, whereas in ($\mathcal{P}$), this term is included in the boundary conditions. The equivalence between ($\mathcal{PS}$) and ($\mathcal{P}$) is proven using mollifiers (${\rho }_n$), smooth functions that handle the singularity of the source term [11]. The solution obtained not only generalizes but also extends previously established solutions. For a single layer and two-dimensional case, it corresponds to the solutions developed by Kumar [15] and Seinfeld [25]. In the context of multi-layer and three-dimensional scenarios, our model aligns with the solutions provided by Farhane et al. [7] and Marie et al. [19] This demonstrates the robustness and validity of our model across various parameterization conditions, confirming its applicability and accuracy in diverse atmospheric settings. The convergence of the proposed solution is demonstrated through rigorous analysis, showing that a minimal number of eigenvalues are sufficient. This efficiency offers a significant advantage over methods like the Generalized Integral Laplace Transform Technique (GILTT) [28], which often require many more eigenvalues for accuracy. Achieving convergence with fewer eigenvalues ensures both accuracy and computational efficiency, enhancing the model’s applicability in real-time atmospheric studies.

Validation against the Hanford experiment data [3] demonstrates the model’s accuracy, showing a strong correlation and minimal bias between observed and projected pollutant concentrations. This confirms the model’s reliability and statistical robustness. Compared with existing models by Marie et al.[19] and Laaouaoucha et al.[16] for cases without deposition, and Moreira et al.[21] and Kumar and Sharan [15] for cases with deposition, our model exhibits superior performance in predicting pollutant dispersion. These findings position our model as a powerful tool for investigating complex environmental phenomena, offering critical insights for both scientific research and practical applications in air quality management. To accomplish the goals of this study, the paper is organized as follows: Section 2 covers the modeling of the advection-diffusion equation. Section 3 explains the construction of the solution and the employed parameterizations. Section 4 showcases the results and provides a statistical performance analysis. Finally, Section 5 concludes the study.

2. Atmospheric Pollutant Advection-Diffusion Equation Modeling and Construction of Analytical Solutions

The advection-diffusion equation related to atmospheric pollution is fundamentally expressed as a principle of conservation of suspended particles. This statement is mathematically formalized by the following expression [31]:

$${\displaystyle \frac{D{c}_{tot}}{Dt}}={\displaystyle \frac{\partial c_{tot}}{\partial t}}+\nabla ·\left({\mathbf{U}}_{tot}{c}_{tot}\right)=S$$

(1)

where:

• ${\displaystyle \frac{D}{Dt}}$: Represents the substantial derivative, capturing the rate of change in pollutant concentration from the perspective of a moving air parcel, embodying both the temporal and spatial variations due to airflow.
• $C_{tot}$: Stands for the total concentration of pollutants, which includes both the average concentration and the transient fluctuations attributable to atmospheric turbulence.
• $U_{tot}$: Denotes the total wind velocity vector represented by ${(U,V,W)}^T$, where U indicates the east-west component, V the north-south component, and W the vertical component, encompassing both the average wind speed and its instantaneous fluctuations.
• $\partial /\partial t$: Indicates the partial derivative with respect to time, denoting how the concentration of pollutants changes at a fixed point as time progresses.
• ∇: Is the gradient operator, employed to calculate the spatial gradient of the pollutant concentration.
• $\nabla ·\left({\mathbf{U}}_{tot}{c}_{tot}\right)$: Describes the divergence of the wind velocity and pollutant concentration product, representing the net movement of pollutants carried by the wind across the spatial domain.
• S: Is the source term within the equation, quantifying the addition or removal of pollutants from the atmosphere via various processes such as industrial emissions, chemical reactions, or the deposition and uptake by terrestrial surfaces and vegetation.

Reynolds decomposition [30] is a pivotal technique for analyzing turbulent flows. It separates the total wind velocity ${\mathbf{U}}_{tot}$ and pollutant concentration $c_{tot}$ into their average components, $\mathbf{U}$ and c, and their fluctuating components, ${\mathbf{U}}^{\prime }$ and $c^{\prime }$, respectively. This separation is mathematically represented as

$${\mathbf{U}}_{tot}=\mathbf{U}+{\mathbf{U}}^{\prime }$$

and

$$c_{tot}=c+c^{\prime }$$

where the prime notation (’) indicates the fluctuating part of the quantity. By doing so, it facilitates the understanding of the complex interplay between the mean flow and the superimposed random fluctuations that characterize turbulent diffusion of contaminants in the atmosphere.

Incorporating the decomposed expressions for $c_{tot}$ and ${\mathbf{U}}_{tot}$, the advection-diffusion equation 1 is reformulated as:

$${\displaystyle \frac{\partial (c+c^{\prime })}{\partial t}}+\nabla ·\left((\mathbf{U}+{\mathbf{U}}^{\prime })(c+c^{\prime })\right)=S$$

Applying the averaging operator across the equation, the mean of the fluctuations alone cancels out, as they are, by definition, zero:

$$\overline{{\displaystyle \frac{\partial c^{\prime }}{\partial t}}}=0,\overline{{\mathbf{U}}^{\prime }}=\mathbf{0},\overline{c^{\prime }}=0$$

Expanding the term $\nabla ·\left(\mathbf{U}c\right)$, we get:

$$\nabla ·\left(\mathbf{U}c\right)=\nabla ·(\mathbf{U}c+\mathbf{U}{c}^{\prime }+{\mathbf{U}}^{\prime }c+{\mathbf{U}}^{\prime }{c}^{\prime })$$

After averaging, only the terms involving the averages and the term ${\mathbf{U}}^{\prime }{c}^{\prime }$ remain:

$$\overline{\nabla ·\left(\mathbf{U}c\right)}=\nabla ·\left(\mathbf{U}c\right)+\nabla ·\overline{{\mathbf{U}}^{\prime }{c}^{\prime }}$$

Assembling the mean terms and subtracting the turbulent diffusion term $\nabla ·\overline{{\mathbf{U}}^{\prime }{c}^{\prime }}$ from the advection term $\nabla ·\left(\mathbf{U}c\right)$, we ultimately derive the advection-diffusion equation for turbulent flow:

$${\displaystyle \frac{\partial c}{\partial t}}+\mathbf{U}·\nabla c=-\nabla ·\overline{{\mathbf{U}}^{\prime }{c}^{\prime }}+S$$

This equation now accurately represents the transport and spread of pollutants in a turbulent atmosphere, accounting for the mean motion of the air as well as the chaotic eddies that enhance diffusion.

In atmospheric modeling, the fluctuating wind component, denoted as ${\mathbf{U}}^{\prime }={(u^{\prime },v^{\prime },w^{\prime })}^T$, captures the erratic behavior of air movement. The gradient transport hypothesis, also known as K-theory, elucidates the mechanism by which turbulence drives materials against their concentration gradients. This is mathematically expressed by :

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

(2)

where the eddy diffusivity matrix $\mathbf{K}$, defined as $\mathrm{diag}(K_x,K_y,K_z)$, incorporates the eddy diffusivity coefficients $K_x$, $K_y$, and $K_z$ for the x-, y-, and z-directions, respectively. This setup enriches our understanding of turbulent diffusion and its impact on the spatial distribution of pollutants.

Furthermore, the equation $\overline{{\mathbf{U}}^{\prime }{c}^{\prime }}={\left(\overline{u^{\prime }{c}^{\prime }},\overline{v^{\prime }{c}^{\prime }},\overline{w^{\prime }{c}^{\prime }}\right)}^T$ accurately models turbulent fluxes, underscoring the model’s effectiveness in simulating atmospheric pollutant transport and distribution, and highlights the pivotal role of turbulent diffusion as described by Fick’s laws.:

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

Integrating Equation 2 with the mass continuity principle yields the advection-diffusion equation 1, as rigorously documented by [31]:

$${\partial }_{t}c+\mathbf{U}·\nabla c=\nabla ·(\mathbf{K}\nabla c)+S$$

(3)

In physics and engineering disciplines, the concept of point forces—actions highly localized in both space and time—is a recurrent theme. To mathematically model such phenomena, the Dirac delta distribution is utilized, enabling the precise identification of point sources. For instance, a point source emitting pollutants at a continuous rate Q at a specific location $(x_s,y_s,z_s)$ is represented by the source term $S(x_s,y_s,z_s)$, formulated as:

$$S(x_s,y_s,z_s)=Q·\delta (x-x_s)\otimes \delta (y-y_s)\otimes \delta (z-z_s),$$

with ⊗ symbolizing the tensor product.

To accurately model the transport of the emitted pollutants, the following boundary conditions are considered:

(i)

The apex of the mixing or inversion layer functions as an impervious boundary preventing the diffusion of contaminants. Consequently, a reflective flux boundary condition is postulated at the height $H_{mix}$ of the inversion layer, articulated mathematically as:

$$K_z(x,z){\displaystyle \frac{\partial c}{\partial z}}=0,\mathrm{at}z=H_{mix}.$$

(ii)

Addressing the scenario where a pollutant is partially absorbed by the ground surface, Calder [2] delineates an applicable boundary condition in flux terms. It is posited that the ground surface may facilitate partial removal or absorption of the pollutant via deposition, characterized by a deposition velocity $V_g$:

$$K_z(x,z){\displaystyle \frac{\partial c}{\partial z}}=V_{g}c,\mathrm{at}z=z_0.$$

where $z_0$ represents the surface roughness length.

(iii)

The lateral diffusive flux diminishes with increasing distance from the source, approaching zero as y extends towards the lateral boundary $L_y$, formalized by:

$$K_y(x,z){\displaystyle \frac{\partial c}{\partial y}}⟶0;\left|y\right|⟶L_y,$$

Within the context of this article, key assumptions include:

• The system reaches a state of equilibrium dispersion ($\partial c/\partial t=0$);
• Contributions from $v(\partial c/\partial y)$ and $w(\partial c/\partial z)$ are omitted, given that the x-axis aligns with the mean wind direction, rendering the w and v components of wind velocity minor;
• Turbulent diffusion along the mean wind direction is disregarded in comparison to advective transport, implying:

  $$u{\displaystyle \frac{\partial c}{\partial x}}\gg {\displaystyle \frac{\partial }{\partial x}}\left(K_x{\displaystyle \frac{\partial c}{\partial x}}\right).$$

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

(4)

Theorem 1.

The problem $\left(\mathcal{PS}\right)$ is equivalent to the problem $\left(\mathcal{PS}\right)$ defined by:

(5)

Proof of Theorem 1.

The steps for establishing equivalence are outlined as follows:

• Consider $c_{\mathcal{P}}(x,y,z)$ as the solution for $x>0$ of the system $\left(\mathcal{P}\right)$, and define $c_{\mathcal{PS}}(x,y,z)$ as follows:

  $$c_{\mathcal{PS}}(x,y,z)=\left\{\begin{array}{cc}{c}_{\mathcal{P}}(x,y,z)\hfill & \mathrm{for}x>0,\hfill \\ 0\hfill & \mathrm{for}x\le 0.\hfill \end{array}\right.$$
  We will establish that $c_{\mathcal{PS}}$ satisfies the requirements of system $\left(\mathcal{PS}\right)$.
• To improve solution regularity, we apply mollifiers ${\rho }_n$, as referenced in [11], over $]0,+\infty [$, defined as follows:

  $$c_{\mathcal{PS},n}(x,y,z)=c_{\mathcal{PS}}(x,y,z)\ast {\rho }_n\left(x\right)=\left\{\begin{array}{cc}{c}_{\mathcal{P},n}(x,y,z);\hfill & x>0,\hfill \\ 0;\hfill & \mathrm{otherwise},\hfill \end{array}\right.=c_{\mathcal{P},n}(x,y,z)·H\left(x\right)$$

  with $c_{\mathcal{P},n}(x,y,z)=c_{\mathcal{P}}(x,y,z)\ast {\rho }_n\left(x\right)$ and using the Heaviside function $H\left(x\right)$, which is 1 for $x>0$ and 0 otherwise, where * denotes convolution.
• We differentiate $c_{\mathcal{PS},n}$ with respect to x, leading to:

  $${\displaystyle \frac{\partial c_{\mathcal{PS},n}}{\partial x}}={\displaystyle \frac{\partial c_{\mathcal{P},n}}{\partial x}}·H\left(x\right)+c_{\mathcal{P},n}·H^{\prime }\left(x\right)={\displaystyle \frac{\partial c_{\mathcal{P},n}}{\partial x}}·H\left(x\right)+c_{\mathcal{P},n}(0,y,z)·\delta \left(x\right),$$
• Integrating the advection and diffusion terms, we find:

  $$\begin{array}{cc}\hfill u\left(z\right){\displaystyle \frac{\partial c_{\mathcal{PS},n}}{\partial x}}-\nabla ·(K\nabla c_{\mathcal{PS},n})& =u\left(z\right){\displaystyle \frac{\partial c_{\mathcal{P},n}}{\partial x}}-\nabla ·(K\nabla c_{\mathcal{P},n})+u\left(z\right)c_{\mathcal{P},n}(0,y,z)·\delta \left(x\right)\hfill \\ \hfill & =Q·\delta \left(x\right)\otimes \delta \left(y\right)\otimes \delta (z-H_s).\hfill \end{array}$$
  This confirms that $c_{\mathcal{P}}$ corresponds to the system $\left(\mathcal{P}\right)$ when n approaches infinity, thus completing the proof of equivalence between the systems $\left(\mathcal{PS}\right)$ and $\left(\mathcal{P}\right)$.

□

The problem $\left(\mathcal{P}\right)$ can be transformed into a series of specific advection-diffusion problems defined for each sub-layer of the atmospheric boundary layer, subdivided into ${\displaystyle H}$ levels. For each level, the formulations of the turbulent diffusivity coefficients and wind profile adopt average values, as cited in [6]. Furthermore, two boundary conditions are applied at the ABL’s extremities, at $z=z_0$ and at $H_{mix}$ (the height of the ABL), supplemented by continuity conditions for both the concentration and its flux at the interfaces between sub-layers. These continuity conditions, as shown in Figure 1, are expressed mathematically as follows

• The concentration at the end of a sub-layer must equal the concentration at the beginning of the subsequent sub-layer to ensure the continuity of concentration across sub-layer interfaces, as denoted by the equation

  $$c_h=c_{h+1},forh=1,2,\dots ,H-1$$
• The concentration flux, which is the product of diffusivity ($K_z$) and the concentration gradient with respect to the vertical axis (z), must also be continuous across the interfaces between sub-layers. This continuity is captured by the equation

  $$K_{z,h}{\displaystyle \frac{\partial c_h}{\partial z}}=K_{z,h+1}{\displaystyle \frac{\partial c_{h+1}}{\partial z}},forh=1,2,\dots ,H-1.$$

We consider the atmospheric boundary layer’s intricate dynamics, where the eddy diffusivities adopt separable formulations that illuminate the interplay between horizontal and vertical mixing processes [7]. Specifically, these formulations are expressed as:

$$\begin{array}{cc}\hfill K_y(x,z)& ={\zeta }_y\left(x\right)u_h,\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill K_z(x,z)& =\xi \left(x\right){\phi }_{z,h}.\hfill \end{array}$$

(7)

With this foundation, the equation $\left({\mathcal{P}}_H\right)$ outlines the advection-diffusion dynamics across each sub-layer of the atmospheric boundary layer, represented through a set of partial differential equations, expressed as follows:

$$\begin{array}{c}\hfill {\displaystyle \left({\mathcal{P}}_H\right):\left\{\begin{array}{c}{\displaystyle u_h{\displaystyle \frac{\partial c_h}{\partial x}}={\zeta }_y\left(x\right)u_h{\displaystyle \frac{{\partial }^2{c}_h}{\partial y^2}}+\xi \left(x\right){\phi }_{z_h}{\displaystyle \frac{{\partial }^2{c}_h}{\partial z^2}};}(x,y,z)\in ]0,L_x[\times ]-L_y,L_y[\times ]z_{h-1},z_h[,\hfill \\ \forall h\in \{1,\cdots ,H\}\hfill \\ Undertheconditions:\left\{\begin{array}{c}{\phi }_{z_1}\frac{\partial c_1(x,y,z_0)}{\partial z}={\displaystyle V_g{c}_1(x,y,z_0);}\hfill \\ \\ \left\{\begin{array}{cc}{c}_{h-1}(x,y,z_{h-1})=c_h(x,y,z_{h-1})\hfill & \\ \hfill & ,\forall h\in \{2,\cdots ,H\},\hfill \\ {\phi }_{z_{h-1}}\frac{\partial c_{h-1}(x,y,z_{h-1})}{\partial z}={\phi }_{z_h}\frac{\partial c_h(x,y,z_{h-1})}{\partial z}\hfill \end{array}\right.\hfill \\ \\ {\phi }_{z_H}\frac{\partial c_H(x,y,z_H)}{\partial z}=0\hfill \end{array}\right.\hfill \end{array}\right.}\end{array}$$

(8)

Theorem 2.

(Solution to the${\mathcal{P}}_H$Problem)

The concentration within the $H^{th}$ sub-layer, denoted by $c_H(x,y,z)$, which is the solution to the problem ${\mathcal{P}}_H$, is given by the following expression:

$$\begin{array}{cc}\hfill c_H(x,y,z)=& {\displaystyle \frac{Q}{2\sqrt{\pi {\int }_0^x{\zeta }_y\left(s\right)ds}}}exp\left(-{\displaystyle \frac{y^2}{4{\int }_0^x{\zeta }_y\left(s\right)ds}}\right)·[\sum _{n=1}^{\infty }exp\left(-{\gamma }_n^2{\int }_0^x\xi \left(s\right)ds\right)P_{H,n}\left(z\right)·\{\sum _{h=2}^H{P}_{h,n}\left(H_s\right)+\hfill \\ \hfill & \sqrt{2}{\displaystyle \frac{cos\left({\lambda }_{1,n}(H_s-z_0)\right)+v_{g}sin\left({\lambda }_{1,n}(H_s-z_0)\right)}{\sqrt{u_1\left(1+v_g^2\right)(z_1-z_0)}}}\}+{\displaystyle \frac{1}{\sum _{h=1}^H{u}_h(z_h-z_{h-1})}}]\hfill \end{array}$$

(9)

where

$$P_{h,n}=\left\{\begin{array}{cc}{\displaystyle \frac{1}{{\left[\sum _{h=1}^H{u}_h(z_h-z_{h-1})\right]}^{1/2}}};\hfill & \mathrm{if}{\gamma }_n=0\hfill \\ \hfill & \\ {\displaystyle \frac{{\alpha }_{h,n}cos\left({\lambda }_{h,n}z\right)+{\beta }_{h,n}sin\left({\lambda }_{h,n}z\right)}{\parallel P_{h,n}\parallel }};\hfill & \mathrm{if}{\gamma }_n\ne 0\hfill \end{array}\right.$$

(10)

with :

• The coefficient ${\alpha }_{h,n}$ is determined by the following relationship:

  $$\begin{array}{cc}\hfill {\alpha }_{h,n}=& {\displaystyle \frac{{\phi }_{z_h}{\lambda }_{h,n}-{\phi }_{z_{h-1}}{\lambda }_{h-1,n}}{2{\phi }_{z_h}{\lambda }_{h,n}}}\left[\right(cos\left(({\lambda }_{h,n}+{\lambda }_{h-1,n})z_{h-1}\right)+{\displaystyle \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}}}\hfill \\ \hfill & cos\left(({\lambda }_{h,n}-{\lambda }_{h-1,n})z_{h-1}\right)){\alpha }_{h-1,n}+(sin\left(({\lambda }_{h,n}+{\lambda }_{h-1,n})z_{h-1}\right)-\hfill \\ \hfill & {\displaystyle \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(({\lambda }_{h,n}-{\lambda }_{h-1,n})z_{h-1}\right)\left){\beta }_{h-1,n}\right];\hfill \end{array}$$

  (11)
• The coefficient ${\beta }_{h,n}$ is determined by the following relationship:

  $$\begin{array}{cc}\hfill {\beta }_{h,n}=& {\displaystyle \frac{{\phi }_{z_h}{\lambda }_{h,n}-{\phi }_{z_{h-1}}{\lambda }_{h-1,n}}{2{\phi }_{z_h}{\lambda }_{h,n}}}\left[\right(sin\left(({\lambda }_{h,n}+{\lambda }_{h-1,n})z_{h-1}\right)+{\displaystyle \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}}}\hfill \\ \hfill & sin\left(({\lambda }_{h,n}-{\lambda }_{h-1,n})z_{h-1}\right)){\alpha }_{h-1,n}-(cos\left(({\lambda }_{h,n}+{\lambda }_{h-1,n})z_{h-1}\right)-\hfill \\ \hfill & {\displaystyle \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(({\lambda }_{h,n}-{\lambda }_{h-1,n})z_{h-1}\right)\left){\beta }_{h-1,n}\right];\hfill \end{array}$$

  (12)
• The norm $\parallel P_{h,n}\parallel$ is defined as:

  $$\begin{array}{cc}\hfill \parallel P_{H,n}{\parallel }^2=& \sum _{h=1}^H{\int }_{z_{h-1}}^{z_h}{u}_h{\left(P_{h,n}\left(s\right)\right)}^{2}ds\hfill \\ \hfill =& \sum _{h=1}^H{\displaystyle \frac{u_h}{2{\lambda }_{h,n}}}[sin\left({\lambda }_{h,n}(z_h-z_{h-1})\right)(({\alpha }_{h,n}^2-{\beta }_{h,n}^2)cos\left({\lambda }_{h,n}(z_h+z_{h-1})\right)\hfill \\ \hfill & +2{\alpha }_{h,n}{\beta }_{h,n}sin\left({\lambda }_{h,n}(z_h+z_{h-1})\right))+{\lambda }_{h,n}({\alpha }_{h,n}^2+{\beta }_{h,n}^2)(z_h-z_{h-1})]\hfill \end{array}$$

  (13)
• The parameter ${\lambda }_{h,n}$ is characterized by the equation:

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

  (14)

  with the associated eigenvalues the eigenvalues ${\gamma }_n$, for $n\in {\mathbb{N}}^{\ast }$, associated with this problem are real and discrete. Additionally, the corresponding eigenfunctions $P_{h,n}$ are mutually orthogonal.

Lemma 1.

(Transformed${\mathcal{P}}_H$Problem Formulation)

The ${\mathcal{P}}_H$ problem, when the Fourier transform is applied to the lateral coordinate y, is reformulated in terms of the transformed variable ${\chi }_h^{\omega }$, leading to the subsequent system :

where ${\chi }_h^{\omega }(x,z)$ is defined as:

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

(15)

and ${\widehat{c}}_h^{\omega }(x,z)$ denotes the Fourier transformation of $c_h$ with respect to y, which is given by:

$${\widehat{c}}_h^{\omega }(x,z)={\int }_{-\infty }^{+\infty }{c}_h(x,y,z)e^{-2i\pi \omega y}dy,h\in \{1,\dots ,H\}$$

(16)

Proof of Lemma 1.

In proving the lemma, we adhere to the following steps:

• We begin by applying the Fourier transform with respect to the lateral coordinate y, as presented in Eqt 16, to the ${\mathcal{P}}_H$ problem. The resulting equation is:

  $$u_h\left[{\displaystyle \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}{\displaystyle \frac{{\partial }^2{\widehat{c}}_h^{\omega }}{\partial z^2}},z\in \left]z_{h-1},z_h\right[.$$

  (17)
• We calculate the derivative of ${\chi }_h^{\omega }(x,z)$ with respect to x, obtaining:

  $${\displaystyle \frac{\partial {\chi }_h^{\omega }}{\partial x}}=\left[{\displaystyle \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)ds\right).$$

  (18)
• We multiply both sides of Eqt 17 by the exponential term $exp\left({\left(2\pi \right)}^2{\omega }^2{\int }_0^x{\zeta }_y\left(s\right)ds\right)$. Given that $u_h$ and ${\phi }_{z_h}$ remain constants within each sub-layer, it simplifies to:

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

  (19)
• We apply the same method to the boundary conditions, expressing them in terms of ${\chi }_h^{\omega }$ to ensure consistency with the transformed system.

□

Having established the lemma, we now proceed to employ its results to demonstrate the theorem concerning the solution to the ${\mathcal{P}}_H$ problem.

Proof of Theorem 2.

The construction of the solution for the ${\mathcal{P}}_H$ problem unfolds through the following steps:

• Applying the lemma 1, we transform the ${\mathcal{P}}_H$ problem into a formulation that depends on ${\chi }_h^{\omega }$;
• Representing ${\chi }_h^{\omega }(x,z)$ as a series, we express it in a separated variables form:

  $${\chi }_h^{\omega }(x,z)=\sum _{n=0}^{\infty }{G}_{h,n}^{\omega }\left(x\right)P_{h,n}\left(z\right),h\in \{1,\cdots ,H\}.$$
• Verifying that the functions $G_{h,n}^{\omega }$ and $P_{h,n}$ satisfy their respective differential equations. We must therefore confirm that $G_{h,n}^{\omega }$ satisfies the equation

  $${\displaystyle \frac{d{G}_{h,n}^{\omega }}{dx}}+{\gamma }_n^2\xi \left(x\right)G_{h,n}^{\omega }=0,x\in ]0,L_x[$$

  (20)

  and $P_{h,n}\left(z\right)$ adheres to the system

  $$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle {\phi }_{z_h}\frac{d^2{P}_{h,n}}{d{z}^2}+{\gamma }_n^2{u}_h{P}_{h,n}=0;}z\in ]z_{h-1},z_h[,h\in \{1,\cdots ,H\}\hfill \\ \hfill \\ Undertheconditions:\left\{\begin{array}{c}{\displaystyle {\phi }_{z_1}\frac{d{P}_{1,n}\left(z_0\right)}{dz}=V_g{P}_{1,n}\left(z_0\right);}\hfill \\ \\ \left\{\begin{array}{cc}{\displaystyle P_{h-1,n}\left(z_{h-1}\right)=P_{h,n}\left(z_{h-1}\right);}\hfill & \\ \hfill & ,h\in \{2,\cdots ,H\}\hfill \\ {\displaystyle {\phi }_{z_{h-1}}\frac{d{P}_{h-1,n}\left(z_{h-1}\right)}{dz}={\phi }_{z_h}\frac{d{P}_{h,n}\left(z_{h-1}\right)}{dz};}\hfill \end{array}\right.\hfill \\ \\ {\phi }_{z_H}{\displaystyle \frac{d{P}_{H,n}\left(z_H\right)}{dz}}=0.\hfill \end{array}\right.\hfill \end{array}\right.\end{array}$$

  (21)
• Determining the specific forms of $G_{h,n}$ and $P_{h,n}$ ensures that they are respectively expressed as:

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

  (22)

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

  $$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)$$

  (23)

  where ${\lambda }_{h,n}$ is expressed in Eqt 14. This formulation must meet the system’s criteria, resulting in the determination of ${\alpha }_{h,n}$ and ${\beta }_{h,n}$, as specified in Eqts 11 and 12 respectively.
• The continuity condition between sub-layers necessitates a recursion of eigenvalues ${\lambda }_{h,n}$ and the eigenfunctions $P_{h,n}$, for $h\in \{1,\cdots ,H\}$. This process begins by determining the eigenvalues associated with the first layer, ${\lambda }_{1,n}$, using the conditions related to $P_{1,n}$. The resulting expressions are enumerated below:

  -

  The function $P_{1,n}\left(z\right)$ is expressed as:

  $$\begin{array}{c}\hfill P_{1,n}\left(z\right)=cos\left({\lambda }_{1,n}(z-z_0)\right)+v_{g}sin\left({\lambda }_{1,n}(z-z_0)\right)\end{array}$$

  (24)

  where ${\displaystyle v_g=\frac{V_g}{{\phi }_{z_1}{\lambda }_{1,n}}}$.

  -

  The norm of $P_{1,n}$ is given by:

  $$\parallel P_{1,n}{\parallel }^2=\left\{\begin{array}{cc}{u}_1(z_1-z_0);\hfill & \mathrm{when}{\gamma }_n=0\hfill \\ u_1\left(1+v_g^2\right){\displaystyle \frac{z_1-z_0}{2}};\hfill & \mathrm{when}{\gamma }_n\ne 0\hfill \end{array}\right.$$

  (25)

  The recursive generation of eigenvalues for the intermediate and final layers depends on the known values of the previous layers, in accordance with their respective conditions. This approach makes it possible to construct the associated eigenfunctions, succinctly maintaining continuity between the layers. For further details, please refer to Appendix Section A.1.
• The Sturm-Liouville theorem [1] allows for the expanded representation of ${\chi }_H^{\omega }$ as:

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

  where the eigenfunctions $P_{h,n}$ form a complete and mutually orthogonal set. Consequently, ${\chi }_H^{\omega }$ is detailed as above. The coefficients $a_n$ are determined leveraging the orthogonality of $P_{h,n}$, calculated as follows:

  $$\begin{array}{ccc}\hfill a_n& =\sum _{\ell =1}^H{\int }_{z_{\ell -1}}^{z_{\ell }}{u}_{\ell }{\chi }_{\ell }^{\omega }(0,z)P_{\ell ,n}\left(z\right)dz\hfill & \hfill =\left\{\begin{array}{c}{\displaystyle \frac{Q}{\sum _{h=1}^H{u}_h(z_h-z_{h-1})}},\mathrm{when}{\gamma }_n=0\hfill \\ Q\sum _{\ell =1}^H{P}_{\ell ,n}\left(H_s\right),\mathrm{when}{\gamma }_n\ne 0\hfill \end{array}\right.\end{array}$$

  (26)
• Applying the inverse Fourier transform to ${\chi }_H^{\omega }(x,z)$ from Equation (16), and considering that the inverse transform of

  $$exp\left({\left(2\pi \right)}^2{\omega }^2{\int }_0^x{\zeta }_y\left(s\right)ds\right)$$

  is given by

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

  leads to the final expression for the solution $c_H$ of the problem ${\mathcal{P}}_H$.

□

The Crosswind-Integrated Concentration (CIC) model is a tool for analyzing the dispersion of pollutants in the atmosphere, accounting for the horizontal wind component crucial for pollutant transport. It integrates the concentration across the lateral wind direction, y, yielding a two-dimensional concentration map along the wind direction (x) and crosswind direction (z). Mathematically, the CIC, $c_H^y(x,z)$, is defined by integrating the pollutant concentration $c_H(x,y,z)$ over y from $-\infty$ to $+\infty$. This integration reflects the dispersal of pollutants due to the horizontal wind and assesses their impact in both windward and lateral directions.

Corollary 1.

The Crosswind-Integrated Concentration (CIC), $c_H^y(x,z)$, is equal to ${\chi }_H^{\omega }(x,z)$ as follows:

$$\begin{array}{cc}\hfill c_H^y(x,z)=& Q·[\sum _{n=1}^{\infty }exp\left(-{\gamma }_n^2{\int }_0^x\xi \left(s\right)ds\right)P_{H,n}\left(z\right)·\{\sqrt{2}{\displaystyle \frac{cos\left({\lambda }_{1,n}(H_s-z_0)\right)+v_{g}sin\left({\lambda }_{1,n}(H_s-z_0)\right)}{\sqrt{u_1\left(1+v_g^2\right)(z_1-z_0)}}}+\hfill \\ \hfill & \sum _{h=2}^H{P}_{h,n}\left(H_s\right)\}+{\displaystyle \frac{1}{\sum _{h=1}^H{u}_h(z_h-z_{h-1})}}]\hfill \end{array}$$

(27)

Proof. 

To confirm the corollary, it suffices to verify the Gaussian integral identity:

$${\int }_{-\infty }^{+\infty }{\displaystyle \frac{1}{2\sqrt{\pi {\int }_0^x{\zeta }_y\left(s\right)ds}}}exp\left(-{\displaystyle \frac{y^2}{4{\int }_0^x{\zeta }_y\left(s\right)ds}}\right)dy=1.$$

□

Corollary 2.

The solution $c_H$ for the ${\mathcal{P}}_H$ problem encompasses and generalizes the previously established solutions under specific parameterizations:

• Single-Layer and Two-Dimensional Solution:

  (a)

  For a single layer without discretization, where ${\zeta }_y\equiv 1$ and $\xi \equiv 1$ , our method yields the solution documented by Kumar[15].

  (b)

  Transitioning to $v_g=0$, ${\zeta }_y\equiv 1$, $\xi \equiv 1$, and $z_0=0$ aligns with the solution documented by Seinfeld [25].

• Multi-Layer and Three-Dimensional Solution:

  (a)

  Setting $v_g=0$ retrieves the solution identified by Farhane et al [7].

  (b)

  Further simplification to $v_g=0$ and $\xi \equiv 1$ leads to the solution presented by Marie et al [19].

Proof. 

A detailed proof of this corollary, please refer to the AppendixB. □

Corollary 3.

In the scenario of a single-layer framework, where $z_0$ approaches 0 and $z_1$ extends to $H_{mix}$, the Crosswind-Integrated Concentration (CIC), $c_H^y(x,z)$, is defined as follows:

$$\begin{array}{cc}\hfill c_1^y(x,z)=& {\displaystyle \frac{Q}{\overline{u}{H}_{mix}}}[1+{\displaystyle \frac{2}{1+v_g^2}}\sum _{n=1}^{\infty }exp\left(-{\displaystyle \frac{{\left(n\pi \right)}^2\overline{{\phi }_z}}{\overline{u}{H}_{mix}^2}}x\right)\{\left(1+v_g^2\right)cos\left({\displaystyle \frac{n\pi }{H_{mix}}}(z-H_s)\right)+\hfill \\ \hfill & \left(1-v_g^2\right)cos\left({\displaystyle \frac{n\pi }{H_{mix}}}(z+H_s)\right)+2sin\left({\displaystyle \frac{n\pi }{H_{mix}}}(z+H_s)\right)\left\}\right]\hfill \end{array}$$

(28)

where:

• $\overline{u}$ represents the average wind speed over the height of the mixing layer, calculated as

  $$\overline{u}={\displaystyle \frac{1}{H_{mix}}}{\int }_0^{H_{mix}}u\left(s\right)ds,$$
• $\overline{{\phi }_z}$ denotes the average eddy diffusivity over the same height, expressed as

  $$\overline{{\phi }_z}={\displaystyle \frac{1}{H_{mix}}}{\int }_0^{H_{mix}}{\phi }_z\left(s\right)ds.$$

Proof. 

To validate the Crosswind-Integrated Concentration (CIC) for a scenario where $H=1$ (a single-layer framework), we follow these steps:

• Substitute Single-Layer Values into the CIC Equation: Begin by inserting the conditions of a single-layer setup into the original CIC formula from Eqt (27):

  $$\begin{array}{cc}\hfill c_1^y(x,z)=& Q\left[{\displaystyle \frac{1}{u_1(z_1-z_0)}}+2\sum _{n=1}^{\infty }exp\left(-{\gamma }_n^2{\int }_0^x\xi \left(s\right)ds\right)\left\{{\displaystyle \frac{cos\left({\lambda }_{1,n}(z-z_0)\right)+v_{g}sin\left({\lambda }_{1,n}(z-z_0)\right)}{u_1(1+v_g^2)(z_1-z_0)}}+\right.\right.\hfill \\ & \left.\left.{\displaystyle \frac{cos\left({\lambda }_{1,n}(H_s-z_0)\right)+v_{g}sin\left({\lambda }_{1,n}(H_s-z_0)\right)}{u_1(1+v_g^2)(z_1-z_0)}}\right\}\right]\hfill \end{array}$$

  where ${\displaystyle {\gamma }_n^2=\frac{{\left(n\pi \right)}^2\overline{{\phi }_z}}{\overline{u}{(z_1-z_0)}^2}}$.
• Apply Limits as $z_0$ Approaches 0 and $z_1$ Extends to $H_{mix}$: Evaluate the limits to reflect the shift to a single-layer framework that spans the entire mixing height $H_{mix}$:

  $$\begin{array}{cc}\hfill c_1^y(x,z)={\displaystyle \frac{Q}{\overline{u}{H}_{mix}}}& \left[1+{\displaystyle \frac{2}{1+v_g^2}}\sum _{n=1}^{\infty }exp\left(-{\displaystyle \frac{{\left(n\pi \right)}^2{\overline{\phi }}_z}{\overline{u}{H}_{mix}^2}}{\int }_0^x\xi \left(s\right)ds\right)\left\{\left(cos\left(n\pi {\displaystyle \frac{z}{H_{mix}}}\right)+v_{g}sin\left(n\pi {\displaystyle \frac{z}{H_{mix}}}\right)\right)\right.\right.\hfill \\ & \left.\left.·\left(cos\left(n\pi {\displaystyle \frac{z}{H_{mix}}}\right)+v_{g}sin\left(n\pi {\displaystyle \frac{z}{H_{mix}}}\right)\right)\right\}\right]\hfill \\ \hfill ={\displaystyle \frac{Q}{\overline{u}{H}_{mix}}}& [1+{\displaystyle \frac{2}{1+v_g^2}}\sum _{n=1}^{\infty }exp\left(-{\displaystyle \frac{{\left(n\pi \right)}^2\overline{{\phi }_z}}{\overline{u}{H}_{mix}^2}}{\int }_0^x\xi \left(s\right)ds\right)\{\left(1+v_g^2\right)cos\left({\displaystyle \frac{n\pi }{H_{mix}}}(z-H_s)\right)+\hfill \\ & \left(1-v_g^2\right)cos\left({\displaystyle \frac{n\pi }{H_{mix}}}(z+H_s)\right)+2sin\left({\displaystyle \frac{n\pi }{H_{mix}}}(z+H_s)\right)\left\}\right]\hfill \end{array}$$

□