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

We start this section by applying Fourier transform to Eq. (11). Let ^*c ω <sup>h</sup>* ð Þ *x*, *z* denote the Fourier transformation of *ch* with respect to *y*.

$$\hat{c}\_h^{\alpha}(\mathbf{x}, \mathbf{z}) = \int\_{-\infty}^{+\infty} c\_h(\mathbf{x}, \mathbf{y}, \mathbf{z}) e^{-2i\pi\alpha\mathbf{y}} d\mathbf{y}, \quad h \in \{1, \cdots, H\}^\*$$

which gives

$$
\mu\_h \left[ \frac{\partial \hat{c}\_h^{\mu}}{\partial \mathbf{x}} + (2\pi)^2 \alpha^2 \zeta\_\gamma(\mathbf{x}) \hat{c}\_h^{\mu} \right] = \xi(\mathbf{x}) \cdot \varrho\_{x\_h} \frac{\partial^2 \hat{c}\_h^{\mu}}{\partial \mathbf{z}^2}, \quad \mathbf{z} \in [z\_{h-1}, z\_h]. \tag{12}
$$

Let

$$\chi\_h^{\boldsymbol{\alpha}}(\boldsymbol{x},\boldsymbol{z}) = \hat{c}\_h^{\boldsymbol{\alpha}}(\boldsymbol{x},\boldsymbol{z}) \cdot \exp\left(\left(2\pi\right)^2 \left.\boldsymbol{\alpha}^2\right]\_0^{\boldsymbol{\alpha}} \zeta\_{\boldsymbol{\mathcal{V}}}(\boldsymbol{s}) d\boldsymbol{s}\right) \tag{13}$$

then

$$\frac{\partial \zeta\_h^{\alpha}}{\partial \mathbf{x}} = \left[ \frac{\partial \hat{c}\_h^{\alpha}}{\partial \mathbf{x}} + \left( 2\pi \right)^2 \alpha^2 \zeta\_\mathcal{\gamma}(\mathbf{x}) \hat{c}\_h^{\alpha} \right] \quad \exp \left( \left( 2\pi \right)^2 \alpha^2 \int\_0^\infty \zeta\_\mathcal{\gamma}(s) ds \right).$$

By multiplying both sides of Eq. (12) by exp 2ð Þ *<sup>π</sup>* <sup>2</sup> *<sup>ω</sup>*<sup>2</sup> Ð *x* <sup>0</sup> *<sup>ζ</sup>y*ð Þ*<sup>s</sup> ds* � � and since *uh* and *φzh* are constants for each interval, we show easily that, for all *h*∈ f g 1, ⋯, *H* :

$$u\_h \frac{\partial \chi\_h^o}{\partial \mathbf{x}} = \xi(\mathbf{x}) \ \begin{array}{c} \varrho\_{z\_h} \ \frac{\partial^2 \chi\_h^o}{\partial \mathbf{z}^2}, \ \quad \mathbf{z} \in [z\_{h-1}, z\_h] \end{array} \tag{14}$$

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

$$\begin{cases} \begin{aligned} &\varrho\_{x\_1}^{\alpha^{\rho}}(\mathbf{x},\mathbf{z}\_0) = \mathbf{0} \\ &\left\{\chi\_{h-1}^{\rho}(\mathbf{x},\mathbf{z}\_{h-1}) = \chi\_h^{\rho}(\mathbf{x},\mathbf{z}\_{h-1}) \right. \\ &\left.\varrho\_{x\_{h-1}}^{\rho}\frac{\partial \chi\_{h-1}^{\rho}(\mathbf{x},\mathbf{z}\_{h-1})}{\partial \mathbf{z}} = \varrho\_{x\_h}\frac{\partial \chi\_h^{\rho}(\mathbf{x},\mathbf{z}\_{h-1})}{\partial \mathbf{z}} \right. \\ &\left.\varrho\_{x\_{H}}^{\rho}\frac{\partial \chi\_H^{\rho}(\mathbf{x},\mathbf{z}\_{H})}{\partial \mathbf{z}} = \mathbf{0} \end{aligned} , h \in \{2,\cdots,H\} \end{cases}$$

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

$$\chi\_h^{\alpha}(\mathbf{x}, \mathbf{z}) = \sum\_{n=0}^{\infty} G\_{h,n}^{\alpha}(\mathbf{x}) \; P\_{h,n}(\mathbf{z}), \qquad h \in \{1, \cdots, H\}$$

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

$$\frac{d}{d\boldsymbol{\omega}}\frac{\mathbf{G}\_{h,n}^{a}}{d\boldsymbol{\omega}} + \boldsymbol{\gamma}\_{n}^{2} \; \; \xi(\boldsymbol{\omega}) \; \; \; \; \mathbf{G}\_{h,n}^{a} = \mathbf{0} \tag{15}$$

and

$$
\rho\_{x\_h} \frac{d^2}{dz^2} \frac{P\_{h,n}}{dz^2} + \gamma\_n^2 \ u\_h \ P\_{h,n} = 0,\tag{16}
$$

where *γ<sup>n</sup>* is a separation constant.

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

$$G\_{h,n}^{\
u}(\infty) = \mu\_n(\alpha) \, . \exp\left(-\gamma\_n^2 \int\_0^\kappa \xi(s)ds\right),$$

where *μ<sup>n</sup>* is an arbitrary function depending on *ω*.

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

$$P\_{h,n}(z) = a\_{h,n} \cos(\lambda\_{h,n} z) + \beta\_{h,n} \sin(\lambda\_{h,n} z),\tag{17}$$

where, *λ<sup>h</sup>*,*<sup>n</sup>* ¼ *γ<sup>n</sup>* ffiffiffiffiffi *uh φzh* q .

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

$$\begin{cases} \begin{aligned} \,\rho\_{x\_{1}}\frac{dP\_{1,\mathfrak{n}}(x\_{0})}{dx} = 0 \quad &(\mathbf{a})\\ \begin{cases} \,P\_{h-1,\mathfrak{n}}(x\_{h-1}) = P\_{h,\mathfrak{n}}(x\_{h-1})\\ \,\rho\_{x\_{h-1}}\frac{dP\_{h-1,\mathfrak{n}}(x\_{h-1})}{dx} = \,\rho\_{x\_{h}}\frac{dP\_{h,\mathfrak{n}}(x\_{h-1})}{dx}, h \in \{2,\cdots,H\}(\mathbf{b}) \end{cases} \end{cases} \end{cases} \tag{18}$$

To calculate the expression of *Ph*,*<sup>n</sup>*, it comes down to calculate the values of *α<sup>h</sup>*,*<sup>n</sup>* and *β<sup>h</sup>*,*<sup>n</sup>*, on each of the sub-layer *zh*�<sup>1</sup> ½ � , *zh* , *h*∈f g 1, ⋯, *H* .

By solving the recursive system resulting from substitution of Eq. (17) in Eq. (18), we obtain respectively the formulations of *α<sup>h</sup>*,*<sup>n</sup>* and *β<sup>h</sup>*,*<sup>n</sup>*. More specifically,

For the first sub-layer, *α*1,*<sup>n</sup>* and *β*1,*<sup>n</sup>* satisfy the equation:

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

from which, we can take *β*1,*<sup>n</sup>* ¼ sinð Þ *λ*1,*nz*<sup>0</sup> , so that *α*1,*<sup>n</sup>* ¼ cosð Þ *λ*1,*nz*<sup>0</sup> , which means:

$$P\_{1,n}(z) = \cos(\lambda\_{1,n}(z - z\_0))\,.$$

*A Novel Development in Three-Dimensional Analytical Solutions for Air Pollution… DOI: http://dx.doi.org/10.5772/intechopen.112225*

For the last sub-layer (*Hth* sub-layer):

$$a\_{H, \mathfrak{n}} = \cot(\lambda\_{H, \mathfrak{n}} \mathfrak{z}\_H) \; \beta\_{H, \mathfrak{n}}.\tag{20}$$

And for the intermediate sub-layers:

$$\begin{split} a\_{h,n} &= \frac{\varrho\_{x\_{k}}\lambda\_{h,n} - \varrho\_{x\_{k-1}}\lambda\_{h-1,n}}{2\varrho\_{x\_{k}}\lambda\_{h,n}} \Bigg[ \Bigl( \cos((\lambda\_{h,n} + \lambda\_{h-1,n})z\_{h-1}) + \frac{\varrho\_{x\_{k}}\lambda\_{h,n} + \varrho\_{x\_{k-1}}\lambda\_{h-1,n}}{\varrho\_{x\_{k}}\lambda\_{h,n} - \varrho\_{x\_{k-1}}\lambda\_{h-1,n}} \\ &\quad \cos((\lambda\_{h,n} - \lambda\_{h-1,n})z\_{h-1}))a\_{h-1,n} + (\sin((\lambda\_{h,n} + \lambda\_{h-1,n})z\_{h-1}) - \\ &\quad \frac{\varrho\_{x\_{k}}\lambda\_{h,n} + \varrho\_{x\_{k-1}}\lambda\_{h-1,n}}{\varrho\_{x\_{k}}\lambda\_{h,n} - \varrho\_{x\_{k-1}}\lambda\_{h-1,n}} \sin((\lambda\_{h,n} - \lambda\_{h-1,n})z\_{h-1}) \Bigg) \beta\_{h-1,n} \Bigg], \end{split} \tag{21}$$

and

$$\begin{split} \beta\_{h,n} &= \frac{\varrho\_{x\_{h}} \lambda\_{h,n} - \varrho\_{x\_{h-1}} \lambda\_{h-1,n}}{2\varrho\_{x\_{h}} \lambda\_{h,n}} \Bigg[ \Bigl( \sin((\lambda\_{h,n} + \lambda\_{h-1,n})x\_{h-1}) + \frac{\varrho\_{x\_{h}} \lambda\_{h,n} + \varrho\_{x\_{h-1}} \lambda\_{h-1,n}}{\varrho\_{x\_{h}} \lambda\_{h,n} - \varrho\_{x\_{h-1}} \lambda\_{h-1,n}} \\ & \quad \sin((\lambda\_{h,n} - \lambda\_{h-1,n})x\_{h-1})) a\_{h-1,n} - (\cos((\lambda\_{h,n} + \lambda\_{h-1,n})x\_{h-1}) - \\ & \quad \frac{\varrho\_{x\_{h}} \lambda\_{h,n} + \varrho\_{x\_{h-1}} \lambda\_{h-1,n}}{\varrho\_{x\_{h}} \lambda\_{h,n} - \varrho\_{x\_{h-1}} \lambda\_{h-1,n}} \cos((\lambda\_{h,n} - \lambda\_{h-1,n})x\_{h-1}) \Big) \boldsymbol{\beta}\_{h-1,n} \Bigg]. \end{split} \tag{22}$$

The eigenvalues *γn*, *n* ∈ℕ<sup>∗</sup> 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 *uh* on each interval *zh*�<sup>1</sup> ½ � , *zh* , *h*∈ f g 1, ⋯, *H* leads to

$$\sum\_{h=1}^{H} \int\_{x\_{h-1}}^{x\_h} u\_h \, P\_{h,m}(s) \, P\_{h,n}(s) \, \, ds = \|P\_{h,m}\|. \quad \|P\_{h,n}\|. \quad \delta\_{m,n},$$

where *δ<sup>m</sup>*,*<sup>n</sup>* is the Kronecker symbol. Then, ∥*Ph*,*<sup>n</sup>*∥<sup>2</sup> is written as follows.

$$\begin{aligned} \left\| \left| P\_{H,n} \right| \right\|^2 &= \sum\_{h=1}^{H} \int\_{z\_{h-1}}^{z\_h} u\_h \left( P\_{h,n} (s) \right)^2 ds \\ &= \begin{cases} \sum\_{h=1}^{H} u\_h (z\_h - z\_{h-1}), & \text{when } \ \gamma\_n = 0 \\\\ \sum\_{h=1}^{H} \frac{u\_h}{2} \left[ \sin \left( \lambda\_{h,n} (z\_h - z\_{h-1}) \right) ( (a\_{h,n}^2 - \beta\_{h,n}^2) \cos \left( \lambda\_{h,n} (z\_h + z\_{h-1}) \right) + \\\\ 2a\_{h,n} \beta\_{h,n} \sin \left( \lambda\_{h,n} (z\_h + z\_{h-1}) \right) \right) + \lambda\_{h,n} (a\_{h,n}^2 + \beta\_{h,n}^2) (z\_h - z\_{h-1}) \right], & \text{when } \ \gamma\_n \neq 0. \end{cases} \end{aligned}$$

The eigenvalues of each sub-layer can be obtained by integrating Eq. (17) on each of the intervals *zh*�<sup>1</sup> ½ � , *zh* , *h*∈f g 1, ⋯, *H* taking into account the boundary conditions Eqs. (19)–(21). Whereas the eigenfunctions *PH*,*<sup>n</sup>* form a complete set, so *χ<sup>ω</sup> <sup>H</sup>* can be developed as:

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

Then, the coefficients *an* are given by:

$$\begin{split} \boldsymbol{a}\_{n} &= \frac{1}{\|\boldsymbol{P}\_{h,n}\|^{2}} \sum\_{\ell=1}^{H} \int\_{\boldsymbol{z}\_{\ell-1}}^{\boldsymbol{z}\_{\ell}} \boldsymbol{u}\_{\ell} \boldsymbol{\chi}\_{\ell}^{\boldsymbol{\alpha}}(\mathbf{0}, \boldsymbol{z}) \boldsymbol{P}\_{\ell,n}(\boldsymbol{z}) d\boldsymbol{z} \\ &= \begin{cases} \frac{\boldsymbol{Q}}{\sum\_{h=1}^{H} \boldsymbol{u}\_{h} (\boldsymbol{z}\_{h} - \boldsymbol{z}\_{h-1})}, & \text{when } \boldsymbol{\chi}\_{\boldsymbol{n}} = \mathbf{0} \\\ \frac{\boldsymbol{Q}}{\|\boldsymbol{P}\_{h,n}\|^{2}} \sum\_{\ell'=1}^{H} \boldsymbol{P}\_{\ell,n}(\boldsymbol{H}\_{\ell}), & \text{when } \boldsymbol{\chi}\_{\boldsymbol{n}} \neq \mathbf{0} \end{cases} \end{split}$$

The inverse Fourier transform of exp 2ð Þ *<sup>π</sup>* <sup>2</sup> *<sup>ω</sup>*<sup>2</sup> Ð *x* <sup>0</sup> *<sup>ζ</sup>y*ð Þ*<sup>s</sup> ds* � � is:

$$\frac{1}{2\sqrt{\pi \int\_0^\infty \zeta\_\mathcal{\nu}(s)ds}} \exp\left(-\frac{\mathcal{y}^2}{4\int\_0^\kappa \zeta\_\mathcal{\nu}(s)ds}\right).$$

Consequently,

$$\begin{split} c\_{H}(\mathbf{x},\boldsymbol{y},\boldsymbol{z}) &= \frac{\boldsymbol{Q}}{2\sqrt{\pi \int\_{0}^{\boldsymbol{x}} \zeta\_{\boldsymbol{y}}(\boldsymbol{s})ds}} \exp\left(-\frac{\boldsymbol{y}^{2}}{4\int\_{0}^{\boldsymbol{x}} \zeta\_{\boldsymbol{y}}(\boldsymbol{s})ds}\right). \\ \sum\_{n=0}^{\infty} \exp\left(-\boldsymbol{\gamma}\_{n}^{2} \int\_{0}^{\boldsymbol{x}} \xi(\boldsymbol{s})ds\right) \frac{P\_{H,n}(\boldsymbol{z})}{||P\_{H,n}||^{2}} \left(\sum\_{h=1}^{H} P\_{h,n}(H\_{\boldsymbol{s}})\right). \end{split} \tag{23}$$

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 y <sup>H</sup>*ð Þ *x*, *z* is obtained by integrating *cH*ð Þ *x*, *y*, *z* in Eq. (20) to *y* from �∞ to þ∞:

$$c\_H^\mathcal{v}(\mathbf{x}, \mathbf{z}) = \int\_{-\infty}^{+\infty} c\_H(\mathbf{x}, \mathbf{y}, \mathbf{z}) d\mathbf{y}$$

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

$$\mathcal{L}\_H^r(\mathbf{x}, \mathbf{z}) = \mathbf{Q} \sum\_{n=0}^{\infty} \exp\left(-\gamma\_n^2 \int\_0^\mathbf{x} \xi(s) ds\right) \frac{P\_{H,n}(\mathbf{z})}{\left\| P\_{H,n} \right\|^2} \sum\_{h=1}^H P\_{h,n}(H\_s). \tag{24}$$

### **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

*A Novel Development in Three-Dimensional Analytical Solutions for Air Pollution… DOI: http://dx.doi.org/10.5772/intechopen.112225*

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(x) = \frac{u\_\*}{k} \left[ \ln \left( \frac{z + x\_0}{z\_0} \right) + 5.75 \left( \frac{z}{H\_{\text{mix}}} \right) - 1.88 \left( \frac{z}{H\_{\text{mix}}} \right)^2 - 1.33 \left( \frac{z}{H\_{\text{mix}}} \right)^3 + 0.25 \left( \frac{z}{H\_{\text{mix}}} \right)^4 \right],\tag{25}$$

where *u*<sup>∗</sup> is the friction velocity, *k* is the von Karman constant, *z*<sup>0</sup> is the roughness length, and *Hmix* 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 *Kz*, given in Eq.(8), is adopted from [11], where the functional *φz*ð Þ*z* is expressed as:

$$\rho\_x(\mathbf{z}) = 0.22 H\_{\text{mix}} w\_\* \left( \frac{\mathbf{z}}{H\_{\text{mix}}} \left( \mathbf{1} - \frac{\mathbf{z}}{H\_{\text{mix}}} \right) \right)^{1/3} \left( \mathbf{1} - \exp \left( -\frac{4\mathbf{z}}{H\_{\text{mix}}} \right) - 0.0003 \exp \left( \frac{8z}{H\_{\text{mix}}} \right) \right), \tag{26}$$

where *w*<sup>∗</sup> is the convective velocity. The integrable correction dimensionless function in Eq. (8) is defined in terms of the along-wind length scale *L*<sup>1</sup> as follows [10]:

$$\xi(\varkappa) = 1 - \exp\left(-\frac{\varkappa}{L\_1}\right). \tag{27}$$

The length *L*<sup>1</sup> is given in terms of *u*, *φ<sup>z</sup>* and *σ<sup>w</sup>* as [10]:

$$L\_1 = \frac{1}{\sigma\_w^2(H\_s)} \,\mu(H\_s) \,\rho\_x(H\_s),$$

where *σ<sup>w</sup>* is the vertical turbulent intensity. Note that, there exist many expressions of *σw*, we adopt here the expression given by [17]:

$$
\sigma\_w = 0.96 w\_\* \left( \frac{3z}{H\_{\rm mix}} + \frac{|L|}{H\_{\rm 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., ∣*L*∣ ≫ 1) indicates a neutral atmosphere, where turbulence is neither generated nor suppressed due to buoyancy effects [18].

It should be noted that when *<sup>L</sup>*<sup>1</sup> ! 0, the term exp � *<sup>x</sup> L*1 � � becomes negligible in (27), and *Kz* will dependent only on *z* (i.e., *Kz* � *φ<sup>z</sup>* in Eq. (8)).

The lateral eddy diffusivity is given by [12]:

$$K\_{\mathcal{Y}}(\mathbf{x}, \mathbf{z}) = \frac{1}{2} \ \!\!\begin{array}{c} \text{ $\mathbf{z}$ } \ \!\begin{array}{c} \text{ $\mathbf{z}$ } \ \text{ $\sigma\_{\mathcal{Y}}^2(\mathbf{x})$ } \\ \text{ $\mathbf{x}$ } \end{array} , \tag{28}$$

where *σ<sup>y</sup>* is the standard deviation in the crosswind direction (depends only on *x*). By identifying Eqs. (7) and (28), we obtain

$$
\zeta\_\circ(\infty) = \frac{1}{2} \left. \frac{d\,\sigma\_\circ^2(\infty)}{d\infty} \right.
$$
