Transformation of a system consisting of plane isotropic source and unit sphere detector into a system consisting of point isotropic source and plane detector in Monte Carlo radiation transport calculation

Abstract

In a nuclear power plant accident, radioactive nuclides may be released which are distributed uniformly on the ground. If estimation of dose rate from such a source by a Monte Carlo calculation is attempted, some difficulty is encountered because the calculation efficiency is very low. To solve this low efficiency problem, we show that a plane isotropic source can be transformed into a point isotropic source by changing the detector shape from a unit sphere to a plane. We verified the validity of this transformation by the numerical comparison of unscattered photon fluence. As an example of this transformation, the ambient dose rate D _i was calculated from the uniform radioactive nuclide distribution on the ground using the EGS5 Monte Carlo code. We also measured the radioactivity and ambient dose rate (M) on the KEK campus within a month after the releases from the Fukushima No. 1 Nuclear Power Plant accident. Using radioactivity data and D _i, we calculated the ambient dose rate (C). The calculated and measured ambient dose rates agreed reasonably well; their ratio (C/M) was 0.62 to 1.28.

Keywords :

• plane isotropic source
• unit sphere detector
• point isotropic source
• plane detector
• Monte Carlo
• EGS5
• transformation
• ambient dose rate, photon
• Fukushima No.1 nuclear power plant accident

1. Introduction

Monte Carlo electron–photon transport codes such as EGS5 [1] are very useful tools for solving many types of radiation transport problems. However, it is difficult to use a Monte Carlo program to solve a problem that involves scoring photons from a large plane source using a unit sphere detector because the calculation efficiency is very low and the calculation time becomes very large. Most of the photons that are generated at position far away from the unit sphere detector have only very small contribution to detector response as the contribution is proportional to the inverse of square of the distance between the source and the detector. Consequently, the size of the plane source is restricted by the limit placed on the calculation time. It is effectively impossible to treat an infinitely large plane source with a unit sphere detector in a Monte Carlo calculation.

If a system consisting of a plane source and unit sphere detector is transformed into a system consisting of a point source and plane detector, the calculation efficiency problem can be solved. In the case of parallel beams, it is well known that this transformation is possible and it has been used in many situations. For example, Chilton et al. [2] state that “it is intuitively obvious that this kind of geometric transformation gives the same expected result.” In addition, Shultis and Faw [3] state that the dose is preserved after the geometric transformation, which involves the exchanging of the source region and detector region, as long as the total source intensity is preserved. However, no specific equations, which can be utilized for programming a practical Monte Carlo code have been provided for individual geometries.

In a nuclear power plant accident, radioactive nuclide may be released which distribute uniformly on the ground. If the relation between the radioactivity and the ambient dose is known, either one can be estimated by measuring the other. However, it is difficult to determine this relation using the Monte Carlo code, as described above. Therefore, in the first part of this article, we show that the transformation of a system consisting of a plane isotropic source and unit sphere detector into a system consisting of a point isotropic source and plane detector is possible. Then, we describe the calculation of the ambient dose from a plane isotropic source for radioactive nuclides, which were released from the Fukushima No. 1 Nuclear Power Plant, as an example of the geometrical transformation. We measured the ambient dose rate and intensity of radioactivity in KEK. We compared the measured and calculated ambient dose rates.

2. Calculation

2.1. Plane isotropic source and unit sphere detector

2.1.1. Without shield

To estimate the photon flux generated from a plane isotropic source, we consider a disc isotropic source of radius r ₁ and central point O, as shown in Figure 1 . We set the source intensity per unit area and total source intensity as I (cm⁻² s⁻¹) and (s⁻¹), respectively. We place a unit sphere detector at point V which is located on the central axis of the disc.

Figure 1. System consisting of plane isotropic source and unit sphere detector.

The partial photon flux dφ at V due to a photon from a tiny area dS of the disc is

The total photon flux at V is

By using the relations, tan θ = r/h and dr/dθ = h/cos²θ, we get

Here, we use the variable u = cos θ:

2.1.2. With shield

We consider that a slab shield of thickness t and attenuation coefficient μ exists between the disc source and point V. As the photon path in this shield is equal to , the photon fluence at V is

In Equation (2), exp{−μt/cos θ} and δ_s on the right-hand-side represent unscattered and scattered components, respectively. B(μt, cos θ) represents the increase of photon fluence due to scattering.

Here, we use the relation tan θ = r/h as in the case in which no shield is used:

Then, we introduce the exponential integral E ₁,

By setting B = 1, we can represent the unscattered part φ as

Most of the equations presented in this subsection have already been described in standard radiation shielding textbooks such as [4]. We have presented these equations so that the reader can easily follow the derivation in the next subsection.

2.2. Point isotropic source and plane detector

2.2.1. Without shield

We exchange the form of the source and the detector, i.e. we consider a system consisting of a unit sphere isotropic source and a plane detector. As a unit sphere source and a point source are equivalent in a Monte Carlo calculation, we change the source shape from a unit sphere to a point to get the system shown in Figure 2 . We use a disc-shaped detector of radius r ₁ and a photon source at point V. The intensity of the photon source is set as identical to the total source intensity of the plane isotropic source, I _tot (S⁻¹). We consider a solid angle dS′, which is defined as a product of sin θ dθ and dψ. As the source is emitted isotropically, the photon intensity emitted from V into the solid angle dS′ is

Figure 2. System consisting of point isotropic source and plane detector.

In the system before the geometry transformation, we use the unit sphere detector. In this case, output of photon fluence does not depend on photon incident direction. In this article, we discuss exchanging the shapes of the source and detector. But the function of the detector should be maintained after the transformation. Accordingly, the plane detector must have a function that output of photon fluence does not depend on photon incident direction. When a photon is incident on the plane detector normally, photon fluence is obtained by dividing the incident number of photons by the area of the plane detector. When the plane detector is placed obliquely to the photon direction, the value which is obtained by dividing the incident number of photons by the area of the plane detector is cos θ times smaller than that in the normal incident case. Accordingly, we have to divide the incident number of photons by cos θ so that output of photon fluence does not depend on incident direction when photons are incident to the detector from an arbitrary direction.

To get photon fluence at the position of the plane detector, we integrate dφ′/cos θ over the detector region and divide it by the detector size :

Here, tan θ₁ = r ₁/h, as shown in Figure 2. We substitute v = cos θ and dv/dθ = −sin θ:

We use the relations

Then,

This equation is identical to Equation (1). Thus, the two systems are mathematically identical.

2.2.2. With shield

The photon intensity that is emitted from V toward dS′ and reaches dS′ is given as

Here, represents the contribution of scattering. After dividing dφ′ by cos θ, we integrate dφ′ and obtain

Here, B′(μt, cos θ) represents the increase of photon fluence due to scattering. As the structure of the shield is not changed in the transformation, B(μt, cos θ) = B′(μt, cos θ). Then this equation is identical to Equation (3). Thus, the two systems are mathematically identical when we have a slab shield between the source and the detector.

So far, we have considered a transformation for a plane source. However, this transformation is also possible for a “modified” plane source that has some distribution in the depth direction because the overlapping of the shield and source regions does not affect the validity of the discussion presented above.

2.3. Numerical comparison

To verify the validity of this transformation, we performed a numerical comparison. We compared the values obtained from an analytical equation (Equation (1) and Equation (4) for the unshielded and shielded cases, respectively) and from the Monte Carlo calculation. In the former, the plane isotropic source with uniform distribution was assumed. In the EGS5 calculation, the point isotropic source was specified and only unscattered photons were counted with the plane detector, which corresponds to the system after the transformation. We used 365 keV photons from I–131 for the comparison, and either air or iron was used as shielding material. As listed in Tables 1 and 2 for the unshielded and shielded cases, respectively, the results of the analytical calculation and those for the unscattered photons from the EGS5 Monte Carlo calculation agreed well. Accordingly, the validity of the transformation was clearly verified by numerical comparison.

Table 1. Comparison of photon fluence (unshielded case).

h (cm)	r 1 (m)	Equation (1)	EGS5
10	1	1.154	1.156 ± 9.7e−4*
	10	2.303	2.299 ± 4.4e−3
	100	3.454	3.436 ± 1.5e−2
50	1	0.402	0.402 ± 3.8e−4
	10	1.498	1.498 ± 1.6e−3
	100	2.649	2.644 ± 6.5e−3
100	1	0.173	0.173 ± 2.7e−4
	10	1.154	1.156 ± 9.7e−4
	100	2.303	2.299 ± 4.4e−3
150	1	0.0919	0.0916 ± 2.1e−4
	10	0.954	0.954 ± 7.2e−4
	100	2.100	2.096 ± 3.5e−3
Note: Read as 9.7 × 10^−4.

Table 2. Comparison of unscattered photon fluence (shielded case).

h (cm)	μh*	Equation (4)	EGS5#
Air shield
10	1.19e−3	3.079	3.073 ± 9.7e−3
50	5.95e−3	2.277	2.275 ± 4.0e−3
100	1.19e−2	1.933	1.932 ± 2.6e−3
150	1.78e−2	1.735	1.733 ± 1.9e−3
Fe shield
1	7.56e−1	0.168	0.168 ± 9.0e−5
2	1.51e + 0	0.0493	0.0491 ± 4.0e−5
3	2.27e + 0	0.0169	0.0169 ± 1.7e−4
4	3.02e + 0	6.36e−3	6.32e−3 ± 7.6e−6
Note: *μ is photon attenuation coefficient. ^#Only unscattered photons were scored.

3. Example application

As an example of the transformation, we calculated the air kerma K _i and the ambient dose rate D _i using the general-purpose electron–photon Monte Carlo transport code EGS5 [1]. The source was radioactive nuclides released by the Fukushima No. 1 Nuclear Power Plant in the accident following the 11 March 2011 earthquake and tsunami; nuclides were distributed uniformly on the ground. The calculation geometry in the original problem and that after transformation is shown in Figure 3(a) and (b), respectively. The calculation geometry consisted of a soil layer of 6 m thickness and an air layer of 2000 m thickness. The soil composition and density (1.6 g/cm³) were taken from [5]. The horizontal size of both the layers was infinite. A radioactive source with a radioactivity of unit strength (1 Bq/cm²) was placed in the soil surface. I-131, I-132, Cs-134, Cs-137, and Te-132 were the sources considered. The gamma ray energies and the emission ratios were taken from Ref. [6]. These were major radioactive nuclides which we observed on the KEK campus after the power plant accident. We ignored Cs-136 because its contribution was small. The number of incident photons was set as 10⁵ for each radioactive nuclide. The calculation results for air kerma K _i and the ambient dose rate D _i for a radioactivity of unit strength (1 Bq/cm²) on the ground surface are listed in Tables 3 and 4, respectively. The heights of the estimator were set as 10, 50, 100, and 150 cm. The standard deviations of obtained air kerma and the ambient dose rate were less than 2% and 1% for 10 cm and 100 cm, respectively. Both air kerma and the ambient dose rate decreased as the estimator height was increased.

Figure 3. Calculation geometry: (a) original problem and (b) transformed problem.

Table 3. Conversion factor for air kerma, K _i (μGy/h per Bq/cm²), calculated by EGS5.

h (cm)	Te-132	I-131	I-132	Cs-134	Cs-137
10	0.0191 ± 2.9e−4	0.0249 ± 5.2e−4	0.119 ± 2.3e−3	0.0987 ± 2.0e−3	0.0381 ± 7.6e−4
50	0.0148 ± 1.4e−4	0.0193 ± 2.0e−4	0.0940 ± 1.1e−3	0.0783 ± 9.3e−4	0.0291 ± 3.3e−4
100	0.0128 ± 1.0e−4	0.0171 ± 1.5e−4	0.0840 ± 9.1e−4	0.0694 ± 7.1e−4	0.0258 ± 2.4e−4
150	0.0119 ± 9.5e−5	0.0158 ± 1.3e−4	0.0764 ± 6.6e−4	0.0636 ± 6.1e−4	0.0239 ± 2.0e−4

Table 4. Conversion factor for ambient dose rate, D _i (μSv/h per Bq/cm²), calculated by EGS5.

h (cm)	Te-132	I-131	I-132	Cs-134	Cs-137
10	0.0249 ± 3.8e−4	0.0323 ± 6.6e−4	0.145 ± 2.8e−3	0.121 ± 2.4e−3	0.0466 ± 9.1e−4
50	0.0195 ± 1.9e−4	0.0252 ± 2.6e−4	0.115 ± 1.3e−3	0.0961 ± 1.1e−3	0.0358 ± 4.0e−4
100	0.0171 ± 1.4e−4	0.0224 ± 1.8e−4	0.103 ± 1.1e−3	0.0853 ± 8.6e−4	0.0318 ± 2.9e−4
150	0.0158 ± 1.3e−4	0.0208 ± 1.6e−4	0.0937 ± 7.9e−4	0.0784 ± 7.4e−4	0.0295 ± 2.5e−4

The ambient dose rates were measured on the KEK campus on 1 April and 11 April 2011, using NaI gamma ray survey meters (Aloka, TCS-161 and TCS-171) near the ground surface and at a height of 1 m. On 1 April, measurements were made at one place and on 11 April, they were made at five places. We calibrated the survey meters using Cs-137 for ambient dose. On the same day and at the same places where the ambient dose rate measurements were made, we also collected soil samples for radioactivity measurement. Each soil sample was 5.7 cm in diameter and 0.5–2.6 g/cm² in depth. We used a portable Ge detector (Canberra) for the radioactivity measurements and ISOCS software for determining the efficiency. The decay correction was applied in a separate calculation. On 1 April, we collected soil samples of two layers, as shown in Figure 4 . Radioactivity in the second layer was 2% of the sum of two layers in the measurement of 1 April. The measured radioactivities are listed in Table 5 .

Figure 4. Geometry of soil sample. The sample size was φ 5.7 cm × 0.5–2.6 g/cm².

Table 5. Measured radioactivity concentration on the ground in the KEK campus (Bq/cm²).

Date	1 April	11 April
Point*	0	1	2	3	4	5
Te-132	0.09
I-131	2.32	0.690	0.458	0.0130	1.16	2.21
I-132	0.07
Cs-134	0.28	0.382			0.181	0.413
Cs-136	0.02
Cs-137	0.27	0.436	0.089	0.0151	0.247	0.452
*Numbers indicate measurement position on KEK campus.

We obtained the calculated value of the ambient dose rate (C) as

Here, D _i and S _i are the ambient dose rate per unit radioactivity (μSv/h per Bq/cm²; calculated) and the radioactivity on the ground (Bq/cm²; measured), respectively. i indicates the type of radioactive nuclide. Bg is the intensity of natural background, which was assumed to be 0.06 μSv/h based on the radiation dose level in operational measurements before the 11 March earthquake. The calculated (C) and measured (M) ambient dose rates are compared in Figure 5 . In the ambient dose rate measurement on 1 April, we did not record the calibration factor and the serial number of the survey meter. Then, we selected smallest and largest calibration factors of the survey meters of the same type in our department to get the possible area of the measured value. The end point of the error bar in Figure 5 corresponds to the data on 1 April with these two calibration factors. The C/M ratio was 0.62–1.28; the agreement between C and M was reasonably good. In addition, the C/M ratios near the ground surface and at 1 m height were close to each other; therefore, the calculated D _i was capable of reproducing the ambient dose rate at both heights.

Figure 5. Comparison of calculated and measured ambient dose rates. Calculated values were obtained using the EGS5 code with a point isotropic source.

4. Conclusion

We showed that a plane isotropic source can be transformed into a point isotropic source by changing the detector shape from a unit sphere to a plane in Monte Carlo radiation transport calculation. We verified the validity of this transformation by the numerical comparison of unscattered photon fluence.

We showed that the relation of radioactivity density on the ground and the ambient dose after the Fukushima No. 1 Nuclear Power Plant accident can be reproduced by EGS5 calculation. The calculated (C) and measured (M) ambient dose rates agreed reasonably well; the ratio C/M was 0.62–1.28.

Acknowledgments

We express our special thanks to Dr. M. Hagiwara, Mr. S. Takahara, and Mr. K. Hozumi who helped us in measuring the ambient dose rates on the KEK campus.