# Probabilistic Dose Assessment from SB-LOCA Accident in Ujung Lemahabang Using TMI-2 Source Term

Authors: Sunarko; Zaki Su’ud; Idam Arif; Syeilendra Pramuditya
Journal: Science and Technology of Nuclear Installations
 (2017)

---

## Abstract

Probabilistic dose assessment and mapping for nuclear accident condition are performed for Ujung Lemahabang site in Muria Peninsula region in Indonesia. Source term is obtained from Three-Mile Island unit 2 (TMI-2) PWR-type SB-LOCA reactor accident inverse modeling. Effluent consisted of Xe-133, Kr-88, I-131, and Cs-137 released from a 50 m stack. Lagrangian Particle Dispersion Method (LPDM) and 3-dimensional mass-consistent wind field are employed to obtain surface-level time-integrated air concentration and spatial distribution of ground-level total dose in dry condition. Site-specific meteorological data is obtained from hourly records obtained during the Site Feasibility Study period in Ujung Lemahabang. Effluent is released from a height of 50 meters in uniform rate during a 6-hour period and the dose is integrated during this period in a neutrally stable atmospheric condition. Maximum dose noted is below regulatory limit of 1 mSv and radioactive plume is spread mostly to the W-SW inland and to N-NE from the proposed plant to Java Sea. This paper has demonstrated for the first time a probabilistic analysis method for assessing possible spatial dose distribution, a hypothetical release, and a set of meteorological data for Ujung Lemahabang region.


---

## Body

## 1. Introduction

The capability to estimate radioactive dose dispersion in a course of an accident in a nuclear facility is important to support emergency planning activities and ensure public safety from ionizing radiation exposure as a consequence of accidental or normal releases. In nuclear power plant site study stage, the information contained in this analysis can be used as an input in emergency response planning, radiation protection activities, and nuclear power plant siting and design. Atmospheric dispersion models are typically applied in site evaluation to measure short-term to long-term normalized concentration and deposition. Other than undergoing advection and diffusion, the effluent may also experience radioactive decay as well as wet deposition and dry deposition. These effects can be expressed mathematically and should be considered in the model appropriately. Whenever possible, site and/or plant specific characteristics should be considered in the model‎ [1].A computer code has been developed at Institut Teknologi Bandung based on Lagrangian Particle Dispersion Method ‎[2–5]. The code is capable of handling spatial and temporal variation in topography, wind, and source term. The inputs for this code include topography, surface wind station data, and source term which are usually available during the site study phase through survey, monitoring, or other means of data collection. The code is intended for power plant site study purpose and covers an area in the scale of up to a few hundreds of kilometers around the plant.Site Feasibility Study for nuclear power plant performed for Muria Peninsula in Central Java has resulted in three candidate sites: Ujung Lemahabang, Ujung Genggrengan, and Ujungwatu. All of the candidate sites are close to population and therefore it is necessary to prepare for emergency planning actions in case of an accident.This paper aims to evaluate, based on site-specific meteorological characteristics, the spatial distribution of doses and their probabilities as a result of an SB-LOCA accident similar to the Three-Mile Island unit 2 (TMI-2) nuclear accident. The source term data is taken from the Tree-Mile Island unit 2 (TMI-2) nuclear accident which according to the Fact Sheet of the US Nuclear Regulatory Committee represents the most serious nuclear accident in the US commercial nuclear power plant operating history [6]. In our previous study we use the worst scenario in which the wind is directed toward the land for several hours consistently [7]. Therefore the resulting dose will become the largest. In this paper we consider about more than 100 possible wind scenarios in stochastic approach. Therefore the current paper shows common condition estimate.


## 2. Lagrangian Particle Dispersion Method

Concentration of radionuclides in the environment is modeled based on the advection-diffusion equation (1) [13]. The equation is numerically solved using random-walk Lagrangian Particle Dispersion Method (LPDM) with a 3-dimensional wind field in Cartesian coordinate that can accommodate spatially variable wind speed, diffusion coefficient, and topographic features of the area under investigation. Empiric dispersion coefficient is used for the advection-diffusion equation (1).ParameterC- is the average concentration, u-, v-, and w- are wind components in x, y, and z directions, Kx, Ky, and Kz are turbulent diffusion coefficients for x, y, and z, and ws is the gravitational settling velocity. Λ is the depletion coefficient, λ is the decay constant, and Q is the source term. In this paper, gravitational settling, depletion, and decay are not taken into account.(1)∂C-∂t=-u-∂C-∂x-v-∂C-∂y-w-∂C-∂z+∂∂xKx∂C-∂x+∂∂yKy∂C-∂y+∂∂zKz∂C-∂z+ws∂C-∂z-ΛC--λC-+Q.Atmospheric turbulence is inherently stochastic and many researchers proposed an approach based on the statistical nature of turbulence. In stochastic representation, marked particles undergo advection process by wind at a certain speed and at the same time experience random movement simulating turbulent fluctuation. Average distribution of particles is determined by averaging particle paths. Since every particle moves independently, simultaneous handling of particles is unnecessary and therefore requires small computer memory‎ [14]. The accuracy of this method increases with smaller computational volume.Lagrangian Particle Dispersion Method (LPDM) is run by tracking a number of particles in a wind field. LPDM uses stochastic differential equation to explain the similar process as the advection-diffusion equation in Lagrangian framework‎ [15]. Stochastic differential equation for the movement of ideal fluid in three dimensions is(2)dxi=u¯dt+2Kx1/2dWx,(3)dyi=v¯dt+2Ky1/2dWy,(4)dzi=w¯dt+2Kz1/2dWz+∂Kz∂z.ParameterdWx,y,z is random numbers from Gaussian distribution with zero average and variance dt; namely, dW¯=0 and dW¯2=dt. In (2) and (3) it is assumed that turbulence is homogeneous in x and y directions. The above equations can be integrated with time to obtain particle path which represents the movement of each individual particle. In numerical calculation of dispersion, the number of released particles is large, each with its own label and carrying a certain characteristic such as mass or radioactivity. Concentration at time t can be calculated from locations of particles and the particle characteristics.Implementation of random-walk method is explained for vertical direction as follows. The equation for horizontal direction has the same form as the vertical direction but without the differential form of eddy diffusivity. The height of a particle after one time stepzi+1 is a summation of four terms, namely, the initial height of the particle zi, movement by wind w-Δt, average diffusive movement Δzi∗, and random diffusive movement Δzi′, as follows:(5)zi+1=zi+w¯Δt+Δzi∗+Δzi′.Average diffusive movement is represented as follows:(6)Δzi∗=∂Kz∂ziΔti.Random diffusive movement has the following average and variance:(7)Δzi′=0,(8)Δzi′2=2KziΔti+∂Kz∂z2Δti2≡σzi2.Random diffusive movementΔzi′ is obtained using uniform distribution function in (9) ‎[15]. This distribution function is used so that (10) can use uniform random number to replace the Gaussian distributed random number.(9)PΔzi′=123σziif-3σzi<Δzi′<3σzi0others.The diffusive displacementΔzi′ is given by(10)Δzi′=23σzir-0.5,where r is now a random number from a pseudorandom number generator with a range of (0,1).Dispersion coefficients for the model are provided in Diehl et al. (1982). Vertical dispersion coefficient is given for unstable and stable condition in (11) and (12), respectively.(11)Kz=0.1u∗hhkL1/3sin⁡z/h+0.5π21-exp⁡-4zh-0.0003exp⁡8zh,(12)Kz=0.065u∗hzh2/31-zh2/31/2.Horizontal dispersion coefficient for unstable and stable condition is given in (13) and (14), respectively.(13)KH=0.13u∗h12+0.5hl1/3unstable,(14)KH=0.14u∗hzh1/21-zh1/21/2stable.Equations (12) and (14) also work for neutral atmospheric condition ‎[14]. In the above equations, h is the boundary layer thickness and L is the Monin-Obukhov length. The value for u∗ is determined according to the surface roughness. L can be found for Pasquill stability classes A–F by the relationship in (15)‎ [8] and coefficients provided in Table 1:(15)1L=a+blog⁡z0.The value of z0 for Ujung Lemahabang area is assumed to be 1 (tree-covered surface) since the site area is located in a rubber plantation.Table 1
Coefficients for Monin-Obukhov length‎ [8].
Pasquillstability class
Coefficient
Description
a
b
A
-
0,096
0,029
Highly unstable
B
-
0,037
0,029
Unstable
C
-
0,002
0,018
Slightly unstable
D
0
0
Neutral
E
+
0,004
-
0,018
Slightly stable
F
+
0,035
-
0,036
StableTime-integrated air concentration (TIAC) is obtained by calculating concentration in one time step and integrating throughout simulation time. Concentration is calculated using kernel-density estimation method for surface level (0–5 m layer). Total Effective Dose Equivalent (TEDE) is then calculated using tabulated values of dose conversion factors provided in Lam et al.‎ [16] using the previously obtained TIAC. The values of Monin-Obukhov lengths used in this simulations are shown in Table 1.


## 3. Mass-Consistent Wind Field

Wind field is provided using mass-consistent model‎ [17]. The method is used to provide wind vectors in a 3-dimensional Cartesian grid of horizontal intervals Δx and Δy of 1 km and vertical intervals Δz of 50 meters covering an area 100 km × 100 km wide and 1500 meter high. Block topography is constructed using data from ground elevation survey. Initial vertical and horizontal wind vector are provided by the metrological data from Ujung Lemahabang station and iterations are performed until convergence is reached. Initial vertical wind distribution is provided by the power law formula ‎[9, 17]:(16)u=u0zz0p.The powerp is determined from least-square fit of multiple-level tower data or from atmospheric stability and terrain conditions as given in Table 2.Table 2
The powerp value for wind-speed profile‎ [9].
Stability category
Village exponent
Urban exponent
A
0.07
0.15
B
0.07
0.15
C
0.10
0.20
D
0.15
0.25
E
0.35
0.30
F
0.55
0.30


## 4. Three-Mile Island Accident Source Term

Three-Mile Island unit 2 underwent an accident on March 28, 1979. The accident was initiated by a loss of feed water resulting in pressure increase in the secondary system and turbine and main feed water pump trip. Pressure increase in the primary system caused the pressurizer pilot operated relief valve (PORV) to open and reactor scram. The PORV valve failed to close causing a small-break loss-of-coolant accident (SB-LOCA). Cooling to the core was performed through forced circulation until the reactor coolant pumps tripped due to high local void ‎[10]. Part of the core was melted and fission products were released from the core to the primary containment which were then transported by the Makeup and Purification System to the auxiliary building. Outgassing from the transported cooling water has caused fission product gasses to enter the auxiliary building atmosphere and ventilation system and allowed some noble gasses and radioiodine to escape to the environment‎ [12].Several attempts have been made to estimate the released radioactive gasses including the usage of reverse modeling to estimate the source strength based on the data from monitors around the site. It is estimated that the radioactive gasses were released during a two-week period with 90% of the noble gasses being released on the first three days after the accident‎ [12]. Approximately 2.4–10 million curies of noble gasses (mainly Xe-133) and about 14 curies of I-131 were released. Simulation with MELCOR/MACCS estimated a total mass release of about 0.1 kg of noble gasses and 2.5 × 10−7kg of CsI‎ [10]. This is an underestimation by 1000-fold and 15-fold, respectively, compared to the actual release estimate. The radionuclide is assumed to be released in a 6-hour period from 8.500 seconds after the onset of the incident to 30.000 seconds after.Simulations are performed using source term for isotopes of Cs-137, I-131, Kr-88, and Xe-133 each with a total strength of 0.011, 15.04, 62.000, and 8.37 (all units in×106 curies)‎ [11]. The plume is released at a constant rate at seconds 8.500 to 30.000. For conservative measure, the maximum values for source strength in Table 3 were chosen as the source term to represent the hypothetical release.Table 3
Source term.
Source strength rate (Bq/sec)
Radionuclide
Haste et al. ‎[10]
McColl and Prosser [11]
Gudiksen and Dickerson ‎[12]
US-NRC
Maximum
Cs-137
2.27
E
+
04
—
—
—
2.27
E
+
04
I-131
3.12
E
+
07
2.31
E
+
07
2.40
E
+
07
—
3.12
E
+
07
Kr-88
—
1.06
E
+
08
—
—
1.06
E
+
08
Xe-133
—
1.44
E
+
10
4.11
E
+
12
4.28
E
+
12
4.28
E
+
12


## 5. Meteorological Data

Meteorological data used in this paper is obtained from on-site monitoring performed in Ujung Lemahabang (ULA) site in Muria Peninsula, Jepara Regency in Central Java province for around 1-year duration (August 1994 to August 1995)‎ [18] during the Site Feasibility Study period as this is the only site-specific and reliable data available for Ujung Lemahabang area. The data is composed of temperature (10 and 50 meters) and wind speed and direction (10 and 40 meters). Data is logged every hour and organized into columns (00:00 hour to 23:00 hour) and rows (days). Temperature data is used to determine atmospheric stability class based on lapse-rate method. In the simulation, a number of time periods are determined by randomly choosing the beginning of the period from the database.


## 6. Methodology

A computer code has been written to perform the analysis. Input data includes meteorology, source term, and topography around the site. Prognostic wind field in Cartesian grid is created using mass-consistent method utilizing 1-station data for a dimension of 100 km × 100 km × 1.5 km (width, length, and height, resp.). A plume consisting of 18000 particles per hour is used in the model. Advection by wind field and dispersion by eddies are simulated using LPDM. Concentrations at grid points are calculated using kernel-density estimator for the lowest 5-meter layer of the atmosphere. Concentrations are then time-integrated to obtain TIAC which is then converted to dose using tabulated conversion factor‎ [16].For a specific release scenario the probability of upper-level dose can be mapped to provide a contour of probability for exceedance of attaining a certain dose limit. The limit may be in form of regulatory limit. The resulting map describes the probability of attaining a certain dose limit around a facility based on wind profile statisticsP[Dx,y≥doselimit∣release]. Release parameter is deterministic. D(x,y) is the dose received at coordinate (x,y) and the dose limit can be based on values determined by regulatory body. For a number of wind profiles Nw the probability P(x,y) at node (x,y) is determined by‎ [19](17)Px,y=1Nw∑i=1Nwni,ni=1,ifDix,y≥dose  limit∣release0,other.i is the index for wind profile data. For each wind profile a dose map is created according to the process flow in Figure 1 based on the source term for the same period. The overall number of wind profiles Nw is equivalent to the number of simulation runs and each simulation is independent of each other since the result of one simulation is not affected by the result of another simulation. The probability map is useful to determining the upper limit of dose received if the release parameter is assumed to be the maximum-release scenario.Figure 1
Process flow.


## 7. Result and Discussion

Meteorological data record from monitoring height of 10 meters shows that wind is distributed into all 16 sectors with the highest frequency of 10.4% from the southwest direction. Wind speed is dominated by 0.5–2.1 m/s class with a frequency of 44.8%. Calm wind condition consists of 3.0% of the total. Average wind speed at this level is 2.63 m/s.On the other hand, at 40-meter level wind direction is dominated by south direction with 9.1%. Wind speed is dominated by 3.6–5.7 m/s class with a frequency of 32.4%. Calm condition is noted at 0.5% while average wind speed is recorded at 5.46 m/s. Wind condition can be summarized by the wind rose plot for the monitoring period of August 23, 1994, to August 31, 1995, in Figure2.Figure 2
Wind rose in ULA at heights of 10 m (a) and 40 m (b).
(a)
(b)In the simulation, the maximum total effective dose never exceeded the regulatory limit of 1 mSv for general public in one year ‎[20] and within the range of field measurement performed around TMI-2 NPP at the time of accident. This is also much less than the equivalent dose for skin for general public in one year of 50 mSv ‎[20]. The peak external dose based on measurement in areas surrounding TMI-2 NPP was around 1 mSv. Therefore, the result is presented for other values to provide knowledge on the pattern of plume dispersion for 0.5 mSv and 0.001 mSv, provided in Figures 3 and 4, respectively, which is below the regulated limit.Figure 3
Dose distribution probability (D≥0.5 mSv).Figure 4
Dose distribution probability (D≥0.001 mSv).It appears that on land the radioactive plume tends to be spread from the plant to N-NE directions and to W-SW direction from the proposed plant location, consistent with wind distribution pattern at both 10 m and 40 m levels. Almost the same amount of plume will be dispersed to the N-NE direction toward the Java Sea. The villages most likely to receive 0.001 mSv or more by an incident in this wind condition are the northern part of Balong and small part of Tuban with a probability of more than 0.57%. Other villages to the west of Balong are partly affected such as Kaliaman, Bondo, Kedungleper, and Bangsri but with smaller probability. To the east of Balong village, the northern part of Bumiharjo and Bandungharjo may receive the same amount of dose. The northern part of Balong village is a state-owned rubber plantation area of PT. Perkebunan Nusantara and only a small number of people are living in the area.In the case of larger amount of radioactive release or in the presence of wet deposition, the magnitude of dose received will be higher, extending to farther areas, and may exceed the regulatory limit. However, it is very likely that the direction of the dispersion will be similar to the dry condition. In order to mitigate the effect from more severe accident, shelters and evacuation routes can be prepared, if they are not readily available, from the shore part of the affected villages to the main road connecting Jepara to the west and Pati regency to the east. Radiation monitoring posts should be placed in the areas most likely to be affected and early warning system should also be prepared to inform the public of possible emergencies in those areas, especially at the northern part of the peninsula to the west-southwest direction from the assumed release point.


## 8. Conclusion

Probabilistic analysis using upper limit scenario has been applied to make individual dose prediction for a no-precipitation case in Muria Peninsula based on source term parameters from Three-Mile Island NPP accident in the USA and meteorological database from on-site monitoring.Regulatory limit of 1 mSv for effective dose and 50 mSv for equivalent skin dose was never reached in the simulation and maps of different values were used to provide an illustration on the trend of dispersion based on the local meteorological behavior. The dispersion pattern is consistent with the meteorological characteristics of the site. On land, effluent has a tendency to be dispersed to the W-SW direction which is consistent with the dominant wind pattern for the area. Northern part of Balong village may receive≥0.5 mSv with a probability of around 18%. The area receiving the same amount of dose is well within a 2.5 km radius from the assumed release point.Better description of plume dispersion pattern is given by the ≥0.001 mSv dose limit. The villages along the shore to the west of Ujung Lemahabang site such as Tuban, Kaliaman, Kedungleper, and Bondo are more likely to be affected than other areas, although the amount of dose received is much less than the regulated dose. The received dose in this case is well below the limit set by Indonesian nuclear regulatory body.



---

*Source: 8761397-2017-08-09.xml*