2 General problem formulation

The turbulent dispersion of pollutants in the PBL is governed by the advection– diffusion equation

(2.1) $$ \begin{align} \dfrac{\partial C}{\partial t}+\nabla\mathbf{.} (\mathbf{U_{w}}\; C)= \nabla\mathbf{ .}( D \;\nabla C)+ S, \end{align} $$

where $\displaystyle {\mathbf {U_{w}} =(U,V,W)^{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; $\displaystyle {D}$ is the molecular diffusion coefficient; $\displaystyle {S}$ is the source term; and $\displaystyle {\nabla }$ is the gradient operator.

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

$$ \begin{align*} \mathbf{U_{w}}=\overline{\mathbf{U_{w}} }+\mathbf{U_{w}}^{^{\prime}}\quad \mbox{with}\quad \overline{\mathbf{U_{w}} }=(u,v,w)^{T}\; \quad\mbox{and}\;\quad \mathbf{ U_{w}}^{^{\prime}}=(u^{\prime},v^{\prime},w^{\prime})^{T}. \end{align*} $$

The application of the Reynolds averaging rules to the vertical mass flow, $\mathbf { U_{w}}\; C$ , leads to [Reference De Visscher6]

(2.2) $$ \begin{align} \overline{\mathbf{ U_{w}}\; C}=\overline{\mathbf{U_{w}} }\; c + \overline{\mathbf{ U_{w}}^{^{\prime}}\;c^{\prime}}\quad \mbox{with}\quad \overline{\mathbf{U_{w}} {}^{^{\prime}}\;c^{\prime}}=(\overline{u^{\prime}c^{\prime}},\overline{v^{\prime}c^{\prime}},\overline{w^{\prime}c^{\prime}})^{T}. \end{align} $$

It turns out that turbulent diffusion can be described with Fick’s laws of diffusion as follows [Reference Pasquill and Smith22]:

$$ \begin{align*}\nonumber \overline{u^{\prime}c^{\prime}}=- K_{x}\,\dfrac{\partial c}{\partial x},\quad\nonumber \overline{v^{\prime}c^{\prime}}=- K_{y}\;\dfrac{\partial c}{\partial y},\quad\nonumber \overline{w^{\prime}c^{\prime}}=- K_{z}\;\dfrac{\partial c}{\partial z},\nonumber \end{align*} $$

where $K_{x}$ , $K_{y}$ and $K_{z}$ are the eddy diffusivity components along the x-, y- and $ z$ -directions, respectively.

Note 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 $\displaystyle {\nabla .(D \;\nabla C)}$ is then to be replaced by an eddy or turbulent diffusivity. The source term could be eliminated from equation (2.1), and should be added to the boundary conditions as a delta function. At the point $(0, 0, H_{s})$ , there is a source rejecting the pollutant with a continuous flow Q,

$$ \begin{align*} u \; c(0,y,z)=Q\delta(y)\delta(z-H_{s}), \end{align*} $$

where $H_{s}$ is the source height. By application of the Reynolds averaging and the divergence operator to equation (2.2), equation (2.1) may be written as

(2.3) $$ \begin{align} \dfrac{\partial c}{\partial t}=- \dfrac{\partial}{\partial x} \overline{u^{\prime}c^{\prime}}- \dfrac{\partial}{\partial y} \overline{v^{\prime}c^{\prime}}- \dfrac{\partial}{\partial z} \overline{w^{\prime}c^{\prime}}-\overline{U_{w}}.\nabla c. \end{align} $$

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

1. (a) the steady state condition (that is, ${\partial c}/{\partial t}=0$ );

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

3. (c) the turbulent diffusion in the direction of the mean wind is neglected compared to the advection transport mechanism, that is,

   $$ \begin{align*}u \dfrac{\partial c}{\partial x}>> \dfrac{\partial }{\partial x} \bigg(K_{x}\dfrac{\partial c}{\partial x}\bigg).\end{align*} $$

These assumptions lead to the steady-state advection–diffusion equation defined as $0<x< L_{x} $ , $ -L_{y}<y< L_{y}$ and $0<z< H_{\text {mix}}$ ,

$$ \begin{align*} u(z) \dfrac{\partial c}{\partial x}= \dfrac{\partial }{\partial y} \bigg(K_{y}\dfrac{\partial c}{\partial y}\bigg)+\dfrac{\partial }{\partial z} \bigg(K_{z}\dfrac{\partial c}{\partial z}\bigg), \end{align*} $$

which is subject to the boundary conditions:

(2.4) $$ \begin{align} \begin{array}{l} \displaystyle{\lim_{L_{y}\longrightarrow \infty}}K_{y}(x,z)\;\dfrac{\partial c}{\partial y} = 0,\quad K_{z}(x,z)\;\dfrac{\partial c}{\partial z} =0\; \quad\mbox{at } z=z_{0}, \quad \mbox{and}\\[10pt] K_{z}(x,z)\;\dfrac{\partial c}{\partial z} =0 \quad \mbox{at } z= H_{\text{mix}}, \end{array} \end{align} $$

where $z_{0}$ is the surface roughness length and $H_{\text {mix}}$ is the PBL height.

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

(2.5) $$ \begin{align} K_{y}(x,z)&= \zeta_{y}(x)\; u(z), \end{align} $$

(2.6) $$ \begin{align} K_{z}(x,z)&= \xi(x)\;\varphi_{z}(z). \end{align} $$

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, \ldots , H$ ,

(2.7) $$ \begin{align} u_{h} &=\frac{1}{z_{h}-z_{h-1}} \int_{z_{h-1}}^{z_{h}} u(s)\, ds, \end{align} $$

(2.8) $$ \begin{align} \varphi_{z_{h}} &= \frac{1}{z_{h}-z_{h-1}} \int_{z_{h-1}}^{z_{h}} \varphi_{z}(s)\, ds. \end{align} $$

Using the formulations of $K_{y}$ and $K_{z}$ in equations (2.5) and (2.6), equation (2.3) can be written as

(2.9) $$ \begin{align} u_{h}\;\dfrac{\partial c_{h}}{\partial x} = \zeta_{y}(x)\; u_{h}\;\dfrac{\partial^{2} c_{h}}{\partial y^{2}}+\xi(x)\;\varphi_{z_{h}}\;\dfrac{\partial^{2} c_{h}}{\partial z^{2}} \end{align} $$

with $u_{h}$ and $\varphi _{z_{h}}$ (given by equations (2.7) and (2.8)) as constants. Equation (2.9) is subject to the first boundary conditions of equation (2.4) on the one hand, and on the other hand, the continuity of both the concentration and the flux at the interface level is applied. For $h\in \{2,\ldots ,H\},$

$$ \begin{align*} \begin{cases} \varphi_{z_{1}}\;\dfrac{\partial c_{1}(x,y,z_{0})}{\partial z}= 0,\quad \varphi_{z_{H}}\;\dfrac{\partial c_{H}(x,y,z_{H})}{\partial z}= 0, \\[6pt] c_{h-1}(x,y,z_{h-1})=\; c_{h}(x,y,z_{h-1}),\\[9pt] \varphi_{z_{h-1}}\;\dfrac{\partial c_{h-1}(x,y,z_{h-1})}{\partial z} = \varphi_{z_{h}}\;\dfrac{\partial c_{h}(x,y,z_{h-1})}{\partial z}. \\ \end{cases} \end{align*} $$

3 Analytical solution

We start this section by applying Fourier transform to equation (2.9). Let $\hat {c}_{h}^{\omega }(x,z)$ denote the Fourier transformation of $c_{h}$ with respect to y. Then,

$$ \begin{align*} \hat{c}_{h}^{\omega}(x,z)=\int_{-\infty}^{+\infty} c_{h}(x,y,z)\: e^{-2 i \pi \omega y} \,dy,\quad h\in\{1,\ldots,H\}, \end{align*} $$

which gives

$$ \begin{align*} u_{h} \bigg[ \frac{\partial \hat{c}_{h}^{\omega}}{\partial x}+(2 \pi)^{2} \omega^{2}\zeta_{y}(x) \hat{c}_{h}^{\omega}\bigg]=\xi(x)\;\varphi_{z_{h}}\;\dfrac{\partial^{2} \hat{c}_{h}^{\omega}}{\partial z^{2}},\quad z\in [z_{h-1},z_{h}]. \end{align*} $$

Let

(3.1) $$ \begin{align} \chi_{h}^{\omega}(x,z)=\hat{c}_{h}^{\omega}(x,z)\exp\!\bigg((2 \pi )^{2}\; \omega^{2} \int_{0}^{x} \zeta_{y}(s)\, ds\bigg), \end{align} $$

then

$$ \begin{align*} \frac{\partial \chi_{h}^{\omega}}{\partial x}=\bigg[ \frac{\partial \hat{c}_{h}^{\omega}}{\partial x}+(2 \pi )^{2}\omega^{2}\zeta_{y}(x) \hat{c}_{h}^{\omega}\bigg]\;\exp\!\bigg((2 \pi )^{2}\; \omega^{2} \int_{0}^{x} \zeta_{y}(s)\, ds\bigg). \end{align*} $$

By multiplying both sides of equation (3.1) by

$$ \begin{align*}\displaystyle \exp\!\bigg((2 \pi )^{2}\; \omega^{2} \int_{0}^{x} \zeta_{y}(s)\, ds\bigg)\end{align*} $$

and, since $u_{h}$ and $\varphi _{z_{h}}$ are constants for each interval, we show easily that for all $h\in \{1,\ldots ,H\}$ ,

(3.2) $$ \begin{align} u_{h} \frac{\partial \chi_{h}^{\omega}}{\partial x}=\xi(x)\;\varphi_{z_{h}}\;\dfrac{\partial^{2} \chi_{h}^{\omega}}{\partial z^{2}},\quad z\in [z_{h-1},z_{h}]. \end{align} $$

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

The solution of equation (3.2) is assumed to be in the form

$$ \begin{align*} \chi_{h}^{\omega}(x, z)=\sum_{n=0}^{\infty} G_{h,n}^{\omega}(x)\; P_{h,n}(z), \quad h\in\{1,\ldots,H\}. \end{align*} $$

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

(3.3) $$ \begin{align} \frac{d\; G_{h,n}^{\omega}}{d x} +\gamma_{n}^{2} \; \xi(x) \; G_{h,n}^{\omega}=0, \end{align} $$

and

(3.4) $$ \begin{align} \varphi_{z_{h}}\;\frac{d^{2}\; P_{h,n}}{d z^{2}} +\gamma_{n}^{2} \; u_{h} \; P_{h,n}=0, \end{align} $$

where $\displaystyle \gamma _{n}$ is a separation constant.

The first-order ordinary differential equation (3.3) has the solution

$$ \begin{align*} G_{h,n}^{\omega}(x)=\mu_{n}(\omega)\exp\!\bigg(-\gamma_{n}^{2} \int_{0}^{x} \xi(s)\, ds\bigg), \end{align*} $$

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

Equation (3.4) represents a Sturm–Liouville problem. Solutions of such problem form an eigenfunction basis of the form

(3.5) $$ \begin{align} P_{h,n}(z)=\alpha_{h,n}\cos(\lambda_{h,n} z)+\beta_{h,n}\sin(\lambda_{h,n} z), \end{align} $$

where $\displaystyle {\lambda _{h,n}=\gamma _{n}\sqrt {(u_{h}/\varphi _{z_{h}})}}$ .

Equation (3.5) satisfies the following boundary conditions:

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

By solving the recursive system resulting from substitution of equation (3.5) in equations (3.6)–(3.8), we obtain respectively the formulations of $\displaystyle {\alpha _{h,n}}$ and $\displaystyle {\beta _{h,n}}$ . More specifically, we have the following.

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

$$ \begin{align*} \alpha_{1,n}\sin(\lambda_{1,n} z_{0})-\beta_{1,n} \cos(\lambda_{1,n} z_{0})=0, \end{align*} $$

from which we can take $\displaystyle {\beta _{1,n} =\sin (\lambda _{1,n} z_{0})}$ , so that $\displaystyle {\alpha _{1,n} =\cos (\lambda _{1,n} z_{0})}$ , which means

$$ \begin{align*} P_{1,n}(z)=\cos(\lambda_{1,n}(z-z_{0})). \end{align*} $$

For the last sub-layer ( $H^{\text {th}}$ sub-layer),

$$ \begin{align*} \alpha_{H,n}=\cot(\lambda_{H,n} z_{H})\;\beta_{H,n}. \end{align*} $$

Additionally for the intermediate sub-layers,

$$ \begin{align*} \alpha_{h,n}&=\dfrac{\varphi_{z_{h}} \lambda_{h,n}-\varphi_{z_{h-1}} \lambda_{h-1,n} }{2 \varphi_{z_{h}} \lambda_{h,n}} \bigg[\bigg\{\cos((\lambda_{h,n}+\lambda_{h-1,n})z_{h-1}) +\dfrac{\varphi_{z_{h}} \lambda_{h,n}+\varphi_{z_{h-1}} \lambda_{h-1,n} }{\varphi_{z_{h}} \lambda_{h,n}-\varphi_{z_{h-1}} \lambda_{h-1,n} } \nonumber\\[3pt] &\quad \times \cos((\lambda_{h,n}-\lambda_{h-1,n})z_{h-1}) \bigg\}\alpha_{h-1,n} +\bigg\{\sin((\lambda_{h,n}+\lambda_{h-1,n})z_{h-1}) \nonumber\\[3pt] &\quad -\dfrac{\varphi_{z_{h}} \lambda_{h,n}+\varphi_{z_{h-1}} \lambda_{h-1,n} }{\varphi_{z_{h}} \lambda_{h,n}-\varphi_{z_{h-1}} \lambda_{h-1,n} } \sin((\lambda_{h,n}-\lambda_{h-1,n})z_{h-1}) \bigg\}\beta_{h-1,n} \bigg], \end{align*} $$

and

$$ \begin{align*} \beta_{h,n}&=\dfrac{\varphi_{z_{h}} \lambda_{h,n}-\varphi_{z_{h-1}} \lambda_{h-1,n} }{2 \varphi_{z_{h}} \lambda_{h,n}}\bigg[\bigg\{\sin((\lambda_{h,n}+\lambda_{h-1,n})z_{h-1})+\dfrac{\varphi_{z_{h}} \lambda_{h,n}+\varphi_{z_{h-1}} \lambda_{h-1,n} }{\varphi_{z_{h}} \lambda_{h,n}-\varphi_{z_{h-1}} \lambda_{h-1,n} } \nonumber\\[3pt] &\quad \times \sin((\lambda_{h,n}-\lambda_{h-1,n})z_{h-1}) \bigg\}\alpha_{h-1,n}-\bigg\{\cos((\lambda_{h,n}+\lambda_{h-1,n})z_{h-1}) \nonumber\\[3pt] &\quad -\dfrac{\varphi_{z_{h}} \lambda_{h,n}+\varphi_{z_{h-1}} \lambda_{h-1,n} }{\varphi_{z_{h}} \lambda_{h,n}-\varphi_{z_{h-1}} \lambda_{h-1,n} }\cos((\lambda_{h,n}-\lambda_{h-1,n})z_{h-1}) \bigg\}\beta_{h-1,n} \bigg]. \end{align*} $$

The eigenvalues $\displaystyle \gamma _{n},\; n\in \mathbb {N}^{*}$ of this problem are real and discrete, and the eigenfunctions are mutually orthogonal. The orthogonality relation developed by Mikhailov and Ozisik [Reference Mikhailov and Ozisik17] for this class of (self-adjoint) problems with respect to the density $u_{h}$ on each interval $[z_{h-1},z_{h}] , h\in \{1,\ldots , H\}$ leads to

$$ \begin{align*} \sum_{h=1}^{H}\: \int_{z_{h-1}}^{z_{h}} u_{h}\: P_{h,m}(s)\: P_{h,n}(s)\, ds= \|P_{h,m}\| \,.\, \|P_{h,n}\| \,.\, \delta_{m,n}, \end{align*} $$

where $\displaystyle {\delta _{m,n}}$ is the Kronecker symbol. Then, we can write

$$ \begin{align*} \|P_{H,n}\|^{2}&=\displaystyle \sum_{h=1}^{H}\:\int_{z_{h-1}}^{z_{h}} u_{h}\: ( P_{h,n}(s))^{2}\, ds\\[3pt] &= \begin{cases} \displaystyle \sum\limits_{h=1}^{H}\;u_{h}(z_{h}-z_{h-1}),\quad \mbox{when } \gamma_{n}=0 , \\[12pt] \displaystyle \sum\limits_{h=1}^{H}\;\dfrac{u_{h}}{2\; \lambda_{h,n}} [\sin(\lambda_{h,n}(z_{h}-z_{h-1}))\{(\alpha_{h,n}^{2}-\beta_{h,n}^{2})\cos(\lambda_{h,n}(z_{h}+z_{h-1}))\\[12pt] +2\alpha_{h,n} \,\beta_{h,n}\sin(\lambda_{h,n}(z_{h}\!+\!z_{h-1}))\}\!+\!\lambda_{h,n} (\alpha_{h,n}^{2}\!+\!\beta_{h,n}^{2})(z_{h}-z_{h-1}) ], \quad \mbox{when } \gamma_{n}\neq 0. \end{cases} \end{align*} $$

The eigenvalues of each sub-layer can be obtained by integrating equation (3.4) on each of the intervals $[z_{h-1},z_{h}], \; h\in \{1,\ldots , H\}$ , taking into account the boundary conditions equations (3.6)–(3.8). However, the eigenfunctions $\displaystyle {P_{H,n}}$ form a complete set, so $\displaystyle {\chi _{H}^{\omega }}$ can be developed as

$$ \begin{align*} \chi_{H}^{\omega}(x, z)=\sum_{n=0}^{\infty} a_{n}\exp\!\bigg(-\gamma_{n}^{2} \int_{0}^{x} \xi(s)\, ds\bigg)\; P_{H,n}(z). \end{align*} $$

Then, the coefficients $a_{n}$ are given by

$$ \begin{align*} a_{n}&=\dfrac{1}{\|P_{h,n}\|^{2}} \sum_{\ell=1}^{H} \int_{z_{\ell-1}}^{z_{\ell}} u_{\ell}\chi_{\ell}^{\omega}(0, z)\:P_{\ell,n}(z)\,dz\nonumber\\[3pt] &= \begin{cases} \dfrac{Q}{{\sum\nolimits_{h=1}^{H}\;u_{h}(z_{h}-z_{h-1})}}, \quad \mbox{when }\gamma_{n}=0, \\[12pt] \dfrac{Q}{{ \|P_{h,n}\|^{2}}}\: \sum\limits_{\ell=1}^{H} {P_{\ell,n}(H_{s})}, \quad \mbox{when }\gamma_{n}\neq0. \end{cases} \end{align*} $$

The inverse Fourier transform of

$$ \begin{align*}\displaystyle{\exp\!\bigg((2 \pi )^{2}\; \omega^{2} \int_{0}^{x} \zeta_{y}(s)\, ds\bigg)}\end{align*} $$

is given by

$$ \begin{align*} \frac{1}{2 \sqrt{\pi\:\int_{0}^{x}\zeta_{y}(s)\, ds} } \exp\!\bigg(-\frac{y^{2}}{4\:\int_{0}^{{x}}\zeta_{y}(s)\, ds } \bigg). \end{align*} $$

Consequently,

(3.9) $$ \begin{align}\nonumber c_{H}(x,y,z)&=\frac{Q}{2 \sqrt{\pi\:\int_{0}^{x}\zeta_{y}(s)\, ds } } \exp\!\bigg(-\frac{y^{2}}{4\:\int_{0}^{{x}}\zeta_{y}(s)\, ds }\bigg) \\[3pt] &\quad \times \sum_{n=0}^{\infty} \exp\!\bigg(-\gamma_{n}^{2}\; \int_{0}^{x} \xi(s)\, ds\bigg)\; \dfrac{P_{H,n}(z)}{\|P_{H,n}\|^{2}}\:\bigg(\sum_{h=1}^{H} {P_{h,n}(H_{s})}\bigg). \end{align} $$

The crosswind-integrated concentration $\displaystyle c_{H}^{y}(x,z)$ is obtained by integrating equation (3.9) with respect to y from $-\infty $ to $+\infty $ , which yields

(3.10) $$ \begin{align} c_{H}^{y}(x,z)=Q \sum_{n=0}^{\infty} \exp\!\bigg(-\gamma_{n}^{2}\; \int_{0}^{x} \xi(s)\, ds\bigg)\; \dfrac{P_{H,n}(z)}{\|P_{H,n}\|^{2}}\:\sum_{h=1}^{H} {P_{h,n}(H_{s})}. \end{align} $$

4 Validation and experimental data

The wind speed profile $u(z)$ is parameterized using the model of Deaves and Harris [Reference Deaves and Harris7]. This model is extrapolated using experimentally measured wind speed profiles. The advantage of this profile over other ones is that it extends the accurate representation of the logarithmic law to a small height, and through better accuracy for the logarithmic and power law models to a moderate height,

$$ \begin{align*} u(z)=\dfrac{u_{*}}{k}\bigg[\ln\bigg(\dfrac{z+z_{0}}{z_{0}}\bigg)+ 5.75\bigg(\dfrac{z}{H_{\text{mix}}}\bigg)-1.88 \bigg(\dfrac{z}{H_{\text{mix}}}\bigg)^{2}-1.33\bigg(\dfrac{z}{H_{\text{mix}}}\bigg)^{3} +0.25\bigg(\dfrac{z}{H_{\text{mix}}}\bigg)^{4} \bigg], \end{align*} $$

where $H_{\text {mix}}$ is the PBL height. The modified form of the vertical eddy diffusivity coefficient $K_{z}$ , given in equation (2.6), is adopted from [Reference Degrazia, Velho and Carvalho8], where the functional

$$ \begin{align*} \varphi_{z}(z)=0.22\:H_{\text{mix}}\:w_{*}\:\bigg(\dfrac{z}{H_{\text{mix}}}\;\bigg(1-\dfrac{z }{ H_{\text{mix}}}\bigg)\bigg)^{1/3} \bigg(1-\exp\!\bigg(-\dfrac{4\: z}{H_{\text{mix}}}\bigg)-0.0003\:\exp\!\bigg(\dfrac{8\: z}{H_{\text{mix}}}\bigg)\bigg), \end{align*} $$

where $w_{*}$ is the convective velocity. The integrable correction dimensionless function in equation (2.6) is defined in terms of the along-wind length scale $L_{1}$ as follows [Reference Mooney and Wilson18]:

(4.1) $$ \begin{align} \xi(x)= 1-\exp{\bigg(-\dfrac{x}{L_{1}}\bigg)}. \end{align} $$

The length $L_{1}$ is given in terms of u, $\varphi _{z}$ and $\displaystyle {\sigma _{w}}$ as [Reference Mooney and Wilson18]

$$ \begin{align*} L_{1}=\dfrac{1}{\sigma_{w}^{2}(H_{s})}\:u(H_{s})\: \varphi_{z}(H_{s}), \end{align*} $$

where $\displaystyle {\sigma _{w}}$ is the vertical turbulent intensity. Note that there exist many expressions of $\displaystyle {\sigma _{w}}$ ; we adopt here the expression given by Hanna et al. [Reference Hanna, Briggs, Rayford and Hosker13]:

$$ \begin{align*} {\sigma}_{w}= 0.96 \: w_{*}\bigg(\dfrac{3\: z}{H_{\text{mix}}}+\dfrac{|L|}{H_{\text{mix}}}\bigg)^{1/3}, \end{align*} $$

where L is the Monin–Obukhov length [Reference Obukhov21]. The Monin–Obukhov length is a parameter used in atmospheric dispersion having the dimension of a length, and describes the atmospheric stability states according to their sign, that is, L is negative (positive) for an unstable (stable) situation whereas $|L|>>1$ indicates neutral situation [Reference Obukhov21].

Note that when $\displaystyle {L_{1} \rightarrow 0}$ , the term $\displaystyle {\exp {(-x/L_{1})}}$ becomes negligible in equation (4.1), and $K_{z}$ depends only on z (that is, $\displaystyle {K_{z}\equiv \varphi _{z}}$ in equation (2.6)).

The lateral eddy diffusivity is given by Huang [Reference Huang14]:

(4.2) $$ \begin{align} K_{y}(x,z)= \dfrac{1}{2}\;u(z)\frac{d\: \sigma_{y}^{2}(x)}{dx}, \end{align} $$

where $\sigma _{y}$ is the standard deviation in the crosswind direction (depends only on x). By identifying equations (2.5) and (4.2), we obtain

$$ \begin{align*} \zeta_{y}(x)= \dfrac{1}{2}\;\dfrac{d\: \sigma_{y}^{2}(x)}{dx}. \end{align*} $$

The model presented as equation (3.10) is assessed and validated using data sets obtained from the Copenhagen diffusion and Prairie Grass experiments.

The Copenhagen experiments were realized in the northern part of Copenhagen between 12/09/1978 and 19/07/1979 as described by Gryning et al. [Reference Gryning, Holtslag, Irwin and Sivertsen10] and Gryning and Lyck [Reference Gryning and Lyck11]. They were carried out using the tracer sulphur hexafluoride (SF6), which was realized without buoyancy from a tower at a height of $\displaystyle H_{s} = 115\ m$ , the roughness length was $\displaystyle {z_{0} = 0.6\ m}$ . The remaining required parameters for these experiments are given in [Reference Gryning, Holtslag, Irwin and Sivertsen10, Reference Gryning and Lyck11].

The Prairie Grass experiments were realized in O’Neil, Nebraska between 03/07/1956 and 30/08/1956 and well detailed by Nieuwstadt [Reference Nieuwstadt20] and the data downloaded from http://www.harmo.org/jsirwin (see [Reference Barad1]). They were carried out using the tracer sulphur dioxide ( $SO_{2}$ ), realized without buoyancy from a tower at a height of $\displaystyle H_{s} = 1.5\ m$ , the roughness length was $\displaystyle {z_{0} = 0.006\ m}$ . All necessary data for these experiments are given in [Reference Barad1, Reference Nieuwstadt20].

For the practical implementation of the analytical solution, we have discretized the PBL into $2$ and $4$ sub-layers. For the $2$ -sub-layers case, the discretization is

$$ \begin{align*}\displaystyle{dz_{1}=\frac{{7(H_{\text{mix}}-z_{0})}}{{13}}}\quad \text{ and} \quad \displaystyle{dz_{2} = \frac{{6(H_{\text{mix}}-z_{0})}}{{13}}};\end{align*} $$

and for the $4$ -sub-layers case, the discretization is

$$ \begin{align*} &\displaystyle{dz_{1}=\frac{{3(H_{\text{mix}}-z_{0})}}{{15}}},\quad \displaystyle{dz_{2} = \frac{{3(H_{\text{mix}}-z_{0})}}{{15}}},\quad \displaystyle{dz_{3} = \frac{{7(H_{\text{mix}}-z_{0})}}{{15}}}, \\[6pt] &\text{and} \quad \displaystyle{dz_{4} = \frac{{2(H_{\text{mix}}-z_{0})}}{{15}}}. \end{align*} $$

The quality and performance of models are usually presented by drawing a scatter diagram using predicted and observed values. Figure 1 represents a scatter diagram between observed and predicted crosswind-integrated concentrations for the two formulations of the vertical diffusivity and the number of sub-layers in the PBL (when $H=2$ and $H=4$ ) for the Copenhagen (top panel) and the Prairie Grass (bottom panel) experiments. The figure shows a good agreement between the predicted and observed values which seem to be a good parameterization of the model. An immediate visualization inspection of the overall model performance shows that the observed crosswind-integrated concentration tended to be slightly larger than those predicted. The quantitative evaluation of models is usually made through statistical performance analysis that is presented and widely used in many works. The statistical indices which we will use here are defined as [Reference Chang and Hanna4]: the normalized mean square error $(NMSE)$ , the mean relative square error $(MRSE)$ , the correlation coefficient $(COR)$ , the fractional bias $(FB)$ , the fractional standard deviations $(FS)$ , the geometric mean bias $(MG)$ , the geometric mean variance $(VG)$ and the factor of two $(FAC2)$ . The perfect models would have $\displaystyle {NMSE=MRSE=FB=FS=0}$ and $\displaystyle {COR=MG=VG=FAC2=1}$ [Reference Chang and Hanna4]. The satisfied numerical results given in Figure 1 are re-affirmed from the values of statistical performance measures summarized in Table 1. The calculated values are mostly within the range of acceptable model performance for both formulations of the vertical eddy diffusivity. Furthermore, Table 1 indicates that similar results are obtained compared to other analytical procedures, previously published in the literature [Reference Ema’a, Ben-Bolie, Patrice, Zarma and Pierre9, Reference Guerrero, Pimentel and Skaggs12, Reference Laaouaoucha, Farhane, Essaouini and Souhar15, Reference Moreira, Tirabassi, Vilhenah and Goulart19].

Figure 1 Scatter plot of observed and predicted crosswind-integrated concentrations of the Copenhagen (top panel) and the Prairie Grass experiments (bottom panel). The solid line is a one-to-one line $(y = x)$ and dotted lines correspond to a factor of two (that is, $y=0.5 x$ and $y=2x$ ).

Table 1 Statistical measures for the Copenhagen and Prairie Grass experiments.