Open access peer-reviewed chapter

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

Written By

Marco Túllio Vilhena, Bardo Bodmann, Umberto Rizza and Daniela Buske

Submitted: November 25th, 2011 Reviewed: July 18th, 2012 Published: October 10th, 2012

DOI: 10.5772/51525

Chapter metrics overview

1,328 Chapter Downloads

View Full Metrics

1. Introduction

While the renaissance of nuclear power was motivated by the increasing energy demand and the related climate problem, the recent history of nuclear power, more specifically two disastrous accidents have forced focus on nuclear safety. Although, experience gathered along nuclear reactor developments has sharpened the rules and regulations that lead to the commissioning of latest generation nuclear technology, an issue of crucial concern is the environmental monitoring around nuclear power plants. These measures consider principally the dispersion of radioactive material that either may be released in control actions or in accidents, where in the latter knowledge from simulations guide the planning of emergency actions. In this line the following contribution focuses on the question of radioactive material dispersion after discharge from a nuclear power plant.

The atmosphere is considered the principal vehicle by which radioactive materials that are either released from a nuclear power plant in experimental or eventually in accidental events could be dispersed in the environment and result in radiation exposure of plants, animals and last not least humans. Thus, the evaluation of airborne radioactive material transport in the atmosphere is one of the requirements for monitoring and planning safety measures in the environment around the nuclear power plant. In order to analyse the (possible) consequences of radioactive discharge atmospheric dispersion models are of need, which have to be tuned using specific meteorological parameters and conditions in the considered region. Moreover, they shall be subject to the local orography and supply with realistic information on radiological consequences of routine discharges and potential accidental releases of radioactive substances.

The present work provides a model that allows to implement afore mentioned simulations by the use of a hybrid system. In a first step the local meteorological parameters are determined using the next-generation mesoscale numerical weather prediction system “Weather Research 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 within a radius of several tenth kilometers.


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

tc-+Uc--TKc-=S  .E1

Here U=(u-,v-,w-)Tis the vector field of the mean wind velocity, the diagonal matrix K= diag (Kx,Ky,Kz)contains the eddy diffusivities and Sis 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)


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 by0<x<Lx, 0<y<Ly(with Lxand Lysufficiently large), 0<z<h(here his the boundary layer height) and subject to the following boundary and initial conditions,

c-(x,y,z,0)=0  .E4

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 rS=(0,y0,HS)is located at the boundary of the domainrSδΓ. 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)0forrΓ\δΓ), 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

S=ωS  dΣE5

where ωSrepresents 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

c-(t,x,y,z)=0tc-˙(t-τ,x,y,z)  dτE6

with instantaneous initial condition

c-˙(0,x,y,z)=c-˙0=limΣ^  dΣ0ωS  dΣ=Qδ(x)δ(y-y0)δ(z-HS)E7

where Qis the emission rate, HSthe 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 yvariable, and expand the pollutant concentration as


where R=(R1,R2,)Tand Y=(Y1,Y2,)Tis a vector in the space of orthogonal eigenfunctions, given by Ym(y)=cos(λmy)with eigenvalues λm=mπLyform=0,1,2,. For convenience we introduce some shorthand notations, 2=(x,0,y)Tand^y=(0,y,0)T, so that equation (1) reads now,

                        =2TK+(K2)T(2RTY)+^yTK+(K^y)T(RT^yY)  .E10

Upon application of the integral operator

0LydyY[F]=0LyFTY  dyE11

here Fis 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

B0=0LydyY[Y]=0LyYTY  dy  ,E12
Z=0LydyY^yY=0Ly^yYTY  dy,E13
Ω1=0LydyY[(2TK)(2RTY)]=0Ly(2TK)(2RTY)TY  dy  ,E14
Ω2=0LydyY[(K2)T(2RTY)]=0Ly(K2)T(2RTY)Y  dy  ,E15
T1=0LydyY[((^yTK)(^yY)]=0Ly((^yTK)(^yY)TY  dy  ,E16
T2=0LydyY[(K^y)T(^yY)]=0Ly(K^y)T(^yY)TY  dy  .E17

Here, B0=Ly2I, where Iis the identity, the elements Zmn=21-n2/m2δ1,jwith δi,jthe Kronecker symbol and j=(m+n) mod 2is the remainder of an integer division (i.e. is one for m+nodd and zero else). Note, that the integrals Ωiand Tidepend on the specific form of the eddy diffusivityK. The integrals (11) are general, but for practical purposes and for application to a case study we truncate the eigenfunction space and consider Mcomponents in Rand Yonly, though continue using the general nomenclature that remains valid. The obtained matrix equation determines now together with initial and boundary condition uniquely the components Rifor i=1,,Mfollowing the procedure introduced in reference [22]:


3.2. A specific case for application

In order to discuss a specific case we introduce a convention and consider the average wind velocity U-=(u-,0,0)Taligned 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 fort=0. By comparison of physically meaningful cases, one finds for the operator norm||xKxx||<<|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 KxandxKx.

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 KK1= diag (0,Ky,Kz)depends in leading order approximation only on the vertical coordinateK1=K1(z). For this specific case the integrals Ωireduce to

Ω1(zKz)(zRT)B  ,E19
Ω2Kz(z2RT)B  ,E20
T10  ,E21
T2-KyΛB  ,E22

whereΛ= diag (λ12,λ22,). The simplified equation system to be solved is then,


which is equivalent to the problem


by virtue of Bbeing 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 Kis assumed to be independent ofc-, whereas in more realistic cases, even if stationary, Kmay 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 Kwhich 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, 2, 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-z(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-z(z)+κz(z,t)and use the asymptotic form ofKy, which is then explored to set-up the structure of the equation that defines the recursive decomposition scheme:


The function R=jRj=1TR(c)is now decomposed into contributions to be determined by recursion. For convenience we introduced the one-vector 1=(1,1,)Tand inflate the vector Rto a vector with each element being itself a vectorRj. Upon inserting the expansion in equation (17) one may regroup terms that obey the recursive equations and starts with the time averaged solution forKz:


The extension to the closed form recursion is then given by

tRj+u-xRj-zK-zzRj+KyΛRj=zκzzRj-1  .E27

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 Kwas 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 R0satisfies 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 xdirection 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 tvariable, (tr) obtaining the following pseudo-steady-state problem:

rR~0+u¯xR~0=zKzzR~0-Λ KyR~0E28
The xand zdependence 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:

where C=(ζ1(z),ζ2(z),)Tare a set of orthogonal eigenfunctions, given byζi(z)=cos(γlz), and γi=iπ/h(fori=0,1,2,) are the set of eigenvalues.

Replacing equation (21) in equation (20) and using the afore introduced projector (10) now for the zdependent degrees of freedom 0hdzC[F]=0hFTC  dzyields a first order differential equation system:

xP+HP=0  ,E30

where P=P(x,r)andH=B1-1B2. The entries of matrices B1and B2are

(B1)i,j=-0hu¯ζi(z)ζj(z)  dzE31
(B2)i,j=0hzKzzζi(z)ζj(z)  dz-γi20hKzζi(z)ζj(z)  dzE32
-r0hζi(z)ζj(z)  dz-λi2Ky0hζi(z)ζj(z)  dz  .E33

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


Following the reasoning of [22] we solve (22) applying Laplace transform and diagonalisation of the matrix H=XDX-1which results in


where P~(s,r)denotes the Laplace Transform ofP(x,r). Here X(-1)is the (inverse) matrix of the eigenvectors of matrix B1-1B2with diagonal eigenvalue matrix Dand the entries of matrix(sI+D)ii=s+di. After performing the Laplace transform inversion of equation (24), we come out with

P(x,r)=XG(x,r)X-1Ω  ,E36

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-1P(0,r).

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

R0(x,z,t)=12πiγ-iγ+iP(x,r)C(z)ert  dr  .E37

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 2M-1in the 1rvariable

R0(x,z,t)=1taTpR0(x,z,pt)  ,E38

where aand pare respectively vectors with the weights and roots of the Gaussian quadrature scheme ([27]), and the argument (x,z,pt)signifies the k-th component of pin the k-th row ofpR0. Note, kis a component from contraction witha.


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 100mheight 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 (10m,60mand100m) together with the temperature gradients between 10mand100m. 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 1mand the Monin-Obukhov length for convective conditions can be written as L=-h/ku*/w*3([35]), where kis the von Karman constant (k=0.4), w*is the convective velocity scale with wind speedU, u*=kU/ln(zr/z0)is the friction velocity, where Uis the wind velocity at the reference heightzr=10m, and h=0.3u*/fcis the height of the boundary layer with the Coriolis coefficientfc=10-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]):

Ky=πσv16(fm)vqv         with         σv2=0.98cv(fm)v23ψεqv23zh23w*2E40
qv=4.16zh  ,        ψε13=1-zh2-zL-23+0.7512         and         (fm)v=0.16E41

where σvis the standard deviation of the longitudinal turbulent velocity component, qvis the stability function, ψεis the dimensionless molecular dissipation rate and (fm)vis the transverse wave peak.

The wind speed profile can be described by a power law uz/u1=z/z1n([25]), where uzand u1are the horizontal mean wind speeds at heights zand z1and nis 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.

Exp.PeriodDistance (m)Observed (Bq/m3)Predictions (Bq/m3)

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.

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.

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 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 oand prefer 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.

Statistical IndicesGILTTG3D-GILTT
NMSE =(Co-Cp)2¯/Cp¯ Co¯1.340.38
COR =(Co-Co¯)(Cp-Cp¯)¯/σoσp0.670.83
FA2 =0.5(Cp/Co)20.530.88
FA5 =0.2(Cp/Co)50.961.00
FB =Co¯-Cp¯/0.5(Co¯+Cp¯)-0.440.13
FS =(σo-σp)/0.5(σo+σp)-0.540.18

Table 3.

Statistical comparisons between GILTTG and 3D-GILTT results.

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-p=1.16C-o+7.01with R2=0.67andκ=0.43, whereas the 3D-GILTT obeys the result C-p=0.69C-o+3.26with R2=0.83andκ=0.36. In order to perform a model validation we introduced an index κ=(a-1)2+b/C-o2withC-o=1ni=1nCoi, which if identical zero indicates a perfect match between the model and the experimental findings. Here ais the slope, bthe intersection, Coiof the experimental data and C-oits 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 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.

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


where α=(x,y,z)andi=u,v,w, FiE(n)is the value of the Eulerian spectrum of energy normalized by the Eulerian velocity variance, and σi2corresponds to the Eulerian variance of the turbulent wind field. Following [33],βi=πU216σi212. For large diffusion travel times (t), the filter function in the integral of eqn. (30) selects FiE(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,


where FiE(0)is the value of the normalised Eulerian energy spectrum atn=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 asSiE(n)=Sib(n)+Sis(n), where the subscripts band sstand for buoyancy and shear, respectively.

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

Sib(n)=0.98cinzu-nfmi*531+1.5nzu-fmi*Ψεb23zzi23w*2  ,E44

while for mechanical turbulence ([10])


where Ψεb=εbhw*3and Φεs=εκzu*3are the dimensional dissipation rate functions, εband εsare the convective and mechanical rate of tkedissipation, fmi*is the normalized frequency of the spectral peaks regardless of stratification and fmiis the reduced frequency with the mean wind speed u-in the mixing layer.

The dimensionless spectrum FiE(n)in eqn. (31) is obtained by normalizing the dimensional spectra with the total variance, σi2=0SiE(n)dn, that is

FiE(n)=SiEσi2=SibE(n)+SisE(n)σi2  .E46

The total wind velocity variance is obtained by the sum of mechanical and convective variancesσi2=0(SibE(n)+SisE(n))  dn=σib2+σis2. Making use of eqns. (30), (32), (33) and eqn. (34) one ends up withKα=βi4SibE(0)+SisE(0), that for the w-component becomes


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 11  th of march, 2011 the Fukushima-Daiichi nuclear power plant accident (coordinates in latitude, longitude: 37 o25’ 17” N, 141 o1’ 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).

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.

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 constant release of radioactive material, the second and third plot correspond to 48 hours and 93 hours after time zero.

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 Kzcoefficient 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.

Figure 4.

The dimensionless eddy diffusivity coefficient dependent on height in multiples of the boundary layer height for a time sequence of four subsequent days.

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.

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.

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.


7. Conclusions

The present work was based on an Eulerian approach to determine dispersion of radioactive contaminants in the PBL. To this end the diffusion equation for the cross-wind integrated concentrations was closed by the relation of the turbulent fluxes to the gradient of the mean concentration by means of eddy diffusivity (K-theory). We are completely aware of the fact that K-closure has its intrinsic limits so that one would like to remove these inconsistencies. However, comparisons of predictions by this approach to experimental data have shown that there are scenarios where this lack is not significantly manifest, which we use as a justification together with its computational simplicity to perform our simulations based on this approach.

Since the consistency of the K-approach depends crucially on the determination of the eddy diffusivity considering the turbulence structure of the PBL in its respective stability regimes, we elaborated parametrisations for the eddy diffusivity coefficients based on the micro-meteorological parameters that were extracted from meso-scale WRF simulations, that allowed to take into account the realistic orography of the larger vicinity of a reactor site in consideration. The approach proposed here for the determination of the eddy-diffusivity coefficient is based on the Taylor statistical diffusion theory and on the spectral properties of turbulence. The assumption of continuous turbulence spectrum and variances, allows the parametrisations to be continuous at all elevations, and in stability conditions ranging from a convective to a neutral condition, and from a neutral to a stable condition so that a simulation 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 that extends from the micro- to the meso-scale.

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 velocity U-and the eddy diffusivityK, means that the quotient of norms ϖ=||K||||U-||defines a length scale for which the pollutant concentration is almost homogeneous. Thus one may conclude that with decreasing length (ϖmand man increasing integer number) variations in the solution become spurious. Upon interpreting ϖ-1as 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 significantly in length scales beyond the mentioned limit.

The square integrable function χ=rc-  dt dx dηL2(η=yorz) with spectrum {λi}which is bounded by mϖ-1has an exact solution for a finite expansion. This statement expresses the Cardinal Theorem of Interpolation Theoryfor our problem. Since the cut-off defines some 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 contain n+1terms, withn= int mLy,z2πϖ+12. For the bounded spectrum and according to the theorem the solution is then exact. In our approximation, if mis properly chosen such that the cut-off part of the spectrum is negligible, then the found solution is almost exact.

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 R0(the seed of the recursive scheme). Let |δZn|=i=n+1Ribe the maximum deviation of the correct from the approximate solutionΓn=i=0nRi, where signifies the maximum norm. Then strong convergence occurs if there exists an n0such that the sign of λis negative for allnn0. Here,λ=1Γnlog|δZn||δZ0|.

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 the observed values, this is not the case for the results of the three dimensional description, mainly because it does not assume turbulence to be homogeneous. Moreover the solution of the advection diffusion equation discussed here is more general than shown in the present context, so that a wider range of applications is possible. Especially other assumptions for the velocity field and the diffusion matrix are possible. In a future work we will focus on a variety of applications and introduce a rigorous proof of convergence from a mathematical point of view, which we indicated in sketched form only in our conclusions.


The authors thank Brazilian CNPq and FAPERGS for the partial financial support of thiswork.


  1. 1. AdomianG.1984A New Approach to Nonlinear Partial Differential EquationsJ. Math. Anal. Appl.,102page numbers (420 EOF434 EOF
  2. 2. AdomianG.1988A Review of the DecompositionMethod in AppliedMathematics. J. Math. Anal. Appl.,135page numbers (501-544).
  3. 3. AdomianG.1994Solving Frontier Problems of Physics: The Decomposition Method,Kluwer, Boston, MA.
  4. 4. BatchelorG. K.1950The application of the similarity theory of turbulence to atmospheric diffusionQuart. J. Royal Meteor. Soc.,76328page numbers (133 EOF146 EOF
  5. 5. BiagioR.GodoyG.NicoliI.NicolliD.ThomasP.1985First atmospheric diffusion experiment campaign at the Angra site.- KfK 3936, Karlsruhe, and CNEN 1201, Rio de Janeiro.
  6. 6. BoichenkoV. A.LeonovG. A.ReitmannV.2005NDimension theory for ordinary equations, Teubner, Stuttgart.
  7. 7. CostaC. P.VilhenaM. T.MoreiraD. M.TirabassiT.2006Semi-analytical solution of the steady three-dimensional advection-diffusion equation in the planetary boundary layerAtmos. Environ.,4029page numbers (5659 EOF
  8. 8. CourantR.HilbertD.(19891989Methods ofMathematical Physics. JohnWiley&Sons, New York.
  9. 9. DeardoffJ. W. .WillisG. E.1975A parameterization of diffusion into the mixed layer. J. Applied Meteor.,14page numbers (1451-1458).
  10. 10. DegraziaG. A.MoraesO. L. L.1992A model for eddy diffusivity in a stable boundary layerBound. Layer Meteor.,58page numbers (205 EOF214 EOF
  11. 11. DegraziaG. A.CamposVelho. H. F.CarvalhoJ. C.1997Nonlocal exchange coefficients for the convective boundary layer derived from spectral propertiesContr. Atmos. Phys.,70page numbers (57 EOF64 EOF
  12. 12. DemuthC.1978A contribution to the analytical steady solution of the diffusion equation for line sourcesAtmos. Environ.,12page numbers (1255 EOF1258 EOF
  13. 13. DjolovG. D.YordanovD. L.SyrakovD. E.2004Baroclinic planetary boundary layer model for neutral and stable stratification conditionsBound. Layer Meteor.,111page numbers (467 EOF490 EOF
  14. 14. HannaS. R.1989Confidence limit for air quality models as estimated by bootstrap and jacknife resampling methods. Atmos. Environ.,23page numbers (1385-1395).
  15. 15. HjstrupJ. H.1982Velocity spectra in the unstable boundary layer. J. Atmos. Sci.,39page numbers (2239-2248).
  16. 16. IrwinJ. S.1979A theoretical variation of the wind profile power-low exponent as a function of surface roughness and stability. Atmos. Environ.,13page numbers (191-194).
  17. 17. MangiaC.DegraziaG. A.RizzaU.2000An integral formulation for the dispersion parameters in a shear/buoyancy driven planetary boundary layer for use in a Gaussian model for tall stacksJ. Applied Meteor.,39page numbers (1913 EOF1922 EOF
  18. 18. MangiaC.MoreiraD. M.SchipaI.DegraziaG. A.TirabassiT.RizzaU.2002Evaluation of a new eddy diffusivity parameterisation fromturbulent eulerian spectra in different stability conditions. Atmos. Environ.,3634page numbers (67-76).
  19. 19. MellorG. L.YamadaT.1982Development of a Turbulence Closure Model for Geophysical Fluid ProblemsReviews of Geo. and Space Phys.,20page numbers (851 EOF
  20. 20. MoengC. H.SullivanP. P.1994A comparison of shear-and buoyancy-driven planetary boundary layer flows.J. Atmos. Sci.,51page numbers (999 EOF1022 EOF
  21. 21. MoreiraD. M.VilhenaM. T.TirabassiT.CostaC.BodmannB.2006Simulation of pollutant dispersion in atmosphere by the Laplace transform: the ADMM approach.Water, Air and Soil PollutionAir and Soil Pollution,177page numbers (411 EOF439 EOF
  22. 22. MoreiraD. M.VilhenaM. T.BuskeD.TirabassiT.2009The state-of-art of the GILTT method to simulate pollutant dispersion in the atmosphereAtmos. Research,92page numbers (1 EOF17 EOF
  23. 23. Nieuwstadt F.T.M. & de Haan B.J.1981An analytical solution of one-dimensional diffusion equation in a nonstationary boundary layer with an application to inversion rise fumigation. Atmos. Environ.,15page numbers (845-851).
  24. 24. OlesenH. R.LarsenS. E.HøjstrupJ.1984Modelling velocity spectra in the lower part of the planetary boundary layer. Bound. Layer Meteor.,29page numbers (285-312).
  25. 25. PanofskyA. H.DuttonJ. A.1988Atmospheric Turbulence. John Wiley & Sons, New York.
  26. 26. SharanM.SinghM. P.YadavA. K.1996Amathematical model for the atmospheric dispersion in low winds with eddy diffusivities as linear functions of downwind distance. Atmos. Environ.,307page numbers (1137 EOF
  27. 27. StroudA. H.SecrestD.1966Gaussian quadrature formulasPrentice Hall Inc., Englewood Cliffs, N.J..
  28. 28. TaylorG. I.1921Diffusion by continuous movement. Proc. Lond.Math. Soc.,2page numbers (196-211).
  29. 29. TirabassiT.2003Operational advanced air pollution modelingPAGEOPH,1601-2page numbers (05 EOF16 EOF
  30. 30. TorresR. H.1991Spaces of sequences, sampling theorem, and functions of exponential type. Studia Mathematica,1001page numbers (51-74).
  31. 31. UlkeA. G.2000New turbulent parameterisation for a dispersion model in the atmospheric boundary layer. Atmos. Environ.,34page numbers (1029-1042).
  32. 32. van UldenA. P.1978Simple estimates for vertical diffusion from sources near the ground. Atmos. Environ.,12page numbers (2125-2129).
  33. 33. WandelC. F.Kofoed-HansenO.1962On the Eulerian-Lagrangian Transform in the Statistical Theory of TurbulenceJ. Geo. Research,67page numbers (3089 EOF3093 EOF
  34. 34. YordanovD.SyrakovD.KolarovaM.1997On the Parameterization of the Planetary Boundary Layer of the Atmosphere: The Determination of the Mixing Height, In: Current Progress and Problems. EURASAPWorkshop Proc..
  35. 35. ZanettiP.(1990). Air Pollution Modeling. Comp. Mech. Publications, Southampton (UK).

Written By

Marco Túllio Vilhena, Bardo Bodmann, Umberto Rizza and Daniela Buske

Submitted: November 25th, 2011 Reviewed: July 18th, 2012 Published: October 10th, 2012