- Email: [email protected]
Data assimilation in the atmospheric dispersion model for nuclear accident assessments
ARTICLE IN PRESS
Atmospheric Environment 41 (2007) 2438–2446 www.elsevier.com/locate/atmosenv
Data assimilation in the atmospheric dispersion model for nuclear accident assessments D.Q. Zhenga, J.K.C. Leunga,, B.Y. Leeb, H.Y. Lamb a
Department of Physics, The University of Hong Kong, Hong Kong b Hong Kong Observatory, Hong Kong
Received 25 January 2006; received in revised form 26 May 2006; accepted 30 May 2006
Abstract Uncertainty factors in atmospheric dispersion models may influence the reliability of model prediction. The ability of a model in assimilating measurement data will be helpful to improve model prediction. In this paper, data assimilation based on ensemble Kalman filter (EnKF) is introduced to a Monte Carlo atmospheric dispersion model (MCADM) designed for assessment of consequences after an accident release of radionuclides. Twin experiment has been performed in which simulated ground-level dose rates have been assimilated. Uncertainties in the source term and turbulence intensity of wind field are considered, respectively. Methodologies and preliminary results of the application are described. It is shown that it is possible to reduce the discrepancy between the model forecast and the true situation by data assimilation. About 80% of error caused by the uncertainty in the source term is reduced, and the value for that caused by uncertainty in the turbulence intensity is about 50%. r 2006 Elsevier Ltd. All rights reserved. Keywords: Data assimilation; Dispersion model; Ensemble Kalman filter; Nuclear accident
1. Introduction Modeling atmospheric dispersion of radioactive plume is one of the important aspects in nuclear accident assessments to predict the radiological consequences to the public. How to give a prediction quickly and accurately is a key issue. Although many models have been developed and some of them can execute well by inputting a source term, meteorological data and dispersion parameters, in an actual nuclear accident, the real source term and dispersion parameters may never be known exactly, while the meteorological data often come from Corresponding author. Tel.: +852 28592858.
E-mail address: [email protected] (J.K.C. Leung). 1352-2310/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.atmosenv.2006.05.076
numerical weather prediction models, not from measurement. As a result, relatively large differences in the dose, for example, may be observed between results of dispersion models and those obtained from observations. A way to overcome the problem is to employ data assimilation (DA) technique. DA is a procedure that combines model and observation data to improve the prediction. This technique has been widely applied in various research fields in the last two decades, especially in the field of numerical weather prediction or meteorological pre-processing which produces meteorological data for the emergency response systems (Kovalets et al., 2003). van Loon et al. (2000) introduced DA in an atmospheric transport
ARTICLE IN PRESS D.Q. Zheng et al. / Atmospheric Environment 41 (2007) 2438–2446
chemistry model to improve modeled ozone concentration. In recent years, some studies have involved the application of DA to the field of decision support in the event of a radiological accident. In the popular real-time on-line decision support system (RODOS) DA tools are under development for several modules, including the atmospheric dispersion, deposition, food chain and hydrological modules. Rojas-Palma et al. (2003) reported some preliminary results of their work. In particular, in the atmospheric dispersion module Risø Mesoscale PUFF (RIMPUFF) model (Mikkelsen et al., 1984), they specified an extended Kalman filter as the DA algorithm. They constructed their state vector with the position and content of each puff particle. The measured dose rates near the plant were assimilated to see if the estimation at sites downstream from the plant could be improved. In our study, a Monte Carlo dispersion model is employed as the forecast model. Simulated dose rate observations are assimilated in the model. The Monte Carlo model uses hundreds of thousands of particles to represent the radioactive cloud. It has been developed as a tool to evaluate the accident consequence for the Daya Bay Nuclear Power Plant in south China. DA is based on the use of the ensemble Kalman filter (EnKF) (Evensen et al., 1997). The release rate of nuclides (one important parameter in source term) and the turbulence of wind field are taken as uncertain, because they are estimated on the basis of some simple assumptions and empirical formulae. No other model errors are taken into account yet. The source term is uncertain because the model uses estimated release rate as the real release rate, and the release rate remains constant in the model, so even though the total emissions are correct, the temporal distribution may be incorrect. The turbulence of wind field comes from empirical parameterization, and their uncertainty can result in significant error in the model. In the next section, EnKF and the methodology of introducing DA into the Monte Carlo dispersion model will be described. In Section 3, some preliminary results of our work will be given. 2. DA methodology This section gives a brief introduction to the methodology of the study, including EnKF, Monte Carlo dispersion model, and the DA procedure. The construction of our DA model is also described.
2439
2.1. Ensemble Kalman filter The standard Kalman filter is a linear, unbiased, and minimum-error variance recursive algorithm to optimally estimate the unknown state of a linear dynamic system from noisy data taken at discrete real-time intervals. There are two stages involved in the Kalman filter: prediction and update. In the prediction stage, one solves the full state evolution and advances the error-covariance matrix in time. The update or analysis stage merges the model prediction and current observations, to provide an analyzed or assimilated state. The Kalman filter can be described by the following equations: xk ¼ Mxk1 þ Z,
(1)
Pfk ¼ MPk1 M T þ Q,
(2)
yo ¼ Hx þ e,
(3)
K ¼ Pf H T ðHPf H T þ RÞ1 ,
(4)
xa ¼ xf þ Kðyo Hxf Þ,
(5)
Pa ¼ ðI KHÞPf .
(6)
Here xk denotes the state vector, and the subscript k or k 1 represents the time step. The superscript f indicates forecast, the superscript o indicates observation, and a means analyzed. The observation vector is y, which contains the measured data. M is the forecast model operator. H is the observation matrix which transforms the state space to observation space. Z is the random noise with a covariance matrix Q, and e is the measurement error vector with a covariance matrix R. P is covariance of the prediction error. K is called the gain matrix. I is the unit matrix. Eqs. (1)–(3) represent the forecast stage, which predicts the state vector of the next time step following the forecast model. Eq. (1) describes the model as a stochastic model instead of a deterministic model which can be described as xk ¼ Mxk1 . The random noise vector Z represents the uncertainty in the model. The extension of a deterministic model to a stochastic model is the basic of Kalman filter. Eqs. (4) and (5) represent the scheme of analysis, which is also an update for the forecast state vector. The gain matrix K determines the extent of the update. The result is an analyzed state vector, which, according to the theory, represents the
ARTICLE IN PRESS 2440
D.Q. Zheng et al. / Atmospheric Environment 41 (2007) 2438–2446
optimal estimation of the state vector (Grewal and Andrews, 1993). The elements in the state vector x can be any variables and parameters we are concerned about in the model. The state vector x is updated over time in two ways: (i) it is allowed to evolve according to the model dynamics; and (ii) it is modified periodically to take into account new observations. The observation vector may or may not include the same variables as the state vector. In our case, the observations are the ground-level dose rates. In the absence of aerial monitoring capability, measuring the nuclide concentration aloft in the atmosphere is very difficult during a nuclear accident. Instead, most measurements are gamma dose rates coming from the ground-level observation stations. The relationship between the dose rate and the concentration will be described in Section 2.3.3. The standard Kalman filter requires a linear system, which includes linear dynamics (Eq. (1)) and a linear measurement equation (Eq. (3)). In a nonlinear system, we are unable to obtain the matrix M and H. An extended Kalman filter partly solves the problem, but it requires the more involved tangent linear and adjoint codes (Ide and Ghil, 1997). Currently very popular, the EnKF proposed by Evensen (1994) takes a statistical approach to the solution of the Kalman filter equations for the covariance matrices of analysis and forecast error. It applies an ensemble of model states to obtain the error statistics of the model estimate, and carries out ensemble integrations to predict the error statistics forward in time. An EnKF uses a random generator to produce noise inputs in order to form an ensemble of state vectors. Each member of the ensemble is a point in the state space. The ensemble of model states is integrated forward in time according to the stochastic model (Eq. (1)). In a practical implementation this becomes just a standard integration of the numerical model but subject to a stochastic noise which resembles the uncertainties in the model. With this ensemble of states, the sample covariance matrices can be obtained by using the following equation: Pf ¼
N 1 X ðxf xf Þðxfi xf ÞT , N 1 i¼1 i
(7)
where the index i refers to sample (ensemble) member, and the x is the sample mean value.
In EnKF, the sample covariance matrices are used as the covariance matrices of the forecast error. When the ensemble size is big enough, we can control the estimate error to within a small value. There are several approaches to handle the analysis scheme in EnKF. In this paper, we followed the approach given by Evensen (2003). The analysis scheme operates directly on the ensemble of model states when observations are assimilated. The analysis ensemble is obtained by assimilating a different set of observations to each member of the state ensemble. The different sets of observations are created by adding random noise to the real observations. The random noise component being generated according to the observational error-covariance matrix. 2.2. Monte Carlo dispersion model The DA process is built on the basis of the Monte Carlo dispersion model (MCADM) which forecasts the nuclide concentration by following the path of hundreds of thousands of particles. In this model, the 3D space is divided into grids. The horizon grid size is 5 km 5 km. In the vertical, there are 13 layers within the 2 km above ground. Each layer has different depth, the height of the layers being 2; 10; 25; 50; 75; 100; 200; 400; 600; 800; 1000; 1500 and 2000 m above ground. For each layer, there are 62 44 grids. At each time step, the movement of each particle is traced and the concentration in each grid is obtained by counting the particles within the grid. This model was designed to simulate the dispersion of radionuclides after an accidental release around Daya Bay nuclear power plants, and has been evaluated by an SF6 field tracing experiment (Zhang et al., 1999). The input data of the model include three parts: source term, terrain data and meteorological data. Source term consists of total release radio activities of the accident, release time of duration, height of release, etc. Terrain data include the height above sea level and surface roughness length of each grid point. The meteorological input consists of 1-hourly 3D fields with the wind, temperature and the 2D field of rainfall. The MCADM releases 2000 particles from the source point per time step and trace the trajectories of the each particles. The location of each individual particle is defined by xj ðt þ DtÞ ¼ xj ðtÞ þ ½uj ðtÞ þ u0j ðtÞDt;
j ¼ 1; 2; . . . . (8)
ARTICLE IN PRESS D.Q. Zheng et al. / Atmospheric Environment 41 (2007) 2438–2446
The variable uj is the mean value of any component of the particle velocity, u0j represents the turbulent fluctuation, subscript j represents each particle. uj is obtained by numerical interpolation of the wind velocity components from the input data. u0j is the random term with mean value zero and standard deviation su . It is generated by a random generator. su is the standard deviation of turbulence wind velocity component. It is a parameter related with turbulence intensity. Usually su cannot be obtained directly from the input data, because it is not the routine observation data of the weather stations. PBL parameterizations were used to obtain it by applying semiempirical formula (Zhang et al., 2000; Zannettip, 1990). Uncertainties in the dispersion model mainly arise from: (1) Uncertainties or errors in model input data: source term, wind field or temperature field. (2) Internal model errors or inadequate model physics. In the MCADM, turbulence parameters estimation is one of the limitations. It is hard to say that PBL parameterizations can give the accurate estimation by using semiempirical formula. This paper has done some attempt to study if the two types of uncertainties can be reduced by introducing DA. For the first type, as has been said, we selected the release rate of radioactivity as the representation; for the second, turbulence intensity is most reasonable to be selected as uncertain. Here, we select the standard deviation of turbulence wind velocity su as uncertain. In this paper, the uncertainty in the release rate and turbulence is considered, respectively. 2.3. Set-up of assimilation model In order to apply EnKF to improve the forecast of the Monte Carlo dispersion model, we have to build the similar construct as the Kalman filter algorithm. It is necessary to first define the state vector x and the observation vector y, and calculate the H matrix (see Eq. (3)). 2.3.1. Selection of state vector As a normal practice, we select the observation vector y to include the gamma dose rate measurements at all observation points. However, there are several choices to construct the state vector. In the Monte Carlo model, the quantities evolving with time are the particle positions, so a direct method of selecting the state vector is to include the positions of all particles. This has not been done because (a)
2441
the gamma dose rate at an observation point is a nonlinear function of the positions of these particles; (b) the computation time used by the analysis procedure involving a large dimension of the state vector is cost prohibitive, because a few hundreds of thousands of particles are used to simulate one accident release. In our case, the state vector is selected to include the activity concentrations of the nuclide at each grid point. There are a total of 36,608 grids in the model, but at any given time, only a few thousands of grids are affected by the plume. To optimize the computation, only the grid points where the concentration is not zero are included in the state vector: x ¼ ðC i Þ
for all C i 40,
(9)
where C i is the concentration of the ith grid. In order to do that, a dynamic selection has been added to the program, and the elements of the state vector have to be selected dynamically whenever the measurement data are assimilated. With this, the computation time has been reduced significantly, which is very important during an emergency situation. The computer program is run on the HPCPOWER cluster in the University of Hong Kong. If eight nodes are used, the program can produce a 24 h forecast in about 4 h. The total computation time depends on the extent of the plume. 2.3.2. Specification of model uncertainty For DA the specification of the model uncertainty is required. In EnKF, random noise can be added to any poorly known model parameters when the ensemble of states is integrated forward in time according to the model equation. In standard Kalman filter, the errors are assumed to be Gaussian white noise. In practice, model errors usually cannot satisfy this condition, even if the underlying uncertainties are Gaussian. For example, if an Gaussian distributed error is added to su , the error of the predicted concentration will not be Gaussian distributed. The limitation for the white noise is also hard to be satisfied. For the MCADM model, it will be more reasonable to assume the error in the release rates as colored noise instead of white noise. Although the release rates are uncertain, they will not fluctuate very fast. The following equation can be used to simulate time correlated model errors: pffiffiffiffiffiffiffiffiffiffiffiffiffi ck ¼ ack1 þ 1 a2 wk1 . (10)
ARTICLE IN PRESS D.Q. Zheng et al. / Atmospheric Environment 41 (2007) 2438–2446
2442
Here ck is a sequence of time correlated errors. wk is a white noise with mean value zero. a determines the time decorrelation between the sequence of the random errors. The covariance in time between ck and ck1 is ck ck1 ¼ a ¼ expðDt=tÞ,
(11)
where Dt is the time step, and t is the time decorrelation length. EnKF allows for a wide range of noise models. It is not restricted to use Gaussian distributed noise, and it is possible to use time correlated noise by transforming it into white noise by introducing a new state vector which consists of x augmented with c (Evensen, 2003). Model bias can also be regarded as one kind of colored noises. In fact, the distinction between model bias and colored noises is not clear. In Eq. (10), if a ¼ 1, ck will be a error which is constant in time. One could also consider the error to be a model bias. 2.3.3. Observation matrix calculation The elements of the observation matrix H (see Eq. (3)) are the derivatives of the measurement prediction with respect to the state vector variables. In our case, the element of the matrix H, hij , represents the contribution of unit concentration at grid i to the gamma dose rate at the jth observation point. Suppose the activity concentration in the ith grid is C i . Its contribution to the gamma effective dose rate of energy E at the jth observation point is (Leung, 1992) ZZZ Y ðEÞCðx; y; zÞBðE; mrij Þ d ij ðEÞ ¼ 4pr2 Vi m=r expðmrÞE dx dy dz,
ð12Þ
where V i is the volume of the grid i, C is the concentration in Bq m3 of the nuclide being considered, Y is the yield to the specified energy, B is the build-up factor, and m is the attenuation coefficient. Refer to Shleien et al. (1998) for further details of these factors. Also, r is the distance between the source point ðx; y; zÞ and the observation point, and mr represents the attenuation length. m=r is the mass energy absorption coefficient for tissue ðcm2 g1 Þ. The total dose rate is obtained by summing over all energy groups. 2.3.4. Parallelization of EnKF The computation power required in the DA precess usually exceeds the available computer
resources. One effective way to solve the problem is to use parallel programs on fast multiprocessor machines. In the context of EnKF (Evensen, 1994), it provided a convenient way to propagate the members of the ensemble independently on different processors. This makes it possible to apply Kalman filter to a model leaving the model intact. Keppenne (2000) applied this method to a large two-layer shallow water model. The main issue is in the analysis step, where the data from all ensemble members are required even though only one measurement is to be assimilated. The time spent on transfer between different processors will be increased. Usually there are two choices in the analysis step: member-decomposition (Keppenne, 2000) and domain-decomposition (Houtekamer and Mitchell, 1998). The latter is usually restricted to domain-decomposition parallelized model. In our case, the member-decomposition is used. Another point that should be noted is the random number generation process. Paralleled generator should be used instead of a serial one in order to avoid the strong correlation between random numbers generated by different processors. 3. Results Twin experiments were used to evaluate the performance of the EnKF as the DA methodology for the Monte Carlo dispersion model. Twin experiment is a sophisticated method for testing a DA algorithm, especially when there are no real observations. Even though there are some observation data, it is worthwhile to test a DA method by using simulated measured data before assimilating real ones, because in twin experiments, the model parameter error and measurement error can easily be controlled. Twin experiments include two parts: ‘truth’ run and simulation run. The ‘truth’ run uses the default meteorological parameter values and source terms, and provides a reference for the ‘true’ situation. Another task of the ‘truth’ run is to produce simulated observations by adding a random measurement error to the ‘true’ value. These simulated observations are to be assimilated in the simulation run. The simulation run is similar to the real assimilation process, except that simulated measurements, not real measurements, are to be assimilated. It differs from the ‘truth’ run in two ways—one is that the wind turbulence or source term is perturbed, which represents the error in the
ARTICLE IN PRESS D.Q. Zheng et al. / Atmospheric Environment 41 (2007) 2438–2446
estimated parameters, and the other is that the simulated measured data are assimilated in the model. Because in the real simulation procedure, the real source term and wind turbulence are not known, what we know is only an estimated value which includes errors in it. After the simulation run, the results can be compared with those results obtained by using the disturbed source term and meteorological parameters, but without DA, to see if the prediction is improved. In this paper, the forecast covers a period of over 24 h. The simulated measurements are assimilated every hour. We used the current observation network in south China and Hong Kong which is shown in Fig. 1. Simulated measured dose rate at the observation stations will be assimilated to see the impact of DA on the results.
2443
The simulated true source term and the source term used in the forecast model (in which a constant release rate is used) is shown in Fig. 2, in which the solid line represents the simulated true source term, while the dotted straight line represents the source term used in the deterministic forecast model. For assimilation, the simulated ‘truth’ is generated with a model using a randomly disturbed source term. The ground-level dose rate for all the observation points are also calculated with the ‘truth’ run, and the measurement are simulated from the ‘truth’ dose rate disturbed with a 5% measurement error. The measured data are assimilated every hour.
3.1. Uncertainty in source term In this experiment, only the nuclide release rate is considered uncertain. With an average release rate r in the deterministic model, the stochastic release rate was modeled according to rk ¼ maxð0; rð1 þ ck ÞÞ,
(13)
where ck is the colored noise which can be produced according to Eq. (10). The standard deviation of the error is 100% of r, and the time decorrelation length t is 2 h.
Release rate (Bq/s)
2.5×106
Forcecast source term True source term
2.0×106 1.5×106 1.0×106 5.0×105 0.0 0
2
4
6
8
10
Time (h) Fig. 2. Simulated true source term and the forecast source term.
Fig. 1. Observation stations in the model region. The star represents the location of power plant.
ARTICLE IN PRESS D.Q. Zheng et al. / Atmospheric Environment 41 (2007) 2438–2446
2444
2.0×10-6
3.5×10-6
True Without data assimilation With data assimilation
Doserate (Gy/s)
Doserate (Gy/s)
3.0×10-6 2.5×10-6 2.0×10-6 1.5×10-6 1.0×10-6
True Without data assimilation With data assimilation
1.5×10-6
1.0×10-6
5.0×10-7
5.0×10-7 0.0
0.0 0
(a)
5
10
15
20
Time(h)
0
5
10
15
20
Time (h)
(b)
Fig. 3. Dose rates from the true source term and the forecast source term (with and without assimilation) for location (a) 101 and (b) 602, respectively.
3.2. Uncertainty in the turbulence intensity of wind field Similar to the test that uses randomly perturbed source term, the perturbed turbulence parameter, here, the standard deviation of the turbulence velocity su , is used this time. As the error in the
6×106 Without assimilation With assimilation
5×106 RMS (Bq/m3)
Fig. 3(a) and (b) shows some assimilation experiment results with a randomly perturbed source term for two locations which are labeled with 101 and 602, respectively, in Fig. 1. The station 101 is at Da Peng town which is about 8 km west of the power plant and is the biggest town near Daya Bay plant. The station 602 is at Kat O island in Hong Kong which is about 25 km away from the power plant. In the deterministic forecast model, the release rate for the entire release period is constant, so the trend of dose rate coming out from the forecast model without assimilation is relatively flat. The forecast source term is corrected by assimilating the measurements which gives a dose rate closer to the true value than that without assimilation. The result in respect of the point 101 follows very closely the true value. In comparison, for the point 602, the estimated value is a little higher than the true value during times of relatively high dose rate, but the trend is in general similar to the true value. In order to see the overall impact of DA to the model prediction, the root mean square (RMS) of the difference between the estimated with the ‘true’ concentration is shown in Fig. 4. From the figure we can see that about 80% of the error that arise from the uncertain in the release rate is overcome. This result is encouraging.
4×106 3×106 2×106 1×106 0 0
5
10
15
20
25
Time (h) Fig. 4. RMS of the estimated vs. ‘true’ concentration with and without DA when release rate is taken as uncertain.
turbulence intensity influences the diffusion parameters directly, it is also an important parameter in the dispersion model. On the other side, because the parameter is estimated by PBL parameterizations, it is prone to involve errors in it. PBL parameterization is an simplification of the model physics. It may overestimate the turbulence intensity under a certain situation and underestimate it under another one. So, in this experiment, it is assumed that there is a bias in the turbulence estimation. In the determinate model, su is perturbed by increasing 30% of the ‘true’ value. The purpose of the experiment is to see if the DA method can work well in case that the forecast turbulence intensity is 30% higher than the true value. The comparison between the concentration distribution with and without DA is shown in Fig. 5. The forecast concentration without DA is found to be lower than the true value. This is
ARTICLE IN PRESS D.Q. Zheng et al. / Atmospheric Environment 41 (2007) 2438–2446
because the estimated turbulence strength is higher than the ‘true’ value. After DA, the difference from the ‘true’ concentration distribution is found to be
2445
1.2×106 Without assimilation With assimilstion
RMS (Bq/s3)
1.0×106 8.0×105 6.0×105 4.0×105 2.0×105 0.0 0
5
10
15 Time (h)
20
25
Fig. 6. RMS of the estimated vs. ’true’ concentration with and without DA when su is taken as uncertain.
reduced, implying that the assimilated measurements have updated the forecast concentration. The RMS of the difference between the estimated state and the ‘true’ state in this case is shown in Fig. 6. It shows that the error is reduced after DA, and the degree is about 50%.
4. Conclusion The results of the twin experiments show that successful application of the ensemble Kalman filter in combination of a nonlinear Monte Carlo dispersion model is possible. Dynamic selection of the state vector seems to be a feasible method to improve the efficiency of assimilation in case of accident. Introducing colored noise instead of white noise is more reasonable, and EnKF can work well under this situation. Although EnKF provides a possibility of combining model and measurement together, much work must be done before the methodology can be applied in a real accident emergency system. Here are some recommendations:
Fig. 5. Contours of radionuclide concentration 10 h after release with and without assimilation. (a) The true situation. (b) The forecast alone with perturbed turbulence. (c) The estimation after the measured dose rate is assimilated.
(i) Other uncertainties except source term and turbulence intensity (i.e. uncertainties in wind direction, source height estimations, etc.) should be taken into account, and the attempt with various uncertainties acting together should be done. (ii) Twin experiment only provides a way for validating the methodology in a simulation way. It is desirable to validate it against realcase measurements.
ARTICLE IN PRESS 2446
D.Q. Zheng et al. / Atmospheric Environment 41 (2007) 2438–2446
(iii) Considering that the real nuclear accident measurement data are difficult to be obtained, wind tunnel experiments can be used to get some real measurements. References Evensen, G., 1994. Sequential data assimilation with a non-linear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. Journal of Geophysical Research 99, 10,143–10,162. Evensen, G., 2003. The ensemble Kalman filter: theoretical formulation and practical implementation. Ocean Dynamics 53, 343–367. Evensen, G., Madsen, H., Gering, F., 1997. Advanced data assimilation for strongly non-linear dynamics. Monthly Weather Review 25, 1342–1354. Grewal, M.S., Andrews, A.P., 1993. Kalman Filtering: Theory and Practice. Prentice-Hall, Englewood Cliffs, NJ. Houtekamer, P.L., Mitchell, H.L., 1998. Data assimilation using an ensemble Kalman filter technique. Monthly Weather Review 126, 796–811. Ide, K., Ghil, M., 1997. Extended Kalman filtering for vortex systems. Part I: methodology and point vortices. Dynamics of Atmospheres and Oceans 27, 301–332. Keppenne, C.L., 2000. Data assimilation into a primitiveequation model with a parallel ensemble Kalman filter. Monthly Weather Review 128, 1971–1981. Kovalets, I., Andronopoulos, S., Bartzis, J.G., Gounaris, N., Kushchan, A., 2003. Introduction of data assimilation
procedures in the meteorological pre-processor of atmospheric dispersion models used in emergency response systems. Atmospheric Environment 38, 457–467. Leung, J.K.C., 1992. Application of shielding factors for protection against gamma radiations during a nuclear accident. IEEE Transactions on Nuclear Science 39, 1512–1518. Mikkelsen, T., Larsen, S.E., Thykier-Nielsen, S., 1984. Description of the RisøPuff diffusion model. Nuclear Technology 67, 56–65. Rojas-Palma, C., Madsen, H., Cering, F., Puch, R., Turcanu, C., Astrup, P., Mu¨ller, H., Richter, K., Zheleznyak, M., Treebushny, D., Kolomeev, M., Kamaev, D., Wynn, H., 2003. Data assimilation in the decision support system RODOS. Radiation Protection Dosimetry 104, 31–40. Shleien, B., Slaback Jr., L.A., Birky, B.K., 1998. Handbook of Health Physics and Radiological Health. Williams & Wilkins, Baltimore, MD. van Loon, M., Builtjes, P.J.H., Segers, A.J., 2000. Data assimilation of ozone in the atmospheric transport chemistry model LOTOS. Environmental Modeling and Software 15, 603–609. Zannettip, P., 1990. Air pollution modelling: theories, computational method and available software. Van Nostrand Reinhold, New York. Zhang, M.G., Han, Z.W., Lei, X.E., 1999. Application of Monte Carlo model for multiple sources to air quality assessment at atomic power plant of Guangdong. Climatic and Environmental Research 14, 203–209. Zhang, M.G., Han, Z.W., Lei, X.E., 2000. A numerical study of atmospheric dispersion in the area of nuclear power plant Dayawan. Climatic and Environmental Research 5, 90–95.
Atmospheric Environment 41 (2007) 2438–2446 www.elsevier.com/locate/atmosenv
Data assimilation in the atmospheric dispersion model for nuclear accident assessments D.Q. Zhenga, J.K.C. Leunga,, B.Y. Leeb, H.Y. Lamb a
Department of Physics, The University of Hong Kong, Hong Kong b Hong Kong Observatory, Hong Kong
Received 25 January 2006; received in revised form 26 May 2006; accepted 30 May 2006
Abstract Uncertainty factors in atmospheric dispersion models may influence the reliability of model prediction. The ability of a model in assimilating measurement data will be helpful to improve model prediction. In this paper, data assimilation based on ensemble Kalman filter (EnKF) is introduced to a Monte Carlo atmospheric dispersion model (MCADM) designed for assessment of consequences after an accident release of radionuclides. Twin experiment has been performed in which simulated ground-level dose rates have been assimilated. Uncertainties in the source term and turbulence intensity of wind field are considered, respectively. Methodologies and preliminary results of the application are described. It is shown that it is possible to reduce the discrepancy between the model forecast and the true situation by data assimilation. About 80% of error caused by the uncertainty in the source term is reduced, and the value for that caused by uncertainty in the turbulence intensity is about 50%. r 2006 Elsevier Ltd. All rights reserved. Keywords: Data assimilation; Dispersion model; Ensemble Kalman filter; Nuclear accident
1. Introduction Modeling atmospheric dispersion of radioactive plume is one of the important aspects in nuclear accident assessments to predict the radiological consequences to the public. How to give a prediction quickly and accurately is a key issue. Although many models have been developed and some of them can execute well by inputting a source term, meteorological data and dispersion parameters, in an actual nuclear accident, the real source term and dispersion parameters may never be known exactly, while the meteorological data often come from Corresponding author. Tel.: +852 28592858.
E-mail address: [email protected] (J.K.C. Leung). 1352-2310/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.atmosenv.2006.05.076
numerical weather prediction models, not from measurement. As a result, relatively large differences in the dose, for example, may be observed between results of dispersion models and those obtained from observations. A way to overcome the problem is to employ data assimilation (DA) technique. DA is a procedure that combines model and observation data to improve the prediction. This technique has been widely applied in various research fields in the last two decades, especially in the field of numerical weather prediction or meteorological pre-processing which produces meteorological data for the emergency response systems (Kovalets et al., 2003). van Loon et al. (2000) introduced DA in an atmospheric transport
ARTICLE IN PRESS D.Q. Zheng et al. / Atmospheric Environment 41 (2007) 2438–2446
chemistry model to improve modeled ozone concentration. In recent years, some studies have involved the application of DA to the field of decision support in the event of a radiological accident. In the popular real-time on-line decision support system (RODOS) DA tools are under development for several modules, including the atmospheric dispersion, deposition, food chain and hydrological modules. Rojas-Palma et al. (2003) reported some preliminary results of their work. In particular, in the atmospheric dispersion module Risø Mesoscale PUFF (RIMPUFF) model (Mikkelsen et al., 1984), they specified an extended Kalman filter as the DA algorithm. They constructed their state vector with the position and content of each puff particle. The measured dose rates near the plant were assimilated to see if the estimation at sites downstream from the plant could be improved. In our study, a Monte Carlo dispersion model is employed as the forecast model. Simulated dose rate observations are assimilated in the model. The Monte Carlo model uses hundreds of thousands of particles to represent the radioactive cloud. It has been developed as a tool to evaluate the accident consequence for the Daya Bay Nuclear Power Plant in south China. DA is based on the use of the ensemble Kalman filter (EnKF) (Evensen et al., 1997). The release rate of nuclides (one important parameter in source term) and the turbulence of wind field are taken as uncertain, because they are estimated on the basis of some simple assumptions and empirical formulae. No other model errors are taken into account yet. The source term is uncertain because the model uses estimated release rate as the real release rate, and the release rate remains constant in the model, so even though the total emissions are correct, the temporal distribution may be incorrect. The turbulence of wind field comes from empirical parameterization, and their uncertainty can result in significant error in the model. In the next section, EnKF and the methodology of introducing DA into the Monte Carlo dispersion model will be described. In Section 3, some preliminary results of our work will be given. 2. DA methodology This section gives a brief introduction to the methodology of the study, including EnKF, Monte Carlo dispersion model, and the DA procedure. The construction of our DA model is also described.
2439
2.1. Ensemble Kalman filter The standard Kalman filter is a linear, unbiased, and minimum-error variance recursive algorithm to optimally estimate the unknown state of a linear dynamic system from noisy data taken at discrete real-time intervals. There are two stages involved in the Kalman filter: prediction and update. In the prediction stage, one solves the full state evolution and advances the error-covariance matrix in time. The update or analysis stage merges the model prediction and current observations, to provide an analyzed or assimilated state. The Kalman filter can be described by the following equations: xk ¼ Mxk1 þ Z,
(1)
Pfk ¼ MPk1 M T þ Q,
(2)
yo ¼ Hx þ e,
(3)
K ¼ Pf H T ðHPf H T þ RÞ1 ,
(4)
xa ¼ xf þ Kðyo Hxf Þ,
(5)
Pa ¼ ðI KHÞPf .
(6)
Here xk denotes the state vector, and the subscript k or k 1 represents the time step. The superscript f indicates forecast, the superscript o indicates observation, and a means analyzed. The observation vector is y, which contains the measured data. M is the forecast model operator. H is the observation matrix which transforms the state space to observation space. Z is the random noise with a covariance matrix Q, and e is the measurement error vector with a covariance matrix R. P is covariance of the prediction error. K is called the gain matrix. I is the unit matrix. Eqs. (1)–(3) represent the forecast stage, which predicts the state vector of the next time step following the forecast model. Eq. (1) describes the model as a stochastic model instead of a deterministic model which can be described as xk ¼ Mxk1 . The random noise vector Z represents the uncertainty in the model. The extension of a deterministic model to a stochastic model is the basic of Kalman filter. Eqs. (4) and (5) represent the scheme of analysis, which is also an update for the forecast state vector. The gain matrix K determines the extent of the update. The result is an analyzed state vector, which, according to the theory, represents the
ARTICLE IN PRESS 2440
D.Q. Zheng et al. / Atmospheric Environment 41 (2007) 2438–2446
optimal estimation of the state vector (Grewal and Andrews, 1993). The elements in the state vector x can be any variables and parameters we are concerned about in the model. The state vector x is updated over time in two ways: (i) it is allowed to evolve according to the model dynamics; and (ii) it is modified periodically to take into account new observations. The observation vector may or may not include the same variables as the state vector. In our case, the observations are the ground-level dose rates. In the absence of aerial monitoring capability, measuring the nuclide concentration aloft in the atmosphere is very difficult during a nuclear accident. Instead, most measurements are gamma dose rates coming from the ground-level observation stations. The relationship between the dose rate and the concentration will be described in Section 2.3.3. The standard Kalman filter requires a linear system, which includes linear dynamics (Eq. (1)) and a linear measurement equation (Eq. (3)). In a nonlinear system, we are unable to obtain the matrix M and H. An extended Kalman filter partly solves the problem, but it requires the more involved tangent linear and adjoint codes (Ide and Ghil, 1997). Currently very popular, the EnKF proposed by Evensen (1994) takes a statistical approach to the solution of the Kalman filter equations for the covariance matrices of analysis and forecast error. It applies an ensemble of model states to obtain the error statistics of the model estimate, and carries out ensemble integrations to predict the error statistics forward in time. An EnKF uses a random generator to produce noise inputs in order to form an ensemble of state vectors. Each member of the ensemble is a point in the state space. The ensemble of model states is integrated forward in time according to the stochastic model (Eq. (1)). In a practical implementation this becomes just a standard integration of the numerical model but subject to a stochastic noise which resembles the uncertainties in the model. With this ensemble of states, the sample covariance matrices can be obtained by using the following equation: Pf ¼
N 1 X ðxf xf Þðxfi xf ÞT , N 1 i¼1 i
(7)
where the index i refers to sample (ensemble) member, and the x is the sample mean value.
In EnKF, the sample covariance matrices are used as the covariance matrices of the forecast error. When the ensemble size is big enough, we can control the estimate error to within a small value. There are several approaches to handle the analysis scheme in EnKF. In this paper, we followed the approach given by Evensen (2003). The analysis scheme operates directly on the ensemble of model states when observations are assimilated. The analysis ensemble is obtained by assimilating a different set of observations to each member of the state ensemble. The different sets of observations are created by adding random noise to the real observations. The random noise component being generated according to the observational error-covariance matrix. 2.2. Monte Carlo dispersion model The DA process is built on the basis of the Monte Carlo dispersion model (MCADM) which forecasts the nuclide concentration by following the path of hundreds of thousands of particles. In this model, the 3D space is divided into grids. The horizon grid size is 5 km 5 km. In the vertical, there are 13 layers within the 2 km above ground. Each layer has different depth, the height of the layers being 2; 10; 25; 50; 75; 100; 200; 400; 600; 800; 1000; 1500 and 2000 m above ground. For each layer, there are 62 44 grids. At each time step, the movement of each particle is traced and the concentration in each grid is obtained by counting the particles within the grid. This model was designed to simulate the dispersion of radionuclides after an accidental release around Daya Bay nuclear power plants, and has been evaluated by an SF6 field tracing experiment (Zhang et al., 1999). The input data of the model include three parts: source term, terrain data and meteorological data. Source term consists of total release radio activities of the accident, release time of duration, height of release, etc. Terrain data include the height above sea level and surface roughness length of each grid point. The meteorological input consists of 1-hourly 3D fields with the wind, temperature and the 2D field of rainfall. The MCADM releases 2000 particles from the source point per time step and trace the trajectories of the each particles. The location of each individual particle is defined by xj ðt þ DtÞ ¼ xj ðtÞ þ ½uj ðtÞ þ u0j ðtÞDt;
j ¼ 1; 2; . . . . (8)
ARTICLE IN PRESS D.Q. Zheng et al. / Atmospheric Environment 41 (2007) 2438–2446
The variable uj is the mean value of any component of the particle velocity, u0j represents the turbulent fluctuation, subscript j represents each particle. uj is obtained by numerical interpolation of the wind velocity components from the input data. u0j is the random term with mean value zero and standard deviation su . It is generated by a random generator. su is the standard deviation of turbulence wind velocity component. It is a parameter related with turbulence intensity. Usually su cannot be obtained directly from the input data, because it is not the routine observation data of the weather stations. PBL parameterizations were used to obtain it by applying semiempirical formula (Zhang et al., 2000; Zannettip, 1990). Uncertainties in the dispersion model mainly arise from: (1) Uncertainties or errors in model input data: source term, wind field or temperature field. (2) Internal model errors or inadequate model physics. In the MCADM, turbulence parameters estimation is one of the limitations. It is hard to say that PBL parameterizations can give the accurate estimation by using semiempirical formula. This paper has done some attempt to study if the two types of uncertainties can be reduced by introducing DA. For the first type, as has been said, we selected the release rate of radioactivity as the representation; for the second, turbulence intensity is most reasonable to be selected as uncertain. Here, we select the standard deviation of turbulence wind velocity su as uncertain. In this paper, the uncertainty in the release rate and turbulence is considered, respectively. 2.3. Set-up of assimilation model In order to apply EnKF to improve the forecast of the Monte Carlo dispersion model, we have to build the similar construct as the Kalman filter algorithm. It is necessary to first define the state vector x and the observation vector y, and calculate the H matrix (see Eq. (3)). 2.3.1. Selection of state vector As a normal practice, we select the observation vector y to include the gamma dose rate measurements at all observation points. However, there are several choices to construct the state vector. In the Monte Carlo model, the quantities evolving with time are the particle positions, so a direct method of selecting the state vector is to include the positions of all particles. This has not been done because (a)
2441
the gamma dose rate at an observation point is a nonlinear function of the positions of these particles; (b) the computation time used by the analysis procedure involving a large dimension of the state vector is cost prohibitive, because a few hundreds of thousands of particles are used to simulate one accident release. In our case, the state vector is selected to include the activity concentrations of the nuclide at each grid point. There are a total of 36,608 grids in the model, but at any given time, only a few thousands of grids are affected by the plume. To optimize the computation, only the grid points where the concentration is not zero are included in the state vector: x ¼ ðC i Þ
for all C i 40,
(9)
where C i is the concentration of the ith grid. In order to do that, a dynamic selection has been added to the program, and the elements of the state vector have to be selected dynamically whenever the measurement data are assimilated. With this, the computation time has been reduced significantly, which is very important during an emergency situation. The computer program is run on the HPCPOWER cluster in the University of Hong Kong. If eight nodes are used, the program can produce a 24 h forecast in about 4 h. The total computation time depends on the extent of the plume. 2.3.2. Specification of model uncertainty For DA the specification of the model uncertainty is required. In EnKF, random noise can be added to any poorly known model parameters when the ensemble of states is integrated forward in time according to the model equation. In standard Kalman filter, the errors are assumed to be Gaussian white noise. In practice, model errors usually cannot satisfy this condition, even if the underlying uncertainties are Gaussian. For example, if an Gaussian distributed error is added to su , the error of the predicted concentration will not be Gaussian distributed. The limitation for the white noise is also hard to be satisfied. For the MCADM model, it will be more reasonable to assume the error in the release rates as colored noise instead of white noise. Although the release rates are uncertain, they will not fluctuate very fast. The following equation can be used to simulate time correlated model errors: pffiffiffiffiffiffiffiffiffiffiffiffiffi ck ¼ ack1 þ 1 a2 wk1 . (10)
ARTICLE IN PRESS D.Q. Zheng et al. / Atmospheric Environment 41 (2007) 2438–2446
2442
Here ck is a sequence of time correlated errors. wk is a white noise with mean value zero. a determines the time decorrelation between the sequence of the random errors. The covariance in time between ck and ck1 is ck ck1 ¼ a ¼ expðDt=tÞ,
(11)
where Dt is the time step, and t is the time decorrelation length. EnKF allows for a wide range of noise models. It is not restricted to use Gaussian distributed noise, and it is possible to use time correlated noise by transforming it into white noise by introducing a new state vector which consists of x augmented with c (Evensen, 2003). Model bias can also be regarded as one kind of colored noises. In fact, the distinction between model bias and colored noises is not clear. In Eq. (10), if a ¼ 1, ck will be a error which is constant in time. One could also consider the error to be a model bias. 2.3.3. Observation matrix calculation The elements of the observation matrix H (see Eq. (3)) are the derivatives of the measurement prediction with respect to the state vector variables. In our case, the element of the matrix H, hij , represents the contribution of unit concentration at grid i to the gamma dose rate at the jth observation point. Suppose the activity concentration in the ith grid is C i . Its contribution to the gamma effective dose rate of energy E at the jth observation point is (Leung, 1992) ZZZ Y ðEÞCðx; y; zÞBðE; mrij Þ d ij ðEÞ ¼ 4pr2 Vi m=r expðmrÞE dx dy dz,
ð12Þ
where V i is the volume of the grid i, C is the concentration in Bq m3 of the nuclide being considered, Y is the yield to the specified energy, B is the build-up factor, and m is the attenuation coefficient. Refer to Shleien et al. (1998) for further details of these factors. Also, r is the distance between the source point ðx; y; zÞ and the observation point, and mr represents the attenuation length. m=r is the mass energy absorption coefficient for tissue ðcm2 g1 Þ. The total dose rate is obtained by summing over all energy groups. 2.3.4. Parallelization of EnKF The computation power required in the DA precess usually exceeds the available computer
resources. One effective way to solve the problem is to use parallel programs on fast multiprocessor machines. In the context of EnKF (Evensen, 1994), it provided a convenient way to propagate the members of the ensemble independently on different processors. This makes it possible to apply Kalman filter to a model leaving the model intact. Keppenne (2000) applied this method to a large two-layer shallow water model. The main issue is in the analysis step, where the data from all ensemble members are required even though only one measurement is to be assimilated. The time spent on transfer between different processors will be increased. Usually there are two choices in the analysis step: member-decomposition (Keppenne, 2000) and domain-decomposition (Houtekamer and Mitchell, 1998). The latter is usually restricted to domain-decomposition parallelized model. In our case, the member-decomposition is used. Another point that should be noted is the random number generation process. Paralleled generator should be used instead of a serial one in order to avoid the strong correlation between random numbers generated by different processors. 3. Results Twin experiments were used to evaluate the performance of the EnKF as the DA methodology for the Monte Carlo dispersion model. Twin experiment is a sophisticated method for testing a DA algorithm, especially when there are no real observations. Even though there are some observation data, it is worthwhile to test a DA method by using simulated measured data before assimilating real ones, because in twin experiments, the model parameter error and measurement error can easily be controlled. Twin experiments include two parts: ‘truth’ run and simulation run. The ‘truth’ run uses the default meteorological parameter values and source terms, and provides a reference for the ‘true’ situation. Another task of the ‘truth’ run is to produce simulated observations by adding a random measurement error to the ‘true’ value. These simulated observations are to be assimilated in the simulation run. The simulation run is similar to the real assimilation process, except that simulated measurements, not real measurements, are to be assimilated. It differs from the ‘truth’ run in two ways—one is that the wind turbulence or source term is perturbed, which represents the error in the
ARTICLE IN PRESS D.Q. Zheng et al. / Atmospheric Environment 41 (2007) 2438–2446
estimated parameters, and the other is that the simulated measured data are assimilated in the model. Because in the real simulation procedure, the real source term and wind turbulence are not known, what we know is only an estimated value which includes errors in it. After the simulation run, the results can be compared with those results obtained by using the disturbed source term and meteorological parameters, but without DA, to see if the prediction is improved. In this paper, the forecast covers a period of over 24 h. The simulated measurements are assimilated every hour. We used the current observation network in south China and Hong Kong which is shown in Fig. 1. Simulated measured dose rate at the observation stations will be assimilated to see the impact of DA on the results.
2443
The simulated true source term and the source term used in the forecast model (in which a constant release rate is used) is shown in Fig. 2, in which the solid line represents the simulated true source term, while the dotted straight line represents the source term used in the deterministic forecast model. For assimilation, the simulated ‘truth’ is generated with a model using a randomly disturbed source term. The ground-level dose rate for all the observation points are also calculated with the ‘truth’ run, and the measurement are simulated from the ‘truth’ dose rate disturbed with a 5% measurement error. The measured data are assimilated every hour.
3.1. Uncertainty in source term In this experiment, only the nuclide release rate is considered uncertain. With an average release rate r in the deterministic model, the stochastic release rate was modeled according to rk ¼ maxð0; rð1 þ ck ÞÞ,
(13)
where ck is the colored noise which can be produced according to Eq. (10). The standard deviation of the error is 100% of r, and the time decorrelation length t is 2 h.
Release rate (Bq/s)
2.5×106
Forcecast source term True source term
2.0×106 1.5×106 1.0×106 5.0×105 0.0 0
2
4
6
8
10
Time (h) Fig. 2. Simulated true source term and the forecast source term.
Fig. 1. Observation stations in the model region. The star represents the location of power plant.
ARTICLE IN PRESS D.Q. Zheng et al. / Atmospheric Environment 41 (2007) 2438–2446
2444
2.0×10-6
3.5×10-6
True Without data assimilation With data assimilation
Doserate (Gy/s)
Doserate (Gy/s)
3.0×10-6 2.5×10-6 2.0×10-6 1.5×10-6 1.0×10-6
True Without data assimilation With data assimilation
1.5×10-6
1.0×10-6
5.0×10-7
5.0×10-7 0.0
0.0 0
(a)
5
10
15
20
Time(h)
0
5
10
15
20
Time (h)
(b)
Fig. 3. Dose rates from the true source term and the forecast source term (with and without assimilation) for location (a) 101 and (b) 602, respectively.
3.2. Uncertainty in the turbulence intensity of wind field Similar to the test that uses randomly perturbed source term, the perturbed turbulence parameter, here, the standard deviation of the turbulence velocity su , is used this time. As the error in the
6×106 Without assimilation With assimilation
5×106 RMS (Bq/m3)
Fig. 3(a) and (b) shows some assimilation experiment results with a randomly perturbed source term for two locations which are labeled with 101 and 602, respectively, in Fig. 1. The station 101 is at Da Peng town which is about 8 km west of the power plant and is the biggest town near Daya Bay plant. The station 602 is at Kat O island in Hong Kong which is about 25 km away from the power plant. In the deterministic forecast model, the release rate for the entire release period is constant, so the trend of dose rate coming out from the forecast model without assimilation is relatively flat. The forecast source term is corrected by assimilating the measurements which gives a dose rate closer to the true value than that without assimilation. The result in respect of the point 101 follows very closely the true value. In comparison, for the point 602, the estimated value is a little higher than the true value during times of relatively high dose rate, but the trend is in general similar to the true value. In order to see the overall impact of DA to the model prediction, the root mean square (RMS) of the difference between the estimated with the ‘true’ concentration is shown in Fig. 4. From the figure we can see that about 80% of the error that arise from the uncertain in the release rate is overcome. This result is encouraging.
4×106 3×106 2×106 1×106 0 0
5
10
15
20
25
Time (h) Fig. 4. RMS of the estimated vs. ‘true’ concentration with and without DA when release rate is taken as uncertain.
turbulence intensity influences the diffusion parameters directly, it is also an important parameter in the dispersion model. On the other side, because the parameter is estimated by PBL parameterizations, it is prone to involve errors in it. PBL parameterization is an simplification of the model physics. It may overestimate the turbulence intensity under a certain situation and underestimate it under another one. So, in this experiment, it is assumed that there is a bias in the turbulence estimation. In the determinate model, su is perturbed by increasing 30% of the ‘true’ value. The purpose of the experiment is to see if the DA method can work well in case that the forecast turbulence intensity is 30% higher than the true value. The comparison between the concentration distribution with and without DA is shown in Fig. 5. The forecast concentration without DA is found to be lower than the true value. This is
ARTICLE IN PRESS D.Q. Zheng et al. / Atmospheric Environment 41 (2007) 2438–2446
because the estimated turbulence strength is higher than the ‘true’ value. After DA, the difference from the ‘true’ concentration distribution is found to be
2445
1.2×106 Without assimilation With assimilstion
RMS (Bq/s3)
1.0×106 8.0×105 6.0×105 4.0×105 2.0×105 0.0 0
5
10
15 Time (h)
20
25
Fig. 6. RMS of the estimated vs. ’true’ concentration with and without DA when su is taken as uncertain.
reduced, implying that the assimilated measurements have updated the forecast concentration. The RMS of the difference between the estimated state and the ‘true’ state in this case is shown in Fig. 6. It shows that the error is reduced after DA, and the degree is about 50%.
4. Conclusion The results of the twin experiments show that successful application of the ensemble Kalman filter in combination of a nonlinear Monte Carlo dispersion model is possible. Dynamic selection of the state vector seems to be a feasible method to improve the efficiency of assimilation in case of accident. Introducing colored noise instead of white noise is more reasonable, and EnKF can work well under this situation. Although EnKF provides a possibility of combining model and measurement together, much work must be done before the methodology can be applied in a real accident emergency system. Here are some recommendations:
Fig. 5. Contours of radionuclide concentration 10 h after release with and without assimilation. (a) The true situation. (b) The forecast alone with perturbed turbulence. (c) The estimation after the measured dose rate is assimilated.
(i) Other uncertainties except source term and turbulence intensity (i.e. uncertainties in wind direction, source height estimations, etc.) should be taken into account, and the attempt with various uncertainties acting together should be done. (ii) Twin experiment only provides a way for validating the methodology in a simulation way. It is desirable to validate it against realcase measurements.
ARTICLE IN PRESS 2446
D.Q. Zheng et al. / Atmospheric Environment 41 (2007) 2438–2446
(iii) Considering that the real nuclear accident measurement data are difficult to be obtained, wind tunnel experiments can be used to get some real measurements. References Evensen, G., 1994. Sequential data assimilation with a non-linear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. Journal of Geophysical Research 99, 10,143–10,162. Evensen, G., 2003. The ensemble Kalman filter: theoretical formulation and practical implementation. Ocean Dynamics 53, 343–367. Evensen, G., Madsen, H., Gering, F., 1997. Advanced data assimilation for strongly non-linear dynamics. Monthly Weather Review 25, 1342–1354. Grewal, M.S., Andrews, A.P., 1993. Kalman Filtering: Theory and Practice. Prentice-Hall, Englewood Cliffs, NJ. Houtekamer, P.L., Mitchell, H.L., 1998. Data assimilation using an ensemble Kalman filter technique. Monthly Weather Review 126, 796–811. Ide, K., Ghil, M., 1997. Extended Kalman filtering for vortex systems. Part I: methodology and point vortices. Dynamics of Atmospheres and Oceans 27, 301–332. Keppenne, C.L., 2000. Data assimilation into a primitiveequation model with a parallel ensemble Kalman filter. Monthly Weather Review 128, 1971–1981. Kovalets, I., Andronopoulos, S., Bartzis, J.G., Gounaris, N., Kushchan, A., 2003. Introduction of data assimilation
procedures in the meteorological pre-processor of atmospheric dispersion models used in emergency response systems. Atmospheric Environment 38, 457–467. Leung, J.K.C., 1992. Application of shielding factors for protection against gamma radiations during a nuclear accident. IEEE Transactions on Nuclear Science 39, 1512–1518. Mikkelsen, T., Larsen, S.E., Thykier-Nielsen, S., 1984. Description of the RisøPuff diffusion model. Nuclear Technology 67, 56–65. Rojas-Palma, C., Madsen, H., Cering, F., Puch, R., Turcanu, C., Astrup, P., Mu¨ller, H., Richter, K., Zheleznyak, M., Treebushny, D., Kolomeev, M., Kamaev, D., Wynn, H., 2003. Data assimilation in the decision support system RODOS. Radiation Protection Dosimetry 104, 31–40. Shleien, B., Slaback Jr., L.A., Birky, B.K., 1998. Handbook of Health Physics and Radiological Health. Williams & Wilkins, Baltimore, MD. van Loon, M., Builtjes, P.J.H., Segers, A.J., 2000. Data assimilation of ozone in the atmospheric transport chemistry model LOTOS. Environmental Modeling and Software 15, 603–609. Zannettip, P., 1990. Air pollution modelling: theories, computational method and available software. Van Nostrand Reinhold, New York. Zhang, M.G., Han, Z.W., Lei, X.E., 1999. Application of Monte Carlo model for multiple sources to air quality assessment at atomic power plant of Guangdong. Climatic and Environmental Research 14, 203–209. Zhang, M.G., Han, Z.W., Lei, X.E., 2000. A numerical study of atmospheric dispersion in the area of nuclear power plant Dayawan. Climatic and Environmental Research 5, 90–95.