Prediction of the 137Cs activity concentration in the atmospheric surface layer of the Chernobyl exclusion zone

Journal of Environmental Radioactivity 110 (2012) 53e58

Contents lists available at SciVerse ScienceDirect

Journal of Environmental Radioactivity journal homepage: www.elsevier.com/locate/jenvrad

Prediction of the 137Cs activity concentration in the atmospheric surface layer of the Chernobyl exclusion zone E.K. Garger a, Yu.I. Kuzmenko a, S. Sickinger b, J. Tschiersch b, * a b

National Academy of Science of Ukraine, Institute for Safety Problems of Nuclear Power Plants, Kiyv 03028, Ukraine Helmholtz Zentrum München, German Research Center for Environmental Health, Institute of Radiation Protection, Ingolstädter Landstr. 1, 85764 Neuherberg, Germany

a r t i c l e i n f o

a b s t r a c t

Article history: Received 25 February 2011 Received in revised form 18 January 2012 Accepted 18 January 2012 Available online 21 February 2012

The time series of the 10-day average 137Cs volumetric activity concentration in the lower atmosphere measured from 1987 to 1991 in the town of Pripyat, close to the Chernobyl nuclear power plant, was used to construct a model to predict the airborne activity concentration inside the 30-km exclusion zone. For that purpose, individual components of the observed time series were separated by regression analysis and the Group Method of Data Handling. The measured data in Pripyat were divided in two periods. The long-term prediction by the model established using the measured data of the first period, has been validated with the data in the second period with good agreement. The behaviour of the model parameters depending on the length of the periods was also analysed, and the first period of 4.5 y was shown as sufficient for estimating the parameters. Further increase in the length will not significantly enhance the model parameters and the predictive power. Ó 2012 Elsevier Ltd. All rights reserved.

Keywords: Atmospheric radioactivity 137 Cs Model development Prognosis Validation

1. Introduction After the termination of the primary release of radionuclides from the destroyed Chernobyl nuclear reactor into the atmosphere, a secondary radioactive atmospheric aerosol was detected in the Chernobyl area (Garger et al., 1998). Strong winds mobilized dust particles from contaminated surfaces such as cultivated fields, meadows and forest. Human activities on fields, roads and construction sites even increased the resuspension of contaminated dust. Additional sources of airborne radioactive particles were forest fires (Holländer and Garger, 1996; Yoschenko et al., 2006; Paatero et al., 2009) and permanent but unsteady emissions from openings of the Chernobyl Shelter (Garger et al., 2006). In this context it is necessary to estimate the resulting inhalation dose, its variation in time and the horizontal redistribution of radioactivity due to migration of radioactive particles in the 30-km exclusion zone and the surrounding area. The knowledge of the activity concentration of the doserelevant radionuclides is required at first for estimation of inhalation dose and related analysis, e.g. the horizontal redistribution of the radionuclides. Since 137Cs is a major radionuclide released from the Chernobyl accident with relatively long half-life, information of the soil contamination density and long-time series data of its

* Corresponding author. Tel.: þ49 89 31872763; fax: þ49 89 31873363. E-mail address: [email protected] (J. Tschiersch). 0265-931X/$ e see front matter Ó 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.jenvrad.2012.01.017

atmospheric concentration have been accumulated for inside and around the exclusion zone. In Garger et al. (1997) prospects and limitations are shown of the approach to forecast the activity concentration in air by using the soil contamination density and the resuspension factor. In particular, Garger et al. (1995) have found that the models based on the resuspension factor does not take into account the fluctuations of the airborne activity concentration, and the predicted values may therefore differ from the measured values several orders of magnitude. In the present work an alternative statistical approach is presented. It uses the measurement of the 137Cs activity concentration in air to construct a long-term predictive model. The aims of this work are i) to report about establishing the predictive model and forecasting the 137Cs air activity based on data measured in the town of Pripyat in the period 1987e1991 and ii) to validate the predicted values with measurement data in the same town in the later period 1991e1997. 2. Material and methods 2.1. Data acquisition Filter samples of aerosol were taken at five locations (Table 1) in the Ukrainian Polesia region (including Pripyat and locations up to 230 km distance from the Chernobyl nuclear power plant, see Fig. 1) from samplers designed by SPA “Typhoon”, Russia, operated

54

E.K. Garger et al. / Journal of Environmental Radioactivity 110 (2012) 53e58

Table 1 Collection period of the filter samples and the averaged sample period for the five sites in the Ukrainian Polesia region.

Kiev

Baryshevka

Chernobyl

Polesskoe Pripyat

Collection period

Averaged sample period (d)

Source

July 1987eDecember 1992 January 1993eDecember 2002 July 1987eAugust 1994 September 1994eDecember 2002 July 1987eDecember 1991 January 1992eDecember 2002 July 1987eMarch 1991 July 1987eDecember 1991 January 1992eAugust 2005

30

Garger et al. (1999)

and compiled for further analysis. In Fig. 2 the 10-year time series of airborne 137Cs activity concentration in the town of Pripyat is presented. Table 1 shows the collection period of the filter samples at the five sites in the Ukrainian Polesia region.

2.2. Data processing

This work 30

Garger et al. (1997) This work

30

Garger et al. (1994) This work

30 10

Garger et al. (1999) Garger et al. (1994) This work

at a flow rate of 100 km3 d1 (1.16 m3 s1) at a sampling height of about 1.5 m above lawn. One sampler was situated in the town of Pripyat in an evacuated residential area in the middle of a large lawn. It must be noticed that the measurement device was relocated closer to the nuclear power plant in 1998. In Pripyat the filters were changed daily, at other places the period of sampling was variable depending on the level of airborne activity concentration (1e10 d). The exposed filters (Russian type Petryanov cloth FPP-1515 with an area of 1.05 m2) were pressed into tablets and analysed by g-spectrometry for 137Cs. The g-spectrometer based on the semi-conductor detector GMX-30190 (relative efficiency 32.5%, energy resolution 1.89 keV on the line 1.33 MeV) and the multichannel buffer Spectrum Master (ORTEC Inc.). Energy and efficiency calibration were performed with the standard source PSGJ No. 13166, traceable to Russian Standards. The measured values of the 137Cs concentration were averaged over consecutive 10 d for Pripyat and consecutive 30 d (monthly average) for the other sites,

The Pripyat data demonstrates the significant temporal variation of the airborne 137Cs activity concentration in the aftermath of the Chernobyl accident (Fig. 2). Therefore, the time series are considered as the result of dynamic processes that include different components such as trend and low-frequency fluctuations and other stochastic fluctuations. A first attempt to calculate a long-time prediction from the Pripyat air concentration data was made by Viswanathan et al. (2000). A technique specially adjusted for a high portion of missing data (about 7%) was applied for the calculations based on the least-square power spectrum analysis of Lomb (Lomb, 1976). To proceed further, an approach is applied of consequent selecting the time trends to estimate the single constituents of the variability of the 137Cs activity concentration. It consists of the following three steps: 1) Identification and reduction of a descending trend, 2) Detection of a periodic component and 3) Calculation of the stochastic components for the prediction. Further analysis is based on the 137Cs time series C(t) collected in Pripyat from 1987 to 1991. 2.2.1. Identification and reduction of descending trend The descending trend of the 137Cs concentration without any periodic component was analysed at first, and the results are shown in this section. The 10-day averages of the 137Cs concentration are used to reduce the influence of measuring errors and high-frequency fluctuations. After some empirical dependence

Fig. 1. The Ukrainian Polesia region with Kiev in the center and the location of the five measurement stations.

55

The GMDH algorithms for modelling polyharmonic processes and fields (Madala and Ivachnenko, 1993; Ivachnenko and Müller, 1984; Yurachkovsky and Popkov, 1986) are widely used and work especially well for short time series with a minimum of a priori information on the modelling process. These algorithms rest upon the idea of the inductive self-organization principle and choice of the best model depending on some external discrimination criterion. During the modelling procedure, the GMDH algorithm involves four heuristics that represent the main features of the GMDH theory:

1.E+05

3

Cs airborne activity density (μBq/m )

E.K. Garger et al. / Journal of Environmental Radioactivity 110 (2012) 53e58

1.E+04

1.E+03

1.E+02

137

1.E+01

1.E+00 1987

1989

1991

1993

1995

1997

Year 137

Fig. 2. Activity concentration of airborne Cs (10-day average) measured in the town of Pripyat. Data until 1991 is from Garger et al. (1994).

functions were tested using regression analysis, the following function was selected as the best one:

C0 ðtÞ ¼ aexpðt=sÞ þ b

(1)

where C0(t), estimated value of 10 d averaged a, b, coefficients, 1/s, decreasing rate coefficient.

137

Cs concentration,

Using the weighted least-squares technique, the coefficients a, s and b were determined to be a ¼ 6725 mBq m3, s ¼ 300 d and b ¼ 346 mBq m3. Then the 137Cs time series of measured C(t) were detrended by dividing by the estimated C0(t), and a dimensionless time series was obtained that represents the relative fluctuations around a descending trend. For further calculations u(t) ¼ log [C(t)/C0(t)] was used (Fig. 3). 2.2.2. Detection of periodic component To obtain a periodic component a method commonly known as the Group Method of Data Handling (GMDH) (Madala and Ivachnenko, 1993) was then applied to u(t) educed above. The Group Method of Data Handling is a set of several algorithms used to solve different problems. It consists of parametric, clusterization, analogues complexing, rebinarization and probability algorithms. This inductive approach is based on sorting-out of gradually complicated models and selecting the best solution based on an external criterion.

In the change of the volumetric activity concentration in the atmospheric surface layer, as well as in most multifactor processes in nature, exact periodicity is absent. Therefore, for the identification of a harmonic trend of the oscillatory process u(t) ¼ un and for the construction of models for long-term forecasting a model structure was introduced in the form of the sum of periodic functions with non-multiple frequencies uk sku1 ; not necessarily related to the period of observation (Ivachnenko, 1980).

un ¼ A0 þ

m P

½Ak sinðuk nÞ þ Bk cosðuk nÞ þ εn ;

k¼1

0 < uk < p;

(2)

k ¼ 1; 2; .; m;

where Ak and Bk are coefficients, minimizing the sum of squares of PN PN 2 * 2 n ¼ 1 εn ¼ n ¼ 1 ðun  un Þ ;

errors

with u*n , computational results of the model, m, present maximum number of harmonic components, N, number of 10-day in the examined interval of time. Unlike the methods of the fast Fourier transform, in this case the frequencies uk are defined from the values of cosu1, cosu2, ., cosum, determined as roots of the algebraic equation

a0 þ

1

Detrended time series u(t)

(1) Collect a set of observations that seems to be relevant to the object. (2) Divide the observations into two groups. The first will be used to estimate the coefficients of the model while the second will separate the information embedded in the data into either useful or harmful. (3) Create a set of elementary functions where its complexity will increase through an iterative procedure producing different models. (4) According to Godel’s incompleteness theorem, apply an external criterion to choose the optimal model.

m1 X

ap cosðpuÞ ¼ cosðmuÞ

(3)

p¼1 0.5

with coefficients ap which must satisfy the system of linear equations

0

m1 P p¼1

-0.5



ði ¼ 1; 2; .; N  2mÞ: (4)

-1

-1.5 1987



ap uiþp1 þ u2mþip1 þ am umþi1 ¼ ui1 þ u2mþi1

1989

1991

1993

1995

1997

Year Fig. 3. Detrended time series u(t) ¼ log [C(t)/C0(t)] of the data given in Fig. 2.

The GMDH harmonic algorithm employs the following scheme. At first, for each fixed m the harmonic components are computed with frequencies uk. The model parameters Ao, Ak, and Bk are estimated for each combination by using the least-squares technique. Then the obtained models are compared by means of the

56

E.K. Garger et al. / Journal of Environmental Radioactivity 110 (2012) 53e58

discrimination criterion. The best model is selected according to the value of the external criterion. The value of the external criterion in GMDH is always calculated on the sample data which has not been used for estimating the model coefficients, so called “test set”, unlike the “training set” used for calculating the coefficients. As an external criterion we used the criterion of minimizing the mean square error SB, calculated on the test set u(t)B.

SB ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ffi PNB  * u  u Bi i¼1 Bi NB

;

(5)

including the NB points of the row explored, not used for estimating the coefficients Ak and Bk in Eq. (2) and consisting of about 30% of the total number of available data points N. For the computation in this section, Fortran-95 was used with a working time usually running from 5 s to 2 min. Statistica 5.0 and Microsoft Excel 2003 were adopted for the other calculations. 2.2.3. Calculation of stochastic components The analysis of the deviations εi ¼ ui  u*i showed, that they are poorly correlated and normally distributed (Fig. 4). Therefore, x ¼ mε þ sεx, was calculated as stochastic component where mε and sε are mean and standard deviation from ε(t), and x w N(0,1) is a random variable following a standard normal distribution. Usually, m ε are very close to zero. The length of the sequence x is equal to the number of 10-day periods in the predictive period.

3. Results and discussion The coefficients of the model obtained using the GMDH on the basis of the trigonometric polynomial Eq. (2) are represented in Table 2. Based on the structure of the model, in Fig. 5 estimations of the spectral density P of the examined row u(t) are displayed in a periodogram as semi-log plot of the function Pðuk Þ ¼ N=2ðA2k þ B2k Þ, where Pðuk Þ is a function of the period 2p=uk . It can be seen that the most relevant contribution to Eq. (2) is by annual low-frequency fluctuations with a period of about 360 d (one year). The annual cycle may be induced by the annual increase in the concentration of 137Cs in April after the withdrawal of the snow cover. Furthermore, there is another peak at a period approximately equal to 70 d indicating seasonal fluctuations. Longterm recurrent seasonal fluctuations are characterized e.g. by a change of the mean wind speed, reduced surface roughness due to the lack of a vegetation cover and dryness/humidity of the soil

Table 2 Coefficients of the model computed using the Group Method of Data Handling (GMDH). Frequency (radian) uk

Period (d) Tk

Amplitude q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A2k þ B2k

Coefficient Ak

Coefficient Bk

Constant term A0

0.174 0.241 0.329 0.528 0.935 1.272 1.748 1.935 2.594

362 261 191 119 67 49 36 32 24

0.237 0.085 0.045 0.055 0.094 0.036 0.025 0.022 0.034

0.048 0.083 0.028 0.052 0.041 0.014 0.012 0.016 0.024

0.232 0.018 0.036 0.018 0.084 0.033 0.022 0.015 0.025

0.112

surface attributable to the local climate. The cycle of 24 d may be associated with human activity at the Shelter site which responds to a particular shift system and the periodicity in the change of the wind direction. Since the method aims to identify hidden periodicities in time series, not all influencing factors or combinations of them may find a clear and unambiguous explanation. The final expression for a long-term prediction of C(t) of the volumetric activity concentration of 137Cs can be written as

CðtÞ* ¼ 10zðtÞ C0 ðtÞ

(6)

Here C0(t) is a descending exponential fit with the constant term calculated in subsection 2.2.1, while z(t) ¼ u*(t) þ x(t) is a random cyclic component, which has the same low-frequency periods of fluctuation, variance and mean value as u(t), with u*(t) and x(t) calculated in subsection 2.2.2. and 2.2.3., respectively. In Fig. 6 the experimental data at Pripyat are shown for the whole measurement period from 1987 till 2005. The descending trend of the volumetric activity can be matched only by a model using a constant term. Finally, the predicted values of the airborne 137 Cs activity concentration as shown in Fig. 6 were computed using Eq. (7)

CðtÞ* ¼ 10zðtÞ ½6725expðt=300Þ þ 346;

(7)

where C(t)* is the required modelling function. The model predictions match the observed volumetric activity concentration of 137Cs pretty well. Two major deviations to the predicted values can be explained by unexpected phenomena; intense forest fires happened in 2002, and a significant release of radionuclides from the Chernobyl Shelter occurred in 2005.

10.0

45

362

Power spectrum P(T)

40

Frequency (%)

35 30 25 20 15

1.0 67

261 119 191

49

24

0.1

36 32

10 5

0.0 0

0

-1.1

-0.9

-0.7

-0.5

-0.3

-0.1

0.1

0.3

0.5

0.7

0.9

1.1

Deviation ε Fig. 4. Histogram of the frequency distribution of the deviations εi ¼ ui  u*i for the measurements at Pripyat.

50

100

150

200

250

300

350

400

450

Tk = 2 π / ω k (d)

Fig. 5. Estimation of the power spectrum P(T) of the time series u(t) of the measurements at Pripyat with GMDH, Tk ¼ 2p=uk (see text). Number in figure shows Tk in day.

E.K. Garger et al. / Journal of Environmental Radioactivity 110 (2012) 53e58 Data Predicted C 6725×EXP(-t/300)+346 C0(t)= EXP(-t/300)+346 0(t) = 6725

3

Cs airborne activity density (μBq/m )

1.E+05

1.E+04

Table 3 Correlation coefficients between the different sampling sites calculated from the mean monthly values (139 values) of the volumetric activity concentrations of 137Cs in the surface air.

Kiev Baryshevka Chernobyl Polesskoe Pripyat

1.E+03

1.E+02

a

137

1.E+01

1.E+00 1987

1989

1991

1993

1995

Fig. 6. Measured and predicted airborne 2005 (10-day average).

137

1997

1999

2001

2003

2005 Year

Cs activity concentration for Pripyat until

Additionally, the behaviour of the models parameters was studied in dependence on the length of the time series. It was observed that the model parameters began to stabilize when the time series of the 137Cs volumetric activity concentration reached a length of about 4.5 y, which corresponds to a time series beginning in 1987 and ending in the middle of 1991. A longer time series could not improve the model significantly and resulted in a variation of the parameters a and s of about 1e3%. At the same time, the change in the constant term b did not exceed 10%. This confirms the reliability of the determined parameters of the model, obtained for the 1987e1991 period, and enhances the confidence of its use for long-term prediction, relying on the estimates from the first 4e5 y after the accident. Furthermore, one is interested in the prediction of the volumetric activity concentration of 137Cs in different locations close to, but outside the exclusion zone. In Fig. 7 the temporal course of the measured activities concentrations of 137Cs at five places in the Ukrainian Polesia is shown. Three of the places (Polesskoe, Kiev, Baryshevka) are all located between 60 and 230 km away from the Chernobyl nuclear power plant (Fig. 1). It can be seen from Fig. 7 that the phase of stabilization of the airborne activity concentration has begun at around 1995 or 1996 for all sites. Table 3 shows the coefficients of correlation of the data measured at the five different sites. The measurements show a high positive correlation. Therefore it is concluded that the volumetric activity concentration of 137Cs at the five places is influenced by the same synoptic processes (wind speed, soil humidity, snow cover, etc.) and hence forecast of the airborne activity concentration can be based on the same prediction method using the observations

Fig. 7. The temporal course of the airborne 137Cs activity concentration at five places in the Polesia region.

57

Kiev

Baryshevka

Chernobyl

Polesskoe

Pripyata

1.0

0.70 1.0

0.71 0.35 1.0

0.82 0.34 0.79 1.0

0.54 0.43 0.70 0.70 1.0

128 values.

from these sites. However, to formulate such a forecast will be the task of a further study. 4. Conclusions A model for long-term forecasting of the activity concentration of 137Cs in the atmospheric surface layer from a limited number of observations in the Pripyat area was developed. The model is based on the identification of individual components of the non stationary process, such as descending and harmonic trends, and stochastic fluctuations. It takes into account the influence of annual and seasonal variations due to basic meteorological processes. According to the sampling and measurements performed from 1987 to 1991, the forecast appeared satisfactory both on average value of the volumetric activity concentration and on values of the amplitudes of fluctuations. It could be shown that the length of the sampling period was sufficient enough for setting up a reliable forecasting model. In general, magnitude and character of fluctuations of the forecasted volumetric activity concentration are similar to those that were measured in the Pripyat area from 1991 to 1997 (validation period). References Garger, E.K., Hoffman, F.O., Thiessen, K.M., 1997. Uncertainty of the long-term resuspension factor. Atmospheric Environment 31, 1647e1656. Garger, E.K., Hoffman, F.O., Thiessen, K.M., Galeriu, D., Kryshev, A.I., Lev, T., Miller, C.W., Nair, S.K., Talerko, N., Watkins, B., 1999. Test of existing mathematical models for atmospheric resuspension of radionuclides. Journal of Environmental Radioactivity 42, 157e175. Garger, E.K., Kashpur, V.A., Gurgula, B.I., Paretzke, H.G., Tschiersch, J., 1994. Statistical characteristics of the activity concentration in the surface layer of the atmosphere in the 30 km zone of Chernobyl. Journal of Aerosol Science 25, 767e777. Garger, E.K., Kashpur, V.A., Li, W.B., Tschiersch, J., 2006. Radioactive aerosols released from the Chernobyl Shelter into the immediate environment. Radiation and Environmental Biophysics 45, 105e114. Garger, E.K., Paretzke, H.G., Tschiersch, J., 1998. Measurement of resuspended aerosol in the Chernobyl area. Part III: size distribution and dry deposition velocity of radioactive particles during anthropogenic enhanced resuspension. Radiation and Environmental Biophysics 37, 201e208. Garger, E.K., Anspaugh, L.R., Shinn, J.H., Hoffman, F.O., 1995. A Test of Resuspension Factor Models against Chernobyl Date. In: Proceedings of an International Symposium on Environmental Impact of Radioactive Releases. International Atomic Energy Agency SM-339/26, Vienna, pp. 369e376. Holländer, W., Garger, E., 1996. Contamination of Surface by Resuspended Material. ECP e 1, Final Report, Rep. EUR 16527. Office for Official Publications of the European Communities, Luxembourg. Ivachnenko, A.G. (Ed.), 1980. The Handbook on Typical Software Programs for Modelling. K. Technica. Ivachnenko, A.G., Müller, J.-A., 1984. Selbstorganisation von Vorhersagemodellen. VEB Verlag Technik, Berlin. Lomb, N.R., 1976. Least-square frequency analysis of unequally spaced data. Astrophysics and Space Science 39, 447e462. Madala, N.R., Ivachnenko, A.G., 1993. Inductive Learning Algorithms for Complex System Modeling. CRC Press, Boca Raton, Florida. Paatero, J., Vesterbacka, K., Makkonen, U., Kyllönen, K., Hellen, H., Hatakka, J., Anttila, P., 2009. Resuspension of radionuclides into the atmosphere due to forest fires. Journal of Radioanalytical and Nuclear Chemistry 282, 473e476.

58

E.K. Garger et al. / Journal of Environmental Radioactivity 110 (2012) 53e58

Viswanathan, G.M., Buldyrev, S.V., Garger, E.K., Kashpur, V.A., Lucena, L.S., Shlyakhter, A., Stanley, H.E., Tschiersch, J., 2000. Quantifying nonstationary radioactivity concentration fluctuations near Chernobyl: a complete statistical description. Physical Review E 62, 4389e4392. Yoschenko, V.I., Kashparov, V.A., Protsak, V.P., Lundin, S.M., Levchuk, S.E., Kadygrib, A.M., Zvarich, S.I., Khomutinin, Yu.V., Maloshtan, I.M., Lanshin, V.P.,

Kovtun, M.V., Tschiersch, J., 2006. Resuspension and redistribution of radionuclides during grassland and forest fires in the Chernobyl exclusion zone. Part I. Fire experiments. Journal of Environmental Radioactivity 86, 143e163. Yurachkovsky, Yu.P., Popkov, N.M., 1986. Estimation of parameters in GMDH algorithms for modelling of polyharmonic processes and fields. Soviet Journal of Automation and Information Sciences 19, 11e17.