**2. On the advection-diffusion approach**

The Eulerian approach is widely used in the field of air pollution studies to model the dispersion properties of the Planetary Boundary Layer (PBL). In this context, the diffusion equation that describes the local mean concentrations *c*¯ = *c*¯(**r**, *t*) at an event point of interest (**r**, *t*)=(*x*, *y*, *z*, *t*) arising from a any contaminant point source, which may be time dependent, can be written as

$$
\partial\_t \overline{\boldsymbol{\varepsilon}} + \mathbf{U} \nabla \overline{\boldsymbol{\varepsilon}} - \nabla^T \mathbf{K} \nabla \overline{\boldsymbol{\varepsilon}} = \mathbf{S} \ . \tag{1}
$$

advection-diffusion equation in each layer is then solved by the Laplace Transform technique. The GIADMT method (Generalized Integral Advection Diffusion Multilayer Technique) ([7]) is a dimensional extension to the previous work, but again assuming the stepwise approximation for the eddy diffusivity coefficient and wind profile. To generalize, a general two-dimensional solution was presented by ([22]). The solving methodology was the Generalized Integral Laplace Transform Technique (GILTT) that is an analytical series solution including the solution of an associate Sturm-Liouville problem, expansion of the pollutant concentration in a series in terms of the attained eigenfunction, replacement of this expansion in the advection-diffusion equation and, finally, taking moments. This procedure leads to a set of differential ordinary equations that is solved analytically by Laplace transform technique. In this work we improve further the solutions of the afore mentioned articles and report on a general analytical solution for the advection-diffusion problem, assuming that eddy diffusivity and wind profiles are arbitrary functions having a continuous dependence on the

Equation (1) is considered valid in the domain (*x*, *y*, *z*) ∈ Γ bounded by 0 < *x* < *Lx*, 0 < *y* < *Ly* (with *Lx* and *Ly* sufficiently large), 0 < *z* < *h* (here *h* is the boundary layer height) and

Instead of specifying the source term as an inhomogeneity of the partial differential equation, we consider a point source located at an edge of the domain, so that the source position **r***<sup>S</sup>* = (0, *y*0, *HS*) is located at the boundary of the domain **r***<sup>S</sup>* ∈ *δ*Γ. Note, that in cases where the source is located in the domain, one still may divide the whole domain in sub-domains, where the source lies on the boundary of the sub-domains which can be solved for each sub-domain separately. Moreover, a set of different sources may be implemented as a superposition of independent problems. Since the source term location is on the boundary, in the domain this term is zero everywhere (*S*(**r**) ≡ 0 for **r** ∈ Γ\*δ*Γ), so that the source influence may be cast in form of a condition rather than a source term of the equation. The source condition for a time

> *S* =

*<sup>c</sup>*¯(*t*, *<sup>x</sup>*, *<sup>y</sup>*, *<sup>z</sup>*) = *<sup>t</sup>*

lim **<sup>Σ</sup>**<sup>ˆ</sup> *<sup>d</sup>***Σ**→<sup>0</sup>

<sup>0</sup> <sup>=</sup>

where ω*<sup>S</sup>* represents a flux across a closed surface that includes the source and is proportional to the source strength. Instead of considering an explicit source term, we implement the solution as a superposition of an infinite number of solutions with instantaneous source represented in an initial condition. The solution for a time dependent source assumes the

0

**K**∇*c*¯|(0,0,0) = **K**∇*c*¯|(*Lx*,*Ly*,*h*) = **0** (3)

*c*¯(*x*, *y*, *z*, 0) = 0 . (4)

On an Analytical Model for the Radioactive Contaminant Release in the Atmosphere from Nuclear Power Plants

221

ω*<sup>S</sup> d***Σ** (5)

*<sup>c</sup>*¯˙(*<sup>t</sup>* − *<sup>τ</sup>*, *<sup>x</sup>*, *<sup>y</sup>*, *<sup>z</sup>*) *<sup>d</sup><sup>τ</sup>* (6)

ω*<sup>S</sup> d***Σ** = *Qδ*(*x*)*δ*(*y* − *y*0)*δ*(*z* − *HS*) (7)

vertical and longitudinal spatial variables.

dependent contamination reads then

with instantaneous initial condition

*c*¯˙(0, *x*, *y*, *z*) = *c*¯˙

following form

subject to the following boundary and initial conditions,

Here **U** = (*u*¯, *v*¯, *w*¯)*<sup>T</sup>* is the vector field of the mean wind velocity, the diagonal matrix **K** = diag(*Kx*, *Ky*, *Kz*) contains the eddy diffusivities and *S* is a source term, to be determined according to the scenario of interest. In equation (1) we tacitly related the turbulent fluxes **U**�*c*� to the gradient of the mean concentration by means of eddy diffusivity (K-theory)

$$\mathbf{U}\overline{\mathbf{U}^{\prime}} = -\mathbf{K}\nabla\overline{\varepsilon} \tag{2}$$

The simplicity of the K-theory has led to the widespread use of this theory as mathematical basis for simulating air pollution phenomena. However, the K-closure has its intrinsic limits: it works well when the dimension of dispersed material is much larger than the size of turbulent eddies involved in the diffusion process. Another crucial point is that the down-gradient transport hypothesis is inconsistent with observed features of turbulent diffusion in the upper portion of the mixed layer ([9]). Despite these well known limits, the K-closure is largely used in several atmospheric conditions because it describes the diffusive transport in an Eulerian framework where almost all measurements are easily cast into an Eulerian form, it produces results that agree with experimental data as well as any other more complex model, and it is not computationally expensive as higher order closures usually are.

For a time dependent regime considered in the present work, we assume that the associated advection-diffusion equation adequately describes a dispersion process of radioactive material. From applications of the approach to tracer dispersion data we saw that our analytical approach does not only yield a solution for the three dimensional advection-diffusion equation but predicts tracer concentrations closer to observed values compared to other approaches from the literature, which is also manifest in better statistical coefficients.

Approaches to the advection-diffusion problem are not new in the literature, that are either based on numerical schemes, stochastic simulations or (semi-)analytical methods as shown in a selection of articles ([12, 23, 26, 29, 32]). Note, that in these works all solutions are valid for scenarios with strong restrictions with respect to their specific wind and vertical eddy diffusivity profiles. A more general approach, the ADMM (Advection Diffusion Multilayer Method) approach solves the two-dimensional advection-diffusion equation with variable wind profile and eddy diffusivity coefficient ([21]). The main idea here relies on the discretisation of the atmospheric boundary layer in a multi-shell domain, assuming in each layer that eddy diffusivity and wind profile take averaged values. The resulting advection-diffusion equation in each layer is then solved by the Laplace Transform technique. The GIADMT method (Generalized Integral Advection Diffusion Multilayer Technique) ([7]) is a dimensional extension to the previous work, but again assuming the stepwise approximation for the eddy diffusivity coefficient and wind profile. To generalize, a general two-dimensional solution was presented by ([22]). The solving methodology was the Generalized Integral Laplace Transform Technique (GILTT) that is an analytical series solution including the solution of an associate Sturm-Liouville problem, expansion of the pollutant concentration in a series in terms of the attained eigenfunction, replacement of this expansion in the advection-diffusion equation and, finally, taking moments. This procedure leads to a set of differential ordinary equations that is solved analytically by Laplace transform technique. In this work we improve further the solutions of the afore mentioned articles and report on a general analytical solution for the advection-diffusion problem, assuming that eddy diffusivity and wind profiles are arbitrary functions having a continuous dependence on the vertical and longitudinal spatial variables.

Equation (1) is considered valid in the domain (*x*, *y*, *z*) ∈ Γ bounded by 0 < *x* < *Lx*, 0 < *y* < *Ly* (with *Lx* and *Ly* sufficiently large), 0 < *z* < *h* (here *h* is the boundary layer height) and subject to the following boundary and initial conditions,

$$\mathbf{K} \nabla \bar{c}|\_{\left(0,0,0\right)} = \mathbf{K} \nabla \bar{c}|\_{\left(L\_x L\_y h\right)} = \mathbf{0} \tag{3}$$

$$
\bar{\varepsilon}(x, y, z, 0) = 0 \,. \tag{4}
$$

Instead of specifying the source term as an inhomogeneity of the partial differential equation, we consider a point source located at an edge of the domain, so that the source position **r***<sup>S</sup>* = (0, *y*0, *HS*) is located at the boundary of the domain **r***<sup>S</sup>* ∈ *δ*Γ. Note, that in cases where the source is located in the domain, one still may divide the whole domain in sub-domains, where the source lies on the boundary of the sub-domains which can be solved for each sub-domain separately. Moreover, a set of different sources may be implemented as a superposition of independent problems. Since the source term location is on the boundary, in the domain this term is zero everywhere (*S*(**r**) ≡ 0 for **r** ∈ Γ\*δ*Γ), so that the source influence may be cast in form of a condition rather than a source term of the equation. The source condition for a time dependent contamination reads then

$$\mathcal{S} = \oint \omega\_{\mathcal{S}} \, d\Sigma \tag{5}$$

where ω*<sup>S</sup>* represents a flux across a closed surface that includes the source and is proportional to the source strength. Instead of considering an explicit source term, we implement the solution as a superposition of an infinite number of solutions with instantaneous source represented in an initial condition. The solution for a time dependent source assumes the following form

$$
\vec{\varepsilon}(t, \mathbf{x}, y, z) = \int\_0^t \dot{\vec{\varepsilon}}(t - \tau, \mathbf{x}, y, z) \, d\tau \tag{6}
$$

with instantaneous initial condition

2 Will-be-set-by-IN-TECH

and Forecasting" (WRF). The forcasting system contains a three dimensional data assimilation system and is suitable for applications from the meso- down to the micro-scale. The second step plays the role of simulating the dispersion process in a micro-scale, i.e. in the environment

The Eulerian approach is widely used in the field of air pollution studies to model the dispersion properties of the Planetary Boundary Layer (PBL). In this context, the diffusion equation that describes the local mean concentrations *c*¯ = *c*¯(**r**, *t*) at an event point of interest (**r**, *t*)=(*x*, *y*, *z*, *t*) arising from a any contaminant point source, which may be time dependent,

Here **U** = (*u*¯, *v*¯, *w*¯)*<sup>T</sup>* is the vector field of the mean wind velocity, the diagonal matrix **K** = diag(*Kx*, *Ky*, *Kz*) contains the eddy diffusivities and *S* is a source term, to be determined according to the scenario of interest. In equation (1) we tacitly related the turbulent fluxes **U**�*c*�

The simplicity of the K-theory has led to the widespread use of this theory as mathematical basis for simulating air pollution phenomena. However, the K-closure has its intrinsic limits: it works well when the dimension of dispersed material is much larger than the size of turbulent eddies involved in the diffusion process. Another crucial point is that the down-gradient transport hypothesis is inconsistent with observed features of turbulent diffusion in the upper portion of the mixed layer ([9]). Despite these well known limits, the K-closure is largely used in several atmospheric conditions because it describes the diffusive transport in an Eulerian framework where almost all measurements are easily cast into an Eulerian form, it produces results that agree with experimental data as well as any other more complex model, and it is

For a time dependent regime considered in the present work, we assume that the associated advection-diffusion equation adequately describes a dispersion process of radioactive material. From applications of the approach to tracer dispersion data we saw that our analytical approach does not only yield a solution for the three dimensional advection-diffusion equation but predicts tracer concentrations closer to observed values compared to other approaches from the literature, which is also manifest in better statistical

Approaches to the advection-diffusion problem are not new in the literature, that are either based on numerical schemes, stochastic simulations or (semi-)analytical methods as shown in a selection of articles ([12, 23, 26, 29, 32]). Note, that in these works all solutions are valid for scenarios with strong restrictions with respect to their specific wind and vertical eddy diffusivity profiles. A more general approach, the ADMM (Advection Diffusion Multilayer Method) approach solves the two-dimensional advection-diffusion equation with variable wind profile and eddy diffusivity coefficient ([21]). The main idea here relies on the discretisation of the atmospheric boundary layer in a multi-shell domain, assuming in each layer that eddy diffusivity and wind profile take averaged values. The resulting

to the gradient of the mean concentration by means of eddy diffusivity (K-theory)

not computationally expensive as higher order closures usually are.

*<sup>∂</sup>tc*¯ <sup>+</sup> **<sup>U</sup>**∇*c*¯ <sup>−</sup> <sup>∇</sup>*T***K**∇*c*¯ <sup>=</sup> *<sup>S</sup>* . (1)

**U**�*c*� = −**K**∇*c*¯ (2)

within a radius of several tenth kilometers.

can be written as

coefficients.

**2. On the advection-diffusion approach**

$$\dot{\hat{\boldsymbol{\varepsilon}}}(0,\mathbf{x},\mathbf{y},\mathbf{z}) = \dot{\hat{\boldsymbol{\varepsilon}}}\_{0} = \lim\_{\oint \hat{\boldsymbol{\Sigma}} \, d\mathbf{\hat{z}} \to 0} \oint \boldsymbol{\omega}\_{\rm S} \, d\mathbf{\hat{z}} = \mathbf{Q} \delta(\mathbf{x}) \delta(\boldsymbol{y} - \boldsymbol{y}\_{0}) \delta(\boldsymbol{z} - \boldsymbol{H}\_{\rm S}) \tag{7}$$

where *Q* is the emission rate, *HS* the height of the source, *δ*(*x*) represents the Cartesian Dirac delta functional and **Σ**ˆ is a unit vector.

### **3. A closed form solution**

In this section we first introduce the general formalism to solve a general problem and subsequently reduce the problem to a more specific one, that is solved and compared to experimental findings.

#### **3.1. The general procedure**

In order to solve the problem (1) we reduce the dimensionality by one and thus cast the problem into a form already solved in reference [22]. To this end we apply the integral transform technique in the *y* variable, and expand the pollutant concentration as

$$\overline{\mathbf{c}}(\mathbf{x}, y, z, t) = \mathbf{R}^T(\mathbf{x}, z, t)\mathbf{Y}(y), \tag{8}$$

Here, **<sup>B</sup>**<sup>0</sup> <sup>=</sup> *Ly*

reference [22]:

and *∂xKx*.

where **Λ** = diag(*λ*<sup>2</sup>

<sup>1</sup>, *<sup>λ</sup>*<sup>2</sup>

which is equivalent to the problem

by virtue of **B** being a diagonal matrix.

(*∂t***R***T*)**B** + **U**¯

**3.2. A specific case for application**

<sup>∇</sup>2**R***T***<sup>B</sup>** <sup>+</sup> **<sup>R</sup>***T***<sup>Z</sup>**

<sup>2</sup> **<sup>I</sup>**, where **<sup>I</sup>** is the identity, the elements (**Z**)*mn* <sup>=</sup> <sup>2</sup>

In order to discuss a specific case we introduce a convention and consider the average wind velocity **U**¯ = (*u*¯, 0, 0)*<sup>T</sup>* aligned with the *x*-axis. Since the variation of the average wind velocity is slow compared to the time intervals for which the meteorological data are extracted from WRF, we superimpose the solution after rotation in the *x* − *y*-plane in order to transform every instantaneous solution into the same coordinate frame, i.e. the coordinate frame for *t* = 0. By comparison of physically meaningful cases, one finds for the operator norm ||*∂xKx∂x*|| << |*u*¯|, which can be understood intuitively because eddy diffusion is observable predominantly perpendicular to the mean wind propagation. As a consequence we neglect the terms with *Kx*

The principal aspect of interest in pollution dispersion is the vertical concentration profile, that responds strongly to the atmospheric boundary layer stratification, so that the simplified eddy diffusivity **K** → **K**<sup>1</sup> = diag(0, *Ky*, *Kz*) depends in leading order approximation only on

**<sup>Ω</sup>**<sup>1</sup> <sup>→</sup> (*∂zKz*)(*∂z***R***T*)**<sup>B</sup>** ,

<sup>2</sup>,...). The simplified equation system to be solved is then,

the vertical coordinate **K**<sup>1</sup> = **K**1(*z*). For this specific case the integrals **Ω***<sup>i</sup>* reduce to

**<sup>Ω</sup>**<sup>2</sup> <sup>→</sup> *Kz*(*∂*<sup>2</sup>

**T**<sup>1</sup> → **0** ,

*∂t***R***T***B** + *u*¯*∂x***R***T***B** = (*∂zKz*)*∂z***R***T***B** + *Kz∂*<sup>2</sup>

*∂t***R** + *u*¯*∂x***R** = (*∂zKz*)*∂z***R** + *Kz∂*<sup>2</sup>

The specific form of the eddy diffusivity determines now whether the problem is a linear or non-linear one. In the linear case the **K** is assumed to be independent of *c*¯, whereas in more realistic cases, even if stationary, **K** may depend on the contaminant concentration and thus renders the problem non-linear. However, until now now specific law is known

Kronecker symbol and *j* = (*m* + *n*)mod2 is the remainder of an integer division (i.e. is one for *m* + *n* odd and zero else). Note, that the integrals **Ω***<sup>i</sup>* and **T***<sup>i</sup>* depend on the specific form of the eddy diffusivity **K**. The integrals (11) are general, but for practical purposes and for application to a case study we truncate the eigenfunction space and consider *M* components in **R** and **Y** only, though continue using the general nomenclature that remains valid. The obtained matrix equation determines now together with initial and boundary condition uniquely the components *Ri* for *i* = 1, . . . , *M* following the procedure introduced in

<sup>1</sup>−*n*2/*m*<sup>2</sup> *<sup>δ</sup>*1,*<sup>j</sup>* with *<sup>δ</sup>i*,*<sup>j</sup>* the

223

= **Ω**1(**R**) + **Ω**2(**R**) + **R***T*(**T**<sup>1</sup> + **T**2) (12)

On an Analytical Model for the Radioactive Contaminant Release in the Atmosphere from Nuclear Power Plants

*<sup>z</sup>***R***T*)**<sup>B</sup>** , (13)

*<sup>z</sup>***R***T***<sup>B</sup>** <sup>−</sup> *Ky***R***T***Λ<sup>B</sup>** (15)

*<sup>z</sup>***R** − *Ky***ΛR** (16)

**T**<sup>2</sup> → −*Ky***ΛB** , (14)

where **R** = (*R*1, *R*2,...)*<sup>T</sup>* and **Y** = (*Y*1,*Y*2,...)*<sup>T</sup>* is a vector in the space of orthogonal eigenfunctions, given by *Ym*(*y*) = cos(*λmy*) with eigenvalues *λ<sup>m</sup>* = *m <sup>π</sup> Ly* for *m* = 0, 1, 2, . . .. For convenience we introduce some shorthand notations, <sup>∇</sup><sup>2</sup> = (*∂x*, 0, *<sup>∂</sup>y*)*<sup>T</sup>* and <sup>ˆ</sup> *∂<sup>y</sup>* = (0, *∂y*, 0)*T*, so that equation (1) reads now,

$$\begin{aligned} \left(\partial\_l \mathbf{R}^T\right) \mathbf{Y} + \bar{\mathbf{U}} \left(\nabla\_2 \mathbf{R}^T \mathbf{Y} + \mathbf{R}^T \hat{\partial}\_{\mathcal{Y}} \mathbf{Y}\right) &= \left(\nabla^T \mathbf{K} + (\mathbf{K} \nabla)^T\right) \left(\nabla\_2 \mathbf{R}^T \mathbf{Y} + \mathbf{R}^T \hat{\partial}\_{\mathcal{Y}} \mathbf{Y}\right) \\ &= \left(\nabla\_2^T \mathbf{K} + (\mathbf{K} \nabla\_2)^T\right) \left(\nabla\_2 \mathbf{R}^T \mathbf{Y}\right) + \left(\hat{\partial}\_{\mathcal{Y}}^T \mathbf{K} + (\mathbf{K} \hat{\partial}\_{\mathcal{Y}})^T\right) \left(\mathbf{R}^T \hat{\partial}\_{\mathcal{Y}} \mathbf{Y}\right) .\end{aligned} \tag{9}$$

Upon application of the integral operator

$$\int\_{0}^{L\_{\text{y}}} dy \, \mathbf{Y}[\mathbf{F}] = \int\_{0}^{L\_{\text{y}}} \mathbf{F}^{T} \wedge \mathbf{Y} \, dy \tag{10}$$

here **F** is an arbitrary function and ∧ signifies the dyadic product operator, and making use of orthogonality renders equation (9) a matrix equation. The appearing integral terms are

**B**<sup>0</sup> = *Ly* 0 *dy***Y**[**Y**] = *Ly* 0 **<sup>Y</sup>***<sup>T</sup>* <sup>∧</sup> **<sup>Y</sup>** *dy* , **Z** = *Ly* 0 *dy***Y**[ˆ *<sup>∂</sup>y***Y**] = *Ly* 0 ˆ *<sup>∂</sup>y***Y***<sup>T</sup>* <sup>∧</sup> **<sup>Y</sup>** *dy* , **Ω**<sup>1</sup> = *Ly* 0 *dy***Y**[(∇*<sup>T</sup>* <sup>2</sup> **<sup>K</sup>**)(∇2**R***T***Y**)] = *Ly* 0 (∇*<sup>T</sup>* <sup>2</sup> **<sup>K</sup>**)(∇2**R***T***Y**) *T* ∧ **Y** *dy* , (11) **Ω**<sup>2</sup> = *Ly* 0 *dy***Y**[(**K**∇2)*T*(∇2**R***T***Y**)] = *Ly* 0 (**K**∇2)*T*(∇2**R***T***Y**) ∧ **Y** *dy* , **T**<sup>1</sup> = *Ly* 0 *dy***Y**[((ˆ *∂T <sup>y</sup>* **<sup>K</sup>**)(<sup>ˆ</sup> *<sup>∂</sup>y***Y**)] = *Ly* 0 ((ˆ *∂T <sup>y</sup>* **<sup>K</sup>**)(<sup>ˆ</sup> *∂y***Y**) *T* ∧ **Y** *dy* , **T**<sup>2</sup> = *Ly* 0 *dy***Y**[(**K**ˆ *∂y*)*T*(ˆ *<sup>∂</sup>y***Y**)] = *Ly* 0 (**K**ˆ *∂y*)*T*(ˆ *∂y***Y**) *T* ∧ **Y** *dy* .

Here, **<sup>B</sup>**<sup>0</sup> <sup>=</sup> *Ly* <sup>2</sup> **<sup>I</sup>**, where **<sup>I</sup>** is the identity, the elements (**Z**)*mn* <sup>=</sup> <sup>2</sup> <sup>1</sup>−*n*2/*m*<sup>2</sup> *<sup>δ</sup>*1,*<sup>j</sup>* with *<sup>δ</sup>i*,*<sup>j</sup>* the Kronecker symbol and *j* = (*m* + *n*)mod2 is the remainder of an integer division (i.e. is one for *m* + *n* odd and zero else). Note, that the integrals **Ω***<sup>i</sup>* and **T***<sup>i</sup>* depend on the specific form of the eddy diffusivity **K**. The integrals (11) are general, but for practical purposes and for application to a case study we truncate the eigenfunction space and consider *M* components in **R** and **Y** only, though continue using the general nomenclature that remains valid. The obtained matrix equation determines now together with initial and boundary condition uniquely the components *Ri* for *i* = 1, . . . , *M* following the procedure introduced in reference [22]:

$$(\partial\_l \mathbf{R}^T) \mathbf{B} + \bar{\mathbf{U}} \left(\nabla\_2 \mathbf{R}^T \mathbf{B} + \mathbf{R}^T \mathbf{Z}\right) = \mathbf{D}\_1(\mathbf{R}) + \mathbf{D}\_2(\mathbf{R}) + \mathbf{R}^T(\mathbf{T}\_1 + \mathbf{T}\_2) \tag{12}$$

#### **3.2. A specific case for application**

4 Will-be-set-by-IN-TECH

where *Q* is the emission rate, *HS* the height of the source, *δ*(*x*) represents the Cartesian Dirac

In this section we first introduce the general formalism to solve a general problem and subsequently reduce the problem to a more specific one, that is solved and compared to

In order to solve the problem (1) we reduce the dimensionality by one and thus cast the problem into a form already solved in reference [22]. To this end we apply the integral

where **R** = (*R*1, *R*2,...)*<sup>T</sup>* and **Y** = (*Y*1,*Y*2,...)*<sup>T</sup>* is a vector in the space of orthogonal

<sup>∇</sup>*T***<sup>K</sup>** + (**K**∇)*<sup>T</sup>*

 ˆ *∂T*

(∇2**R***T***Y**) +

 *Ly* 0

here **F** is an arbitrary function and ∧ signifies the dyadic product operator, and making use of orthogonality renders equation (9) a matrix equation. The appearing integral terms are

*<sup>∂</sup>y***Y***<sup>T</sup>* <sup>∧</sup> **<sup>Y</sup>** *dy* ,

 *Ly* 0

 *Ly* 0

 *Ly* 0

> *Ly* 0

 ((ˆ *∂T <sup>y</sup>* **<sup>K</sup>**)(<sup>ˆ</sup> *∂y***Y**) *T*

 (**K**ˆ *∂y*)*T*(ˆ *∂y***Y**) *T*

 (∇*<sup>T</sup>*

For convenience we introduce some shorthand notations, <sup>∇</sup><sup>2</sup> = (*∂x*, 0, *<sup>∂</sup>y*)*<sup>T</sup>* and <sup>ˆ</sup>

*dy***Y**[**F**] =

**<sup>Y</sup>***<sup>T</sup>* <sup>∧</sup> **<sup>Y</sup>** *dy* ,

*c*¯(*x*, *y*, *z*, *t*) = **R***T*(*x*, *z*, *t*)**Y**(*y*), (8)

<sup>2</sup> **<sup>K</sup>**)(∇2**R***T***Y**)

(**K**∇2)*T*(∇2**R***T***Y**)

*<sup>y</sup>* **<sup>K</sup>** + (**K**<sup>ˆ</sup>

<sup>∇</sup>2**R***T***<sup>Y</sup>** <sup>+</sup> **<sup>R</sup>***<sup>T</sup>* <sup>ˆ</sup>

*∂y*)*<sup>T</sup>* (**R***<sup>T</sup>* ˆ

**<sup>F</sup>***<sup>T</sup>* <sup>∧</sup> **<sup>Y</sup>** *dy* (10)

*T*

∧ **Y** *dy* ,

∧ **Y** *dy* .

∧ **Y** *dy* ,

*Ly* for *m* = 0, 1, 2, . . ..

*∂y***Y**  *∂<sup>y</sup>* =

*∂y***Y**) . (9)

∧ **Y** *dy* , (11)

transform technique in the *y* variable, and expand the pollutant concentration as

eigenfunctions, given by *Ym*(*y*) = cos(*λmy*) with eigenvalues *λ<sup>m</sup>* = *m <sup>π</sup>*

*∂y***Y** = 

<sup>2</sup> **<sup>K</sup>** + (**K**∇2)*<sup>T</sup>*

 *Ly* 0

 *Ly* 0

> *Ly* 0 ˆ

*dy***Y**[(**K**∇2)*T*(∇2**R***T***Y**)] =

*∂y*)*T*(ˆ

<sup>2</sup> **<sup>K</sup>**)(∇2**R***T***Y**)] =

*∂y***Y**)] =

*∂y***Y**)] =

delta functional and **Σ**ˆ is a unit vector.

**3. A closed form solution**

**3.1. The general procedure**

(0, *∂y*, 0)*T*, so that equation (1) reads now,

Upon application of the integral operator

 *Ly* 0

 *Ly* 0

 *Ly* 0

 *Ly* 0

 *Ly* 0

 *Ly* 0

*dy***Y**[**Y**] =

*dy***Y**[(∇*<sup>T</sup>*

*dy***Y**[((ˆ *∂T <sup>y</sup>* **<sup>K</sup>**)(<sup>ˆ</sup>

*dy***Y**[(**K**ˆ

*∂y***Y**] =

*dy***Y**[ˆ

= ∇*T*

<sup>∇</sup>2**R***T***<sup>Y</sup>** <sup>+</sup> **<sup>R</sup>***<sup>T</sup>* <sup>ˆ</sup>

(*∂t***R***T*)**Y** + **U**¯

**B**<sup>0</sup> =

**Z** =

**Ω**<sup>1</sup> =

**Ω**<sup>2</sup> =

**T**<sup>1</sup> =

**T**<sup>2</sup> =

experimental findings.

In order to discuss a specific case we introduce a convention and consider the average wind velocity **U**¯ = (*u*¯, 0, 0)*<sup>T</sup>* aligned with the *x*-axis. Since the variation of the average wind velocity is slow compared to the time intervals for which the meteorological data are extracted from WRF, we superimpose the solution after rotation in the *x* − *y*-plane in order to transform every instantaneous solution into the same coordinate frame, i.e. the coordinate frame for *t* = 0. By comparison of physically meaningful cases, one finds for the operator norm ||*∂xKx∂x*|| << |*u*¯|, which can be understood intuitively because eddy diffusion is observable predominantly perpendicular to the mean wind propagation. As a consequence we neglect the terms with *Kx* and *∂xKx*.

The principal aspect of interest in pollution dispersion is the vertical concentration profile, that responds strongly to the atmospheric boundary layer stratification, so that the simplified eddy diffusivity **K** → **K**<sup>1</sup> = diag(0, *Ky*, *Kz*) depends in leading order approximation only on the vertical coordinate **K**<sup>1</sup> = **K**1(*z*). For this specific case the integrals **Ω***<sup>i</sup>* reduce to

$$\begin{aligned} \mathbf{^0 \mathbf{D}\_1 \rightarrow (\partial\_{\overline{z}} K\_{\overline{z}}) (\partial\_{\overline{z}} \mathbf{R}^T) \mathbf{B} \ , \\ \mathbf{^0 \mathbf{D}\_2 \rightarrow K\_{\overline{z}} (\partial\_{\overline{z}}^2 \mathbf{R}^T) \mathbf{B} \ , \end{aligned} \tag{13}$$

$$\begin{aligned} \mathbf{T}\_1 &\rightarrow \mathbf{0} \\ \mathbf{T}\_2 &\rightarrow -\mathbf{K}\_\mathcal{Y} \mathbf{A} \mathbf{B} \end{aligned} \tag{14}$$

where **Λ** = diag(*λ*<sup>2</sup> <sup>1</sup>, *<sup>λ</sup>*<sup>2</sup> <sup>2</sup>,...). The simplified equation system to be solved is then,

$$
\partial\_l \mathbf{R}^T \mathbf{B} + \vec{u} \partial\_\mathbf{x} \mathbf{R}^T \mathbf{B} = (\partial\_z K\_z) \partial\_{\bar{z}} \mathbf{R}^T \mathbf{B} + K\_{\bar{z}} \partial\_{\bar{z}}^2 \mathbf{R}^T \mathbf{B} - K\_y \mathbf{R}^T \mathbf{A} \mathbf{B} \tag{15}
$$

which is equivalent to the problem

$$
\hat{\boldsymbol{\partial}}\_{l}\mathbf{R} + \bar{\boldsymbol{u}}\hat{\boldsymbol{\partial}}\_{\mathbf{x}}\mathbf{R} = (\hat{\boldsymbol{\partial}}\_{2}\boldsymbol{K}\_{2})\hat{\boldsymbol{\partial}}\_{2}\mathbf{R} + \boldsymbol{K}\_{2}\hat{\boldsymbol{\partial}}\_{2}^{2}\mathbf{R} - \boldsymbol{K}\_{\mathcal{Y}}\boldsymbol{\Lambda}\mathbf{R} \tag{16}
$$

by virtue of **B** being a diagonal matrix.

The specific form of the eddy diffusivity determines now whether the problem is a linear or non-linear one. In the linear case the **K** is assumed to be independent of *c*¯, whereas in more realistic cases, even if stationary, **K** may depend on the contaminant concentration and thus renders the problem non-linear. However, until now now specific law is known

that links the eddy diffusivity to the concentration so that we hide this dependence using a phenomenologically motivated expression for **K** which leaves us with a partial differential equation system in linear form, although the original phenomenon is non-linear. In the example below we demonstrate the closed form procedure for a problem with explicit time dependence, which is novel in the literature.

The solution is generated making use of the decomposition method ([1–3]) which was originally proposed to solve non-linear partial differential equations, followed by the Laplace transform that renders the problem a pseudo-stationary one. Further we rewrite the vertical diffusivity as a time average term *K*¯ *<sup>z</sup>*(*z*) plus a term representing the variations *κz*(*z*, *t*) around the average for the time interval of the measurement *Kz*(*x*, *z*, *t*) = *K*¯ *<sup>z</sup>*(*z*) + *κz*(*z*, *t*) and use the asymptotic form of *Ky*, which is then explored to set-up the structure of the equation that defines the recursive decomposition scheme:

$$
\partial\_t \mathbf{R} + \bar{u} \partial\_{\mathbf{x}} \mathbf{R} - \partial\_{\mathbf{z}} \left( \bar{K}\_{\mathbf{z}} \partial\_{\mathbf{z}} \mathbf{R} \right) + K\_{\mathcal{Y}} \mathbf{A} \mathbf{R} = \partial\_{\mathbf{z}} \left( \kappa\_{\mathbf{z}} \partial\_{\mathbf{z}} \mathbf{R} \right) \tag{17}
$$

The *x* and *z* dependence may be separated using the same reasoning as already introduced in

where **C** = (*ζ*1(*z*), *ζ*2(*z*),...)*<sup>T</sup>* are a set of orthogonal eigenfunctions, given by *ζi*(*z*) =

Replacing equation (21) in equation (20) and using the afore introduced projector (10) now for

<sup>1</sup> **B**2. The entries of matrices **B**<sup>1</sup> and **B**<sup>2</sup> are

*i h* 0

*<sup>i</sup> Ky h* 0

*dy***Y**[*δ*(*<sup>y</sup>* <sup>−</sup> *<sup>y</sup>*0)] = *<sup>Q</sup>***B**−<sup>1</sup>

Following the reasoning of [22] we solve (22) applying Laplace transform and diagonalisation

where **P**˜(*s*,*r*) denotes the Laplace Transform of **P**(*x*,*r*). Here **X**(−1) is the (inverse) matrix of

(*s***I** + **D**)*ii* = *s* + *di*. After performing the Laplace transform inversion of equation (24), we

where **G**(*x*,*r*) is the diagonal matrix with components (**G**)*ii* = *e*−*dix*. Further the still

The analytical time dependence for the recursion initialisation (20) is obtained upon applying

 *γ*+*i*∞ *γ*−*i*∞

To overcome the drawback of evaluating the line integral appearing in the above solution, we perform the calculation of this integral by the Gaussian quadrature scheme, which is exact if

> *t* **<sup>a</sup>***<sup>T</sup>*

<sup>0</sup> *dz***C**[**F**] = *<sup>h</sup>*

*uζi*(*z*)*ζj*(*z*) *dz*

*<sup>∂</sup>zKz∂zζi*(*z*)*ζj*(*z*) *dz* <sup>−</sup> *<sup>γ</sup>*<sup>2</sup>

*<sup>ζ</sup>i*(*z*)*ζj*(*z*) *dz* <sup>−</sup> *<sup>λ</sup>*<sup>2</sup>

**R**˜ <sup>0</sup> = **PC** (21)

On an Analytical Model for the Radioactive Contaminant Release in the Atmosphere from Nuclear Power Plants

*∂x***P** + **HP** = **0** , (22)

*Kzζi*(*z*)*ζj*(*z*) *dz*

*ζi*(*z*)*ζj*(*z*) *dz* .

**P**˜(*s*,*r*) = **X**(*s***I** + **D**)−1**X**−1**P**(0,*r*) (24)

<sup>1</sup> **B**<sup>2</sup> with diagonal eigenvalue matrix **D** and the entries of matrix

**P**(*x*,*r*)**C**(*z*)*e*

**pR**0(*x*, *z*,

*<sup>r</sup>* variable

**p** *t* ) 

**P**(*x*,*r*) = **XG**(*x*,*r*)**X**−1**Ω** , (25)

<sup>0</sup> **<sup>F</sup>***<sup>T</sup>* <sup>∧</sup> **<sup>C</sup>** *dz* yields a first order differential

225

<sup>1</sup> (**C**(*HS*) ∧ **1**) (**1** ∧ **Y**(*y*0)) (23)

*rt dr* . (26)

, (27)

(8). To this end we pose the solution of problem (20) in the form:

cos(*γlz*), and *γ<sup>i</sup>* = *iπ*/*h* (for *i* = 0, 1, 2, . . .) are the set of eigenvalues.

 *h* 0

 *h* 0

−*r h* 0

A similar procedure leads to the source condition for (22):

*dz***C**[*δ*(*<sup>z</sup>* <sup>−</sup> *HS*)]

unknown arbitrary constant matrix is given by **Ω** = **X**−1**P**(0,*r*).

**<sup>R</sup>**0(*x*, *<sup>z</sup>*, *<sup>t</sup>*) = <sup>1</sup>

the integrand is a polynomial of degree 2*<sup>M</sup>* <sup>−</sup> 1 in the <sup>1</sup>

2*πi*

**<sup>R</sup>**0(*x*, *<sup>z</sup>*, *<sup>t</sup>*) = <sup>1</sup>

the *z* dependent degrees of freedom *<sup>h</sup>*

(**B**1)*i*,*<sup>j</sup>* = −

(**B**2)*i*,*<sup>j</sup>* =

of the matrix **H** = **XDX**−<sup>1</sup> which results in

the inverse Laplace transform definition

where **P** = **P**(*x*,*r*) and **H** = **B**−<sup>1</sup>

equation system:

**P**(0,*r*) = *Q***B**−<sup>1</sup>

come out with

1 

the eigenvectors of matrix **B**−<sup>1</sup>

The function **R** = ∑*<sup>j</sup>* **R***<sup>j</sup>* = **1***T***R**(*c*) is now decomposed into contributions to be determined by recursion. For convenience we introduced the one-vector **1** = (1, 1, . . .)*<sup>T</sup>* and inflate the vector **R** to a vector with each element being itself a vector **R***j*. Upon inserting the expansion in equation (17) one may regroup terms that obey the recursive equations and starts with the time averaged solution for *Kz*:

$$
\partial\_t \mathbf{R}\_0 + \bar{u} \partial\_{\bar{x}} \mathbf{R}\_0 - \partial\_{\bar{z}} \left( \bar{K}\_{\bar{z}} \partial\_{\bar{z}} \mathbf{R}\_0 \right) + K\_{\bar{y}} \mathbf{A} \mathbf{R}\_0 = 0 \tag{18}
$$

The extension to the closed form recursion is then given by

$$
\partial\_l \mathbf{R}\_j + \mathbf{\bar{u}} \partial\_{\mathbf{x}} \mathbf{R}\_j - \partial\_z \left( \bar{K}\_2 \partial\_{\mathbf{z}} \mathbf{R}\_j \right) + K\_y \mathbf{\bar{R}}\_j = \partial\_z \left( \kappa\_z \partial\_{\mathbf{z}} \mathbf{R}\_{j-1} \right) \,. \tag{19}
$$

From the construction of the recursion equation system it is evident that other schemes are possible. The specific choice made here allows us to solve the recursion initialisation using the procedure described in reference [22], where a stationary **K** was assumed. For this reason the time dependence enters as a known source term from the first recursion step on.

#### **3.3. Recursion initialisation**

The boundary conditions are now used to uniquely determine the solution. In our scheme the initialisation solution that contains **R**<sup>0</sup> satisfies the boundary conditions (equations (3)) while the remaining equations satisfy homogeneous boundary conditions. Once the set of problems (19) is solved by the GILTT method, the solution of problem (1) is well determined. It is important to consider that we may control the accuracy of the results by a proper choice of the number of terms in the solution series.

In reference [22] a two dimensional problem with advection in the *x* direction in stationary regime was solved which has the same formal structure than (19) except for the time dependence. Upon rendering the recursion scheme in a pseudo-stationary form problem and thus matching the recursive structure of [22], we apply the Laplace Transform in the *t* variable, (*t* → *r*) obtaining the following pseudo-steady-state problem:

$$r\tilde{\mathbf{R}}\_0 + \overline{\boldsymbol{\pi}} \partial\_{\mathbf{x}} \tilde{\mathbf{R}}\_0 = \partial\_{\mathbf{z}} \left( \mathcal{K}\_{\mathbf{z}} \partial\_{\mathbf{z}} \tilde{\mathbf{R}}\_0 \right) - \Lambda \, \mathcal{K}\_{\mathbf{y}} \tilde{\mathbf{R}}\_0 \tag{20}$$

The *x* and *z* dependence may be separated using the same reasoning as already introduced in (8). To this end we pose the solution of problem (20) in the form:

$$
\mathbf{R}\_0 = \mathbf{P} \mathbf{C} \tag{21}
$$

where **C** = (*ζ*1(*z*), *ζ*2(*z*),...)*<sup>T</sup>* are a set of orthogonal eigenfunctions, given by *ζi*(*z*) = cos(*γlz*), and *γ<sup>i</sup>* = *iπ*/*h* (for *i* = 0, 1, 2, . . .) are the set of eigenvalues.

Replacing equation (21) in equation (20) and using the afore introduced projector (10) now for the *z* dependent degrees of freedom *<sup>h</sup>* <sup>0</sup> *dz***C**[**F**] = *<sup>h</sup>* <sup>0</sup> **<sup>F</sup>***<sup>T</sup>* <sup>∧</sup> **<sup>C</sup>** *dz* yields a first order differential equation system:

$$
\partial\_{\mathbf{x}} \mathbf{P} + \mathbf{H} \mathbf{P} = \mathbf{0} \,, \tag{22}
$$

where **P** = **P**(*x*,*r*) and **H** = **B**−<sup>1</sup> <sup>1</sup> **B**2. The entries of matrices **B**<sup>1</sup> and **B**<sup>2</sup> are

$$\begin{split} (\mathbf{B}\_{1})\_{i,j} &= -\int\_{0}^{h} \overline{u} \zeta\_{i}(z) \zeta\_{j}(z) \, dz \\ (\mathbf{B}\_{2})\_{i,j} &= \int\_{0}^{h} \partial\_{z} K\_{z} \partial\_{z} \zeta\_{i}(z) \zeta\_{j}(z) \, dz - \gamma\_{i}^{2} \int\_{0}^{h} K\_{z} \zeta\_{i}(z) \zeta\_{j}(z) \, dz \\ &\quad - r \int\_{0}^{h} \zeta\_{i}(z) \zeta\_{j}(z) \, dz - \lambda\_{i}^{2} K\_{y} \int\_{0}^{h} \zeta\_{i}(z) \zeta\_{j}(z) \, dz \, . \end{split}$$

A similar procedure leads to the source condition for (22):

6 Will-be-set-by-IN-TECH

that links the eddy diffusivity to the concentration so that we hide this dependence using a phenomenologically motivated expression for **K** which leaves us with a partial differential equation system in linear form, although the original phenomenon is non-linear. In the example below we demonstrate the closed form procedure for a problem with explicit time

The solution is generated making use of the decomposition method ([1–3]) which was originally proposed to solve non-linear partial differential equations, followed by the Laplace transform that renders the problem a pseudo-stationary one. Further we rewrite the vertical diffusivity as a time average term *K*¯ *<sup>z</sup>*(*z*) plus a term representing the variations *κz*(*z*, *t*) around the average for the time interval of the measurement *Kz*(*x*, *z*, *t*) = *K*¯ *<sup>z</sup>*(*z*) + *κz*(*z*, *t*) and use the asymptotic form of *Ky*, which is then explored to set-up the structure of the equation that

The function **R** = ∑*<sup>j</sup>* **R***<sup>j</sup>* = **1***T***R**(*c*) is now decomposed into contributions to be determined by recursion. For convenience we introduced the one-vector **1** = (1, 1, . . .)*<sup>T</sup>* and inflate the vector **R** to a vector with each element being itself a vector **R***j*. Upon inserting the expansion in equation (17) one may regroup terms that obey the recursive equations and starts with the

From the construction of the recursion equation system it is evident that other schemes are possible. The specific choice made here allows us to solve the recursion initialisation using the procedure described in reference [22], where a stationary **K** was assumed. For this reason

The boundary conditions are now used to uniquely determine the solution. In our scheme the initialisation solution that contains **R**<sup>0</sup> satisfies the boundary conditions (equations (3)) while the remaining equations satisfy homogeneous boundary conditions. Once the set of problems (19) is solved by the GILTT method, the solution of problem (1) is well determined. It is important to consider that we may control the accuracy of the results by a proper choice

In reference [22] a two dimensional problem with advection in the *x* direction in stationary regime was solved which has the same formal structure than (19) except for the time dependence. Upon rendering the recursion scheme in a pseudo-stationary form problem and thus matching the recursive structure of [22], we apply the Laplace Transform in the *t* variable,

> *Kz∂z***R**˜ <sup>0</sup>

the time dependence enters as a known source term from the first recursion step on.

*<sup>∂</sup>t***<sup>R</sup>** <sup>+</sup> *<sup>u</sup>*¯*∂x***<sup>R</sup>** <sup>−</sup> *<sup>∂</sup><sup>z</sup>* (*K*¯ *<sup>z</sup>∂z***R**) <sup>+</sup> *Ky***Λ<sup>R</sup>** <sup>=</sup> *<sup>∂</sup><sup>z</sup>* (*κz∂z***R**) (17)

*<sup>∂</sup>t***R**<sup>0</sup> <sup>+</sup> *<sup>u</sup>*¯*∂x***R**<sup>0</sup> <sup>−</sup> *<sup>∂</sup><sup>z</sup>* (*K*¯ *<sup>z</sup>∂z***R**0) <sup>+</sup> *Ky***ΛR**<sup>0</sup> <sup>=</sup> <sup>0</sup> (18)

*<sup>κ</sup>z∂z***R***j*−<sup>1</sup>

<sup>−</sup> <sup>Λ</sup> *Ky***R**˜ <sup>0</sup> (20)

. (19)

+ *Ky***ΛR***<sup>j</sup>* = *∂<sup>z</sup>*

dependence, which is novel in the literature.

defines the recursive decomposition scheme:

The extension to the closed form recursion is then given by

 *K*¯ *<sup>z</sup>∂z***R***<sup>j</sup>* 

*∂t***R***<sup>j</sup>* + *u*¯*∂x***R***<sup>j</sup>* − *∂<sup>z</sup>*

time averaged solution for *Kz*:

**3.3. Recursion initialisation**

of the number of terms in the solution series.

(*t* → *r*) obtaining the following pseudo-steady-state problem:

*r***R**˜ <sup>0</sup> + *u∂x***R**˜ <sup>0</sup> = *∂<sup>z</sup>*

$$\mathbf{P}(0,r) = Q\mathbf{B}\_1^{-1} \int dz \mathbf{C}[\delta(z - H\_\mathcal{S})] \int dy \mathbf{Y}[\delta(y - y\_0)] = Q\mathbf{B}\_1^{-1} \left(\mathbf{C}(H\_\mathcal{S}) \wedge \mathbf{1}\right) \left(\mathbf{1} \wedge \mathbf{Y}(y\_0)\right) \tag{23}$$

Following the reasoning of [22] we solve (22) applying Laplace transform and diagonalisation of the matrix **H** = **XDX**−<sup>1</sup> which results in

$$\mathbf{P}(s,r) = \mathbf{X}(s\mathbf{I}+\mathbf{D})^{-1}\mathbf{X}^{-1}\mathbf{P}(0,r) \tag{24}$$

where **P**˜(*s*,*r*) denotes the Laplace Transform of **P**(*x*,*r*). Here **X**(−1) is the (inverse) matrix of the eigenvectors of matrix **B**−<sup>1</sup> <sup>1</sup> **B**<sup>2</sup> with diagonal eigenvalue matrix **D** and the entries of matrix (*s***I** + **D**)*ii* = *s* + *di*. After performing the Laplace transform inversion of equation (24), we come out with

$$\mathbf{P}(\mathbf{x},r) = \mathbf{X}\mathbf{G}(\mathbf{x},r)\mathbf{X}^{-1}\mathbf{D}\ ,\tag{25}$$

where **G**(*x*,*r*) is the diagonal matrix with components (**G**)*ii* = *e*−*dix*. Further the still unknown arbitrary constant matrix is given by **Ω** = **X**−1**P**(0,*r*).

The analytical time dependence for the recursion initialisation (20) is obtained upon applying the inverse Laplace transform definition

$$\mathbf{R}\_0(\mathbf{x}, z, t) = \frac{1}{2\pi i} \int\_{\gamma - i\infty}^{\gamma + i\infty} \mathbf{P}(\mathbf{x}, r) \mathbf{C}(z) e^{rt} \, dr \,. \tag{26}$$

To overcome the drawback of evaluating the line integral appearing in the above solution, we perform the calculation of this integral by the Gaussian quadrature scheme, which is exact if the integrand is a polynomial of degree 2*<sup>M</sup>* <sup>−</sup> 1 in the <sup>1</sup> *<sup>r</sup>* variable

$$\mathbf{R}\_0(\mathbf{x}, z, t) = \frac{1}{t} \mathbf{a}^T \left( \mathbf{p} \mathbf{R}\_0(\mathbf{x}, z, \frac{\mathbf{p}}{t}) \right) \; , \tag{27}$$

where **a** and **p** are respectively vectors with the weights and roots of the Gaussian quadrature scheme ([27]), and the argument (*x*, *z*, **<sup>p</sup>** *<sup>t</sup>* ) signifies the *k*-th component of **p** in the *k*-th row of **pR**0. Note, *k* is a component from contraction with **a**.

*Ky* =

1 3 *�* = <sup>1</sup> <sup>−</sup> *<sup>z</sup> h* <sup>2</sup> − *z L* <sup>−</sup> <sup>2</sup> 3

*qv* <sup>=</sup> 4.16 *<sup>z</sup>*

transverse wave peak.

*<sup>h</sup>* , *<sup>ψ</sup>*

related to the intensity of turbulence ([16]).

<sup>√</sup>*πσ<sup>v</sup>* 16(*fm*)*vqv*

The wind speed profile can be described by a power law *uz*/*u*<sup>1</sup> = (*z*/*z*1)

source point in an unstable/neutral atmospheric boundary layer.

this simulation did not make use of the terrain's realistic complexity.

0,0 0,2 0,4 0,6 0,8 1,0

Experiment 2

Cp (Bq/m <sup>3</sup> )

approach for the experiment; dotted lines indicate a factor of two.

0,0

0,2

0,4

Co (Bq/m 3

)

0,6

0,8

1,0

with *σ*<sup>2</sup>

where *σv* is the standard deviation of the longitudinal turbulent velocity component, *qv* is the stability function, *ψ�* is the dimensionless molecular dissipation rate and (*fm*)*<sup>v</sup>* is the

and *u*<sup>1</sup> are the horizontal mean wind speeds at heights *z* and *z*<sup>1</sup> and *n* is an exponent that is

Thus, in this study we introduce the vertical and lateral eddy diffusivities (eq. (35) and eq. (29)) and the power law wind profile in the 3D-GILTT model (eq. (16) or equivalently eq. (20)) to calculate the ground-level concentration of emissions released from an elevated continuous

The validation of the 3D-GILTT model predictions against experimental data from the Angra site together with a two dimensional model (GILTTG) are shown in table 2. While the present approach (3D-GILTT) is based on a genuine three dimensional description an earlier analytical approach (GILTTG) uses a Gaussian assumption for the horizontal transverse direction ([22]). Figure 1 shows the comparison of predicted concentrations against observed ones for the three dimensional approach, which reproduces acceptably the observed concentrations, although

In the further we use the standard statistical indices in order to compare the quality of the two approaches. Note, that we present the two analytical model approaches, since the earlier

**Figure 1.** Observed and predicted scatter diagram of ground-level concentrations using the 3D-GILTT

Co (Bq/m 3

)

<sup>+</sup> 0.75 <sup>1</sup> 2

*<sup>v</sup>* <sup>=</sup> 0.98*cv* (*fm*) 2 3 *v*

 *ψ� qv*

On an Analytical Model for the Radioactive Contaminant Release in the Atmosphere from Nuclear Power Plants

> 2 <sup>3</sup> *z h* 2 3 *w*2 ∗

0 10 20 30 40 50 60 70 80

Experiment 3

Cp (Bq/m <sup>3</sup> )

and (*fm*)*<sup>v</sup>* = 0.16 (29)

*<sup>n</sup>* ([25]), where *uz*

227
