**4. Experimental data and turbulent parameterisation**

For model validation we chose a controlled release of radioactive material performed in 1985 at the Itaorna Beach, close to the nuclear reactor site Angra dos Reis in the Rio de Janeiro state, Brazil. Details of the dispersion experiment is described elsewhere ([5]). The experiment consisted in the controlled releases of radioactive tritiated water vapour from the meteorological tower at 100*m* height during five days (28 November to 4 December 1984). During the whole experiment, four meteorological towers collected the relevant meteorological data. Wind speed and direction were measured at three levels (10*m*, 60*m* and 100*m*) together with the temperature gradients between 10*m* and 100*m*. Some additional data of relative humidity were available in some of the sampling sites, and were used to calculate the concentration of radioactive tritiated water in the air (after measuring the radioactivity of the collected samples). All relevant details, as well as the synoptic meteorological conditions during the dispersion campaign are described in ref. [5]. The data from experiments 2 and 3 were used to obtain the numerical results and are presented in table 1.


**Table 1.** Micro-meteorological parameters and emission rate for experiments 2 and 3 at third period.

The micro-meteorological parameters shown in table 1 are calculated from equations obtained in the literature. The roughness length utilized was 1*m* and the Monin-Obukhov length for convective conditions can be written as *L* = −*h*/*k* (*u*∗/*w*∗) <sup>3</sup> ([35]), where *k* is the von Karman constant (*k* = 0.4), *w*<sup>∗</sup> is the convective velocity scale with wind speed *U*, *u*<sup>∗</sup> = *kU*/ln(*zr*/*z*0) is the friction velocity, where *U* is the wind velocity at the reference height *zr* = 10*m*, and *<sup>h</sup>* <sup>=</sup> 0.3*u*∗/ *fc* is the height of the boundary layer with the Coriolis coefficient *fc* <sup>=</sup> <sup>10</sup>−4.

In the atmospheric diffusion problems the choice of a turbulent parameterisation represents a fundamental aspect for contaminant dispersion modelling. From the physical point of view a turbulence parameterisation an approximation for the natural phenomenon, where details are hidden in the parameters used, that have to be adjusted in order to reproduce experimental findings. The reliability of each model strongly depends on the way the turbulent parameters are calculated and related to the current understanding of the planetary boundary layer. In terms of the convective scaling parameters the vertical and lateral eddy diffusivities can be formulated as follows ([11]):

$$K\_2 = 0.22 w\_\* h \left(\frac{z}{h}\right)^{\frac{1}{3}} \left(1 - \frac{z}{h}\right)^{\frac{1}{3}} \left(1 - e^{\frac{4\tilde{\omega}}{h}} - 0.0003 e^{\frac{8\tilde{\omega}}{h}}\right) \tag{28}$$

226 Nuclear Power – Practical Aspects On an Analytical Model for the Radioactive Contaminant Release in the Atmosphere from Nuclear Power Plants <sup>9</sup> 227 On an Analytical Model for the Radioactive Contaminant Release in the Atmosphere from Nuclear Power Plants

$$K\_{\mathcal{Y}} = \frac{\sqrt{\pi}\sigma\_{\upsilon}}{16(f\_{\mathcal{W}})\_{\upsilon}q\_{\upsilon}} \qquad \text{with} \qquad \sigma\_{\upsilon}^{2} = \frac{0.98c\_{\upsilon}}{(f\_{\mathcal{W}})\_{\upsilon}^{\frac{2}{3}}} \left(\frac{\psi\_{\varepsilon}}{q\_{\upsilon}}\right)^{\frac{2}{3}} \left(\frac{z}{h}\right)^{\frac{2}{3}} w\_{\ast}^{2}$$

$$q\_{\upsilon} = 4.16 \frac{z}{h}, \qquad \psi\_{\varepsilon}^{\frac{1}{3}} = \left(\left(1 - \frac{z}{h}\right)^{2} \left(-\frac{z}{L}\right)^{-\frac{2}{3}} + 0.75\right)^{\frac{1}{2}} \qquad \text{and} \qquad (f\_{\mathcal{W}})\_{\upsilon} = 0.16 \tag{29}$$

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 transverse wave peak.

8 Will-be-set-by-IN-TECH

where **a** and **p** are respectively vectors with the weights and roots of the Gaussian quadrature

For model validation we chose a controlled release of radioactive material performed in 1985 at the Itaorna Beach, close to the nuclear reactor site Angra dos Reis in the Rio de Janeiro state, Brazil. Details of the dispersion experiment is described elsewhere ([5]). The experiment consisted in the controlled releases of radioactive tritiated water vapour from the meteorological tower at 100*m* height during five days (28 November to 4 December 1984). During the whole experiment, four meteorological towers collected the relevant meteorological data. Wind speed and direction were measured at three levels (10*m*, 60*m* and 100*m*) together with the temperature gradients between 10*m* and 100*m*. Some additional data of relative humidity were available in some of the sampling sites, and were used to calculate the concentration of radioactive tritiated water in the air (after measuring the radioactivity of the collected samples). All relevant details, as well as the synoptic meteorological conditions during the dispersion campaign are described in ref. [5]. The data from experiments 2 and 3

> Exp Period *U*(*m*/*s*) *h*(*m*) *u*∗(*m*/*s*) *L*(*m*) *w*∗(*m*/*s*) *Q*(*MBq*/*s*) 2 3 2.2 1134 0.4 −951 0.6 25.3 3 3 2.6 1367 0.5 −1147 0.7 20.5

**Table 1.** Micro-meteorological parameters and emission rate for experiments 2 and 3 at third period.

The micro-meteorological parameters shown in table 1 are calculated from equations obtained in the literature. The roughness length utilized was 1*m* and the Monin-Obukhov length for

constant (*k* = 0.4), *w*<sup>∗</sup> is the convective velocity scale with wind speed *U*, *u*<sup>∗</sup> = *kU*/ln(*zr*/*z*0) is the friction velocity, where *U* is the wind velocity at the reference height *zr* = 10*m*, and *<sup>h</sup>* <sup>=</sup> 0.3*u*∗/ *fc* is the height of the boundary layer with the Coriolis coefficient *fc* <sup>=</sup> <sup>10</sup>−4.

In the atmospheric diffusion problems the choice of a turbulent parameterisation represents a fundamental aspect for contaminant dispersion modelling. From the physical point of view a turbulence parameterisation an approximation for the natural phenomenon, where details are hidden in the parameters used, that have to be adjusted in order to reproduce experimental findings. The reliability of each model strongly depends on the way the turbulent parameters are calculated and related to the current understanding of the planetary boundary layer. In terms of the convective scaling parameters the vertical and lateral eddy diffusivities can be

*<sup>t</sup>* ) signifies the *k*-th component of **p** in the *k*-th row of

<sup>3</sup> ([35]), where *k* is the von Karman

*<sup>h</sup>* − 0.0003*e*

8*z h* 

(28)

scheme ([27]), and the argument (*x*, *z*, **<sup>p</sup>**

**pR**0. Note, *k* is a component from contraction with **a**.

**4. Experimental data and turbulent parameterisation**

were used to obtain the numerical results and are presented in table 1.

convective conditions can be written as *L* = −*h*/*k* (*u*∗/*w*∗)

*Kz* = 0.22*w*∗*h*

 *z h* 1 <sup>3</sup> <sup>1</sup> <sup>−</sup> *<sup>z</sup> h* 1 <sup>3</sup> 1 − *e* 4*z*

formulated as follows ([11]):

The wind speed profile can be described by a power law *uz*/*u*<sup>1</sup> = (*z*/*z*1) *<sup>n</sup>* ([25]), where *uz* 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 related to the intensity of turbulence ([16]).

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 source point in an unstable/neutral atmospheric boundary layer.

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 this simulation did not make use of the terrain's realistic complexity.

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 approach for the experiment; dotted lines indicate a factor of two.


Statistical Indices GILTTG 3D-GILTT

NMSE = (*Co* − *Cp*)2/*Cp Co* 1.34 0.38 COR = (*Co* − *Co*)(*Cp* − *Cp*)/*σoσ<sup>p</sup>* 0.67 0.83 FA2 = 0.5 ≤ (*Cp*/*Co*) ≤ 2 0.53 0.88 FA5 = 0.2 ≤ (*Cp*/*Co*) ≤ 5 0.96 1.00 FB = *Co* − *Cp*/0.5(*Co* + *Cp*) −0.44 0.13 FS = (*σ<sup>o</sup>* − *σp*)/0.5(*σ<sup>o</sup>* + *σp*) −0.54 0.18

with *R*<sup>2</sup> = 0.67 and *κ* = 0.43, whereas the 3D-GILTT obeys the result *C*¯ *<sup>p</sup>* = 0.69*C*¯*<sup>o</sup>* + 3.26 with *R*<sup>2</sup> = 0.83 and *κ* = 0.36. In order to perform a model validation we introduced an index

match between the model and the experimental findings. Here *a* is the slope, *b* the intersection, *Coi* of the experimental data and *C*¯*<sup>o</sup>* its arithmetic mean. Since the experiment is of stochastic character whereas the stochastic properties are hidden in the model parameters, considerable fluctuations are present. Nevertheless, by comparison one observes that the present approach

0 10 20 30 40 50 60

Co

The consistency of the K-approach strongly depends on the way the eddy diffusivity is determined on the basis of the turbulence structure of the PBL and on the model ability

**Figure 2.** Linear regression for the GILTTG and 3D-GILTT. The bisector was added as an eye guide.

*<sup>i</sup>*=<sup>1</sup> *Coi*, which if identical zero indicates a perfect

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

229

 GILTTG 3D-GILTT

*<sup>n</sup>* <sup>∑</sup>*<sup>n</sup>*

**Table 3.** Statistical comparisons between GILTTG and 3D-GILTT results.

<sup>2</sup> with *<sup>C</sup>*¯*<sup>o</sup>* = <sup>1</sup>

*κ* = 

(*<sup>a</sup>* <sup>−</sup> <sup>1</sup>)<sup>2</sup> <sup>+</sup> (*b*/*C*¯*o*)

yields the better description of the data.

0

**5. Meso-scale simulation for** *K***-closure**

10

20

30

Cp

40

50

60

**Table 2.** Concentrations of nine runs with various positions of the Angra dos Reis experiment and model prediction by the approaches GILTTG and 3D-GILTT.

one was found to be acceptable in comparison to other approaches found in the literature and both give a solution in closed form. The standard statistical indices are NMSE, the normalized mean square error; COR, the correlation coefficient; FA2 and FA5, the fraction of data (in %) in the cones determined by a factor of two and five, respectively; FB, the fractional bias and FS, the fractional standard deviation. The subscripts *o* and *p* refer to observed and predicted quantities, respectively, and *C*¯ indicates the averaged values. Table 3 presents the results of the statistical indices used to evaluate the model performance ([14]) and further compare our model to the GILTTG approach. The statistical index FB indicates weather the predicted quantities (*Cp*) under- or overestimates the observed ones (*Co*). The statistical index NMSE represents the quadratic error of the predicted quantities in relation to the observed ones. Best results are indicated by values compatible with zero for NMSE, FB and FS, and compatible with unity for COR, FA2 and FA5. The statistical indices point out that a reasonable agreement is obtained between experimental data and the 3D-GILTT model.

In order to validate the two models we fit the predicted versus observed values by a linear regression (see figure 2), where the closer their intersect to the origin and the closer the slope is to unity the better is the approach. The GILTTG approach results in *C*¯ *<sup>p</sup>* = 1.16*C*¯*<sup>o</sup>* + 7.01


**Table 3.** Statistical comparisons between GILTTG and 3D-GILTT results.

10 Will-be-set-by-IN-TECH

Exp. Period Distance (*m*) Observed (*Bq*/*m*3) Predictions (*Bq*/*m*3)

 3 610 0.58 0.20 0.40 3 600 0.50 0.19 0.40 3 700 0.53 0.29 0.44 3 815 0.61 0.38 0.47 3 970 0.54 0.47 0.48 3 1070 0.86 0.51 0.48 3 750 0.39 0.33 0.46 3 935 0.40 0.45 0.48 3 705 38.89 47.18 31.13 3 700 24.09 46.53 31.02 3 815 48.95 59.98 32.66 3 970 36.22 73.03 32.95 3 1070 33.50 78.65 32.44 3 500 50.26 17.74 22.58 3 375 26.86 2.57 11.67 3 960 19.61 72.35 32.97 3 915 18.02 69.04 33.03

**Table 2.** Concentrations of nine runs with various positions of the Angra dos Reis experiment and

one was found to be acceptable in comparison to other approaches found in the literature and both give a solution in closed form. The standard statistical indices are NMSE, the normalized mean square error; COR, the correlation coefficient; FA2 and FA5, the fraction of data (in %) in the cones determined by a factor of two and five, respectively; FB, the fractional bias and FS, the fractional standard deviation. The subscripts *o* and *p* refer to observed and predicted quantities, respectively, and *C*¯ indicates the averaged values. Table 3 presents the results of the statistical indices used to evaluate the model performance ([14]) and further compare our model to the GILTTG approach. The statistical index FB indicates weather the predicted quantities (*Cp*) under- or overestimates the observed ones (*Co*). The statistical index NMSE represents the quadratic error of the predicted quantities in relation to the observed ones. Best results are indicated by values compatible with zero for NMSE, FB and FS, and compatible with unity for COR, FA2 and FA5. The statistical indices point out that a reasonable agreement

In order to validate the two models we fit the predicted versus observed values by a linear regression (see figure 2), where the closer their intersect to the origin and the closer the slope is to unity the better is the approach. The GILTTG approach results in *C*¯ *<sup>p</sup>* = 1.16*C*¯*<sup>o</sup>* + 7.01

model prediction by the approaches GILTTG and 3D-GILTT.

is obtained between experimental data and the 3D-GILTT model.

GILTTG 3D-GILTT

with *R*<sup>2</sup> = 0.67 and *κ* = 0.43, whereas the 3D-GILTT obeys the result *C*¯ *<sup>p</sup>* = 0.69*C*¯*<sup>o</sup>* + 3.26 with *R*<sup>2</sup> = 0.83 and *κ* = 0.36. In order to perform a model validation we introduced an index *κ* = (*<sup>a</sup>* <sup>−</sup> <sup>1</sup>)<sup>2</sup> <sup>+</sup> (*b*/*C*¯*o*) <sup>2</sup> with *<sup>C</sup>*¯*<sup>o</sup>* = <sup>1</sup> *<sup>n</sup>* <sup>∑</sup>*<sup>n</sup> <sup>i</sup>*=<sup>1</sup> *Coi*, which if identical zero indicates a perfect match between the model and the experimental findings. Here *a* is the slope, *b* the intersection, *Coi* of the experimental data and *C*¯*<sup>o</sup>* its arithmetic mean. Since the experiment is of stochastic character whereas the stochastic properties are hidden in the model parameters, considerable fluctuations are present. Nevertheless, by comparison one observes that the present approach yields the better description of the data.

**Figure 2.** Linear regression for the GILTTG and 3D-GILTT. The bisector was added as an eye guide.

#### **5. Meso-scale simulation for** *K***-closure**

The consistency of the K-approach strongly depends on the way the eddy diffusivity is determined on the basis of the turbulence structure of the PBL and on the model ability

#### 12 Will-be-set-by-IN-TECH 230 Nuclear Power – Practical Aspects On an Analytical Model for the Radioactive Contaminant Release in the Atmosphere from Nuclear Power Plants <sup>13</sup>

to reproduce experimental diffusion data. Keeping the K-theory limitations in mind many efforts have been made to develop turbulent parametrisations for practical applications in air pollution modelling which reveals the essential features of turbulent diffusion, but which as far as possible preserves the simplicity and flexibility of the K-theory formulation. The aim of this step is to elaborate parametrisations for the eddy diffusivity coefficients in the PBL based on the micro-meteorological parameters that were extracted from mesoscale WRF simulations. The WRF model is based on the Taylor's statistical theory and a model for Eulerian spectra ([11, 24]). The main idea of the proposed spectral model relies on considering the turbulent spectra as a superposition of a buoyant produced part (with a convective peak wavelength) and a shear produced part (with a mechanical peak wavelength). By such a model, the plume spreading rate is directly connected with the spectral distribution of eddies in the PBL, that is with the energy containing eddies of the turbulence.

where *F<sup>E</sup>*

[11]

where <sup>Ψ</sup>*�<sup>b</sup>* <sup>=</sup> *�bh*

variances *σ*<sup>2</sup>

plant.

dimensional spectrum as *S<sup>E</sup>*

buoyancy and shear, respectively.

while for mechanical turbulence ([10])

*w*<sup>3</sup> ∗

wind speed *u*¯ in the mixing layer. The dimensionless spectrum *F<sup>E</sup>*

spectra with the total variance, *σ*<sup>2</sup>

*<sup>i</sup>* <sup>=</sup> <sup>∞</sup>

<sup>0</sup> (*S<sup>E</sup>*

*Kz* <sup>=</sup> *<sup>β</sup><sup>i</sup>* 4

and eqn. (34) one ends up with *<sup>K</sup><sup>α</sup>* = *<sup>β</sup><sup>i</sup>*

and Φ*�<sup>s</sup>* = *�κ<sup>z</sup>*

*<sup>i</sup>* (0) is the value of the normalised Eulerian energy spectrum at *n* = 0. In this

*<sup>i</sup>* (*n*) = *Sib*(*n*) + *Sis*(*n*), where the subscripts *b* and *s* stand for

 2 3 *w*2

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

are the dimensional dissipation rate functions, *�<sup>b</sup>* and *�<sup>s</sup>*

*<sup>i</sup>* (*n*) in eqn. (31) is obtained by normalizing the dimensional

*is*(*n*)

1.5*cw z u*¯ 

*n*(*fmw*) 5 3 Φ2 3 *�su*<sup>2</sup> ∗ 

<sup>∗</sup> , (32)

231

<sup>∗</sup> (33)

*mi* is the normalized frequency of

. (34)

(35)

*is*. Making use of eqns. (30), (32), (33)

, that for the *w*-component becomes

2 3 *�b z zi*

way the eddy diffusivity is directly associated to the energy-containing eddies which are the principal contribution to turbulent transport. In order to use eqn. (31) we have to find an analytical form for the dimensionless Eulerian spectrum. We assume here that the spectral distribution of turbulent kinetic energy is a superposition of buoyancy and shear components. Such a TKE model may be evaluated as a good approximation for a real PBL, where turbulent production is due to both mechanisms ([15, 20]). In these conditions we may write the Eulerian

An analytical form for the dimensional spectra in convective turbulence has been reported in

 *nz u*¯ 

<sup>1</sup> <sup>+</sup> 1.5 *nz u*¯ *f* ∗ *mi* Ψ

the spectral peaks regardless of stratification and *fmi* is the reduced frequency with the mean

<sup>=</sup> *<sup>S</sup><sup>E</sup>*

The total wind velocity variance is obtained by the sum of mechanical and convective

*ib*(0) + *<sup>S</sup><sup>E</sup>*

In order to illustrate the suitability of the discussed formulation to simulate contaminant dispersion in the atmospheric boundary layer, we evaluate the performance of the new solution and simulate radioactive substance dispersion around the Fukushima-Daiichi power

*ib* <sup>+</sup> *<sup>σ</sup>*<sup>2</sup>

 2 3 *w*2 <sup>∗</sup> +

*<sup>i</sup>* (*n*)*dn*, that is

*ib*(*n*) + *<sup>S</sup><sup>E</sup>*

*σ*2 *i*

*is*(0) 

 *nz u*¯ 

<sup>1</sup> <sup>+</sup> 1.5 *nz u*¯ *fmi* Φ2 3 *�su*<sup>2</sup>

*Sib*(*n*) = 0.98*ci*

*n f* ∗ *mi* 5 3 

*u*3 ∗

*FE*

*n*(*f* ∗ *mw*) 5 3 Ψ 2 3 *�b z zi*

**6. Application to the Fukushima-Daiichi accident**

*ib*(*n*) + *<sup>S</sup><sup>E</sup>*

 0.98*cw*

*<sup>i</sup>* <sup>=</sup> <sup>∞</sup> <sup>0</sup> *<sup>S</sup><sup>E</sup>*

*<sup>i</sup>* (*n*) = *<sup>S</sup><sup>E</sup> i σ*2 *i*

*is*(*n*)) *dn* <sup>=</sup> *<sup>σ</sup>*<sup>2</sup>

4 *SE*

 *z u*¯ 

are the convective and mechanical rate of *tke* dissipation, *f* ∗

*Sis*(*n*) = 1.5*ci*

*n* (*fmi*) 5 3 

The WRF Simulator is a meso-scale numerical weather prediction system that features multiple dynamical cores and a 3-dimensional variational data assimilation system. The simulator offers multiple physics options that can be combined in various ways. Since this study focusses on the implementation of an interface with a model for the PBL, orography related features of WRF were of importance, more specifically the Land-Surface and PBL physics options were chosen for the present study. In WRF, when a PBL scheme is activated, a specific vertical diffusion is de-activated with the assumption that the PBL scheme will handle this process. The Mellor-Yamada-Janjic PBL scheme derives the eddy diffusivities coefficients and the boundary layer height from the estimations of the Turbulent Kinetic Energy (TKE) through the full range of atmospheric turbulent regimes ([19]).

Two grids were used for the WRF meso-scale simulation. The outer grid has an extension of the order of half the earth radius so that a significant part of the large scale geological domain of interest is included. The inner grid is centred at the point of interest, i.e. the centre of the power plant where typically the nuclear reactor is located. The simulation may in principal contain a sequence of days or even months. The micro-meteorological data are extracted at the centre point of the inner WRF grid. The spectral model needs these quantities to calculate the eddy diffusivity coefficients.

On the basis of Taylor's theory, Taylor proposed that under the hypothesis of homogeneous turbulence, the eddy diffusivities may be expressed as

$$K\_{\underline{a}} = \frac{d}{dt} \left( \frac{\sigma\_{\underline{a}}^2}{2} \right) = \frac{\sigma\_{\underline{i}}^2 \beta\_{\underline{i}}}{2\pi} \int\_0^\infty F\_{\underline{i}}^E(n) \frac{\sin(2\pi nt \beta\_{\underline{i}}^{-1})}{n} dn,\tag{30}$$

where *α* = (*x*, *y*, *z*) and *i* = *u*, *v*, *w*, *F<sup>E</sup> <sup>i</sup>* (*n*) is the value of the Eulerian spectrum of energy normalized by the Eulerian velocity variance, and *σ*<sup>2</sup> *<sup>i</sup>* corresponds to the Eulerian variance of the turbulent wind field. Following [33], *β<sup>i</sup>* = *<sup>π</sup>U*<sup>2</sup> 16*σ*<sup>2</sup> *i* 1 2 . For large diffusion travel times (*t* → ∞), the filter function in the integral of eqn. (30) selects *F<sup>E</sup> <sup>i</sup>* (*n*) at the origin of the frequency space, such that the rate of dispersion becomes independent of the travel time from the source and can be expressed as a function of local properties of turbulence,

$$K\_{\mathfrak{A}} = \frac{\sigma\_i^2 \beta\_i F\_i^E(0)}{4} \tag{31}$$

where *F<sup>E</sup> <sup>i</sup>* (0) is the value of the normalised Eulerian energy spectrum at *n* = 0. In this way the eddy diffusivity is directly associated to the energy-containing eddies which are the principal contribution to turbulent transport. In order to use eqn. (31) we have to find an analytical form for the dimensionless Eulerian spectrum. We assume here that the spectral distribution of turbulent kinetic energy is a superposition of buoyancy and shear components. Such a TKE model may be evaluated as a good approximation for a real PBL, where turbulent production is due to both mechanisms ([15, 20]). In these conditions we may write the Eulerian dimensional spectrum as *S<sup>E</sup> <sup>i</sup>* (*n*) = *Sib*(*n*) + *Sis*(*n*), where the subscripts *b* and *s* stand for buoyancy and shear, respectively.

An analytical form for the dimensional spectra in convective turbulence has been reported in [11]

$$S\_{ib}(n) = \frac{0.98c\_i\left(\frac{nz}{a}\right)}{n\left(f\_{mi}^\*\right)^{\frac{5}{7}}\left(1 + 1.5\frac{\frac{nz}{a}}{f\_{mi}^\*}\right)}\Psi\_{\epsilon b}^{\frac{2}{3}}\left(\frac{z}{z\_i}\right)^{\frac{2}{3}}w\_\*^2\,\,\,\tag{32}$$

while for mechanical turbulence ([10])

12 Will-be-set-by-IN-TECH

to reproduce experimental diffusion data. Keeping the K-theory limitations in mind many efforts have been made to develop turbulent parametrisations for practical applications in air pollution modelling which reveals the essential features of turbulent diffusion, but which as far as possible preserves the simplicity and flexibility of the K-theory formulation. The aim of this step is to elaborate parametrisations for the eddy diffusivity coefficients in the PBL based on the micro-meteorological parameters that were extracted from mesoscale WRF simulations. The WRF model is based on the Taylor's statistical theory and a model for Eulerian spectra ([11, 24]). The main idea of the proposed spectral model relies on considering the turbulent spectra as a superposition of a buoyant produced part (with a convective peak wavelength) and a shear produced part (with a mechanical peak wavelength). By such a model, the plume spreading rate is directly connected with the spectral distribution of eddies in the PBL, that is

The WRF Simulator is a meso-scale numerical weather prediction system that features multiple dynamical cores and a 3-dimensional variational data assimilation system. The simulator offers multiple physics options that can be combined in various ways. Since this study focusses on the implementation of an interface with a model for the PBL, orography related features of WRF were of importance, more specifically the Land-Surface and PBL physics options were chosen for the present study. In WRF, when a PBL scheme is activated, a specific vertical diffusion is de-activated with the assumption that the PBL scheme will handle this process. The Mellor-Yamada-Janjic PBL scheme derives the eddy diffusivities coefficients and the boundary layer height from the estimations of the Turbulent Kinetic Energy (TKE)

Two grids were used for the WRF meso-scale simulation. The outer grid has an extension of the order of half the earth radius so that a significant part of the large scale geological domain of interest is included. The inner grid is centred at the point of interest, i.e. the centre of the power plant where typically the nuclear reactor is located. The simulation may in principal contain a sequence of days or even months. The micro-meteorological data are extracted at the centre point of the inner WRF grid. The spectral model needs these quantities to calculate

On the basis of Taylor's theory, Taylor proposed that under the hypothesis of homogeneous

 ∞ 0 *FE <sup>i</sup>* (*n*)

> *<sup>π</sup>U*<sup>2</sup> 16*σ*<sup>2</sup> *i* 1 2

space, such that the rate of dispersion becomes independent of the travel time from the source

*<sup>i</sup> <sup>β</sup>iF<sup>E</sup> <sup>i</sup>* (0)

*<sup>K</sup><sup>α</sup>* <sup>=</sup> *<sup>σ</sup>*<sup>2</sup>

sin(2*πntβ*−<sup>1</sup>

*n*

*<sup>i</sup>* )

*<sup>i</sup>* (*n*) is the value of the Eulerian spectrum of energy

*<sup>i</sup>* corresponds to the Eulerian variance of

<sup>4</sup> (31)

. For large diffusion travel times (*t* →

*<sup>i</sup>* (*n*) at the origin of the frequency

*dn*, (30)

with the energy containing eddies of the turbulence.

through the full range of atmospheric turbulent regimes ([19]).

turbulence, the eddy diffusivities may be expressed as

 *σ*<sup>2</sup> *α* 2 <sup>=</sup> *<sup>σ</sup>*<sup>2</sup> *<sup>i</sup> β<sup>i</sup>* 2*π*

*<sup>K</sup><sup>α</sup>* <sup>=</sup> *<sup>d</sup> dt*

normalized by the Eulerian velocity variance, and *σ*<sup>2</sup>

∞), the filter function in the integral of eqn. (30) selects *F<sup>E</sup>*

and can be expressed as a function of local properties of turbulence,

the turbulent wind field. Following [33], *β<sup>i</sup>* =

where *α* = (*x*, *y*, *z*) and *i* = *u*, *v*, *w*, *F<sup>E</sup>*

the eddy diffusivity coefficients.

$$S\_{\rm is}(n) = \frac{1.5c\_i \left(\frac{nz}{\hbar}\right)}{n \left(f\_{mi}\right)^{\frac{\pi}{\hbar}} \left(1 + 1.5 \frac{\frac{n}{\pi}}{f\_{mi}}\right)} \Phi\_{\rm es}^{\frac{2}{3}} u\_\*^2 \tag{33}$$

where <sup>Ψ</sup>*�<sup>b</sup>* <sup>=</sup> *�bh w*<sup>3</sup> ∗ and Φ*�<sup>s</sup>* = *�κ<sup>z</sup> u*3 ∗ are the dimensional dissipation rate functions, *�<sup>b</sup>* and *�<sup>s</sup>* are the convective and mechanical rate of *tke* dissipation, *f* ∗ *mi* is the normalized frequency of the spectral peaks regardless of stratification and *fmi* is the reduced frequency with the mean wind speed *u*¯ in the mixing layer.

The dimensionless spectrum *F<sup>E</sup> <sup>i</sup>* (*n*) in eqn. (31) is obtained by normalizing the dimensional spectra with the total variance, *σ*<sup>2</sup> *<sup>i</sup>* <sup>=</sup> <sup>∞</sup> <sup>0</sup> *<sup>S</sup><sup>E</sup> <sup>i</sup>* (*n*)*dn*, that is

$$F\_i^E(n) = \frac{S\_i^E}{\sigma\_i^2} = \frac{S\_{ib}^E(n) + S\_{is}^E(n)}{\sigma\_i^2} \,. \tag{34}$$

The total wind velocity variance is obtained by the sum of mechanical and convective variances *σ*<sup>2</sup> *<sup>i</sup>* <sup>=</sup> <sup>∞</sup> <sup>0</sup> (*S<sup>E</sup> ib*(*n*) + *<sup>S</sup><sup>E</sup> is*(*n*)) *dn* <sup>=</sup> *<sup>σ</sup>*<sup>2</sup> *ib* <sup>+</sup> *<sup>σ</sup>*<sup>2</sup> *is*. Making use of eqns. (30), (32), (33) and eqn. (34) one ends up with *<sup>K</sup><sup>α</sup>* = *<sup>β</sup><sup>i</sup>* 4 *SE ib*(0) + *<sup>S</sup><sup>E</sup> is*(0) , that for the *w*-component becomes

$$K\_{z} = \frac{\beta\_{i}}{4} \left( \frac{0.98 c\_{w} \left( \frac{z}{\hbar} \right)}{n \left( f\_{mw}^{\*} \right)^{\frac{5}{3}}} \Psi\_{\epsilon b}^{\frac{2}{3}} \left( \frac{z}{z\_{i}} \right)^{\frac{2}{3}} w\_{\*}^{2} + \frac{1.5 c\_{w} \left( \frac{z}{\hbar} \right)}{n \left( f\_{mw} \right)^{\frac{5}{3}}} \Phi\_{\epsilon s}^{\frac{4}{3}} u\_{\*}^{2} \right) \tag{35}$$

#### **6. Application to the Fukushima-Daiichi accident**

In order to illustrate the suitability of the discussed formulation to simulate contaminant dispersion in the atmospheric boundary layer, we evaluate the performance of the new solution and simulate radioactive substance dispersion around the Fukushima-Daiichi power plant.

At the 11th of march, 2011 the Fukushima-Daiichi nuclear power plant accident (coordinates in latitude, longitude: 37*<sup>o</sup>* 25' 17" N, 141*<sup>o</sup>* 1' 57" E) caused considerable radiation leakage into the atmosphere and into the sea. The radioactive pollution of the environment and sea was caused principally by the direct release of contaminated water from the power station. To a lesser extent atmospheric release of the radio-nuclide from the atmospheric plume are carried by the winds over the sea during and after the accident sequence. Shorter–lived radioactive elements, such as Iodine-131 were detectable for a few months (half-live of approximately 80 days). Others, such as Ruthenium-106 and Caesium-134 will still persist in the environment for several years (Caesium-137 has a half-life of approximately 30 years).

constant release of radioactive material, the second and third plot correspond to 48 hours and

From the meso-scale meteorological data one may determine the eddy diffusivity coefficients for each specific hour. In the *z*-time plot of fig. 4 we report the dimensionless vertical eddy diffusivity coefficients as calculated by eq. (35) for four subsequent days. The figure shows in a simple way the spatial and time structure of this coefficient. In this context, it is important to point out that the largest values of the *Kz* coefficient correspond to the strongest mixing and likely the minimum level of contamination at ground level. It is evident analysing fig. 4 that the maximum values are reached during the day as a consequence of the strong diurnal convective mixing and in the range of a dimensionless height between [0.4, 0.7], that is in the bulk of the convective boundary layer. During the night the mixing is reduced as a consequence of the formation of the stable boundary layer due to the inversion of the heat

Kz/(w\* zi)

0 10 20 30 40 50 60 70 80 90 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

In the further we show the radioactive substance concentrations close to the surface around the nuclear power plant. Figures 5 show the distributions for 3 hours, 48 hours and 93 hours

The centre of the nuclear power plant is located in the centre of the plot, the cost line is almost in the north south direction, that is parallel to the *y*-axis in the plot with the ocean to the right side. Shortly after the beginning of the release the mean wind pointed towards the ocean,

**Figure 4.** The dimensionless eddy diffusivity coefficient dependent on height in multiples of the

whereas after three days the wind blew towards the south in the direction of Tokyo.

z/zi

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

233

Time/[h]

boundary layer height for a time sequence of four subsequent days.

after the beginning of the substance release with a logarithmic scale.

93 hours after time zero.

 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18

flux.

In the following we show the results for a sequence of four days from the 12th to the 15th of march. Figure 3 shows some meso-scale meteorological information, that was obtained from WRF. The first plot in fig. 3 corresponds to the situation three hours after the beginning of

**Figure 3.** Temperature and mean wind profile from WRF for 3 hours, 48 hours and 93 hours after the beginning of constant release of radioactive substances.

constant release of radioactive material, the second and third plot correspond to 48 hours and 93 hours after time zero.

14 Will-be-set-by-IN-TECH

At the 11th of march, 2011 the Fukushima-Daiichi nuclear power plant accident (coordinates in latitude, longitude: 37*<sup>o</sup>* 25' 17" N, 141*<sup>o</sup>* 1' 57" E) caused considerable radiation leakage into the atmosphere and into the sea. The radioactive pollution of the environment and sea was caused principally by the direct release of contaminated water from the power station. To a lesser extent atmospheric release of the radio-nuclide from the atmospheric plume are carried by the winds over the sea during and after the accident sequence. Shorter–lived radioactive elements, such as Iodine-131 were detectable for a few months (half-live of approximately 80 days). Others, such as Ruthenium-106 and Caesium-134 will still persist in the environment

In the following we show the results for a sequence of four days from the 12th to the 15th of march. Figure 3 shows some meso-scale meteorological information, that was obtained from WRF. The first plot in fig. 3 corresponds to the situation three hours after the beginning of

**Figure 3.** Temperature and mean wind profile from WRF for 3 hours, 48 hours and 93 hours after the

beginning of constant release of radioactive substances.

for several years (Caesium-137 has a half-life of approximately 30 years).

From the meso-scale meteorological data one may determine the eddy diffusivity coefficients for each specific hour. In the *z*-time plot of fig. 4 we report the dimensionless vertical eddy diffusivity coefficients as calculated by eq. (35) for four subsequent days. The figure shows in a simple way the spatial and time structure of this coefficient. In this context, it is important to point out that the largest values of the *Kz* coefficient correspond to the strongest mixing and likely the minimum level of contamination at ground level. It is evident analysing fig. 4 that the maximum values are reached during the day as a consequence of the strong diurnal convective mixing and in the range of a dimensionless height between [0.4, 0.7], that is in the bulk of the convective boundary layer. During the night the mixing is reduced as a consequence of the formation of the stable boundary layer due to the inversion of the heat flux.

In the further we show the radioactive substance concentrations close to the surface around the nuclear power plant. Figures 5 show the distributions for 3 hours, 48 hours and 93 hours after the beginning of the substance release with a logarithmic scale.

The centre of the nuclear power plant is located in the centre of the plot, the cost line is almost in the north south direction, that is parallel to the *y*-axis in the plot with the ocean to the right side. Shortly after the beginning of the release the mean wind pointed towards the ocean, whereas after three days the wind blew towards the south in the direction of Tokyo.

of a full diurnal cycle is possible. Simulating micro-meteorology for a short period for the Fukushima Nuclear Power Station Accident may be considered a first step into a direction where the impact of the contamination of radioactive material in the site may be simulated and evaluated for the whole period of the accident until today. Thus the present work may be understood as one tile in a larger program development that simulates radioactive material dispersion using analytical resources, i.e. solutions. In a longer term we intend to build a library that allows to predict radioactive material transport in the planetary boundary layer

The quality of the solution may be estimated by the following considerations. Recalling, that the structure of the pollutant concentration is essentially determined by the mean wind

a length scale for which the pollutant concentration is almost homogeneous. Thus one may

the solution become spurious. Upon interpreting *�*−<sup>1</sup> as a sampling density, one may now employ the Cardinal Theorem of Interpolation Theory ([30]) in order to find the truncation that leaves the analytical solution almost exact, i.e. introduces only functions that vary

is bounded by *m�*−<sup>1</sup> has an exact solution for a finite expansion. This statement expresses

sort of sampling density, its introduction is an approximation and is related to convergence of the approach and Parseval's theorem may be used to estimate the error. In order to keep the solution error within a prescribed error, the expansion in the region of interest has to

the theorem the solution is then exact. In our approximation, if *m* is properly chosen such that

Further, the Cauchy-Kowalewski theorem ([8]) guarantees that the proposed solution is a valid solution of the discussed problem, since this problem is a special case of the afore mentioned theorem, so that existence and uniqueness are guaranteed. It remains to justify convergence of the decomposition method. In general convergence by the decomposition method is not guaranteed, so that the solution shall be tested by an appropriate criterion. Since standard convergence criteria do not apply in a straight forward manner for the present case, we resort to a method which is based on the reasoning of Lyapunov ([6]). While Lyapunov introduced this conception in order to test the influence of variations of the initial condition on the solution, we use a similar procedure to test the stability of convergence while starting from an approximate (initial) solution **R**<sup>0</sup> (the seed of the recursive scheme). Let

*<sup>i</sup>*=*n*+<sup>1</sup> **R***i*� be the maximum deviation of the correct from the approximate solution

*<sup>i</sup>*=<sup>0</sup> **R***i*, where �·� signifies the maximum norm. Then strong convergence occurs if

For model validation one faces the drawback, that the majority of measurements are at ground level, so that one could think that a two dimensional description would suffice, however the present analysis clearly shows the influence of the additional dimension. While in the two dimensional approach the tendency of the predicted concentrations is to overestimate

*<sup>m</sup>* and *m* an increasing integer number) variations in

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

*<sup>r</sup> c dt dx d* ¯ *<sup>η</sup>* <sup>∈</sup> *<sup>L</sup>*<sup>2</sup> (*<sup>η</sup>* <sup>=</sup> *<sup>y</sup>* or *<sup>z</sup>*) with spectrum {*λi*} which

for our problem. Since the cut-off defines some

. For the bounded spectrum and according to


235

�Γ*n*� log <sup>|</sup>*δZn*<sup>|</sup>


velocity **U**¯ and the eddy diffusivity **K**, means that the quotient of norms *�* = ||**K**||

that extends from the micro- to the meso-scale.

conclude that with decreasing length ( *�*

The square integrable function *χ* =

 of

contain *<sup>n</sup>* <sup>+</sup> 1 terms, with *<sup>n</sup>* <sup>=</sup> int *mLy*,*<sup>z</sup>*

the

Cardinal

<sup>|</sup>*δZn*<sup>|</sup> <sup>=</sup> � <sup>∑</sup><sup>∞</sup>

Γ*<sup>n</sup>* = ∑*<sup>n</sup>*

 Theorem

significantly in length scales beyond the mentioned limit.

Interpolation

 Theory

<sup>2</sup>*π�* <sup>+</sup> <sup>1</sup> 2 

there exists an *<sup>n</sup>*<sup>0</sup> such that the sign of *<sup>λ</sup>* is negative for all *<sup>n</sup>* <sup>≥</sup> *<sup>n</sup>*0. Here, *<sup>λ</sup>* <sup>=</sup> <sup>1</sup>

the cut-off part of the spectrum is negligible, then the found solution is almost exact.

**Figure 5.** The logarithmic concentration distribution of radioactive substances released from the nuclear power plant for 3 hours, 48 hours and 93 hours after the beginning of the release.
