1- Introduction

Primordial radionuclides have been present in the crust of the Earth since its formation. The natural decay-chains of the uranium and thorium series and of ⁴⁰K radionuclide are usually accompanied by the emission of gamma ray radiation gathering information about the primary source. The release of anthropogenic radionuclides from atmospheric testing of nuclear explosive devices and releases from nuclear reactors presents another source of environmental gamma ray radiation.

Depending on gamma ray interaction mechanisms with different materials, many gamma radiation detectors have been conceived and studied [1, 2]. These detectors are used for several purposes, ranging from astrophysical observations [3] to gamma ray field measurement in radiation processing and radioprotection [4-9]. In the present work, we are concerned specifically with two types of detectors; the semiconductor-based detectors and the scintillation detectors. Their principles are described in ref. [10] and [6]. An advantage of such detectors is the possibility of manufacturing light and portable prototypes. In the present work, a semi-conductor based portable detector has been studied, but the same method could be applied to a scintillation portable detector. Portable detectors are widely used for in-situ measurements to collect information about the source distribution in the ground. The response of a portable detector to incoming photons is measured as a pulse height distribution, which essentially represents the spectral distribution of energy depositions in its active volume. Full energy depositions of primary photons lead to full energy peaks. But a substantial fraction of the recorded events corresponds to photons for which only a partial deposition of their incoming energy took place, due to the escape of Compton-scattered photons and of secondary radiation (fluorescence photons, annihilation gamma-rays, bremsstrahlung and electrons) from the active volume, as well as to energy loss after scattering from the non-active parts (housing and cryostat) of the detector.

In conventional applications of in-situ gamma-ray spectrometry with portable detectors, only counts of full energy peaks are determined in the measured detector response. The photopeaks coupled with careful and precise energy calibration, are used for the determination of the nuclides in a sample. The net area of the peak is directly related to the activity of each isotope after correcting for the efficiency of the detector, the branching ratio of the source and its half life. But it is of considerable interest to obtain from the detector response the full spectral distribution of the incident photon fluence, including the continuum component coming from photons that have scattered in the environment. This component contains valuable information about the source distribution [11]. For example, given the energy dependence of the attenuation of gamma rays in the soil and the ratio of the scattered radiation to the unscattered one, it is possible to extract information about the depth distribution of the source [11].

In-situ Gamma ray field measurement requires a prior knowledge of the detector response. However, the large scale geometry of the soil and air media involved and the reduced dimensions of the detector make it difficult to perform a field calibration of the portable detectors. To deal with this problem, a semi-empirical method has been widely applied in previous works [6-9, 12, 13]. In the semi-empirical method, the response of a portable gamma detector is described by a calibration factor. The calibration factor for in-situ measurements can be expanded into a term corresponding to the intrinsic efficiency, which deals with the partial absorption of the energy of the incident events, and a term corresponding to the angular response of the detector. According to the literature, these two terms are obtained experimentally in the laboratory by studying the response of the portable detector to the gamma ray flux emitted by point-like radioactive sources of known activity and known emission energy distant from the portable detector and positioned at different angles. In certain works, the same experience has been simulated by numerical techniques [9]. This protocol attempts to mimic the semi-infinite geometry which the detector will be exposed to in real measurement conditions. Details of this method can be found in ref. [10] and in ref. [13].

Despite being effective, the semi-empirical method does not take into account the real distribution of the gamma ray field in the air emitted by sources in the ground. This represents a source of uncertainty, especially when one tries to identify sources of very low activity. Nowadays, Monte Carlo methods are widely used to solve complex physical and mathematical problems. In Monte Carlo simulation of radiation transport, the history of a particle is viewed as a random sequence of free flights that end with an interactive event. The Monte Carlo model is viewed as a good reproduction of reality. Better description of the distribution of the gamma ray field could be achieved by applying a convenient Monte Carlo simulation. Thereby the response of a portable detector could be determined as accurately as possible, avoiding the uncertainty hat arises from the application of the semi-empirical method. It is worth to mention that in the literature, it was not possible for authors to use only Monte Carlo simulation for the determination of the entire response of portable detectors for an open field geometry such as in-situ measurements. Rather Monte Carlo technique was a support for a semi-empirical protocol of calibration.

In the present work, a Monte Carlo based approach is proposed to determine the portable gamma detector response for in-situ measurement. This approach is based on the assumption that each photon track reaching the detector should originate from an event starting at a spherical surface enclosing this detector. We denote this surface by the detector phase space shell. The approach consists in two-stage Monte Carlo simulation. In the first stage, an optimized geometry of the soil-air medium is applied to determine the probability distribution of the phase space coordinates of the events that will originate from the detector phase space shell. In the second stage, events starting from the phase space shell are propagated and the detector response to those events is studied.The full energy absorption spectrum and the partial energy spectrum are obtained. This approach permits to generate events in the detector phase space shell that contribute effectively to the gamma ray field distribution around the detector volume. It also permits to reduce considerably the amount of computation time needed to calculate the portable detector response in the open soil-air geometry. For validation, this method is applied to the Canbera p type HpGe detector used in ref. [13]. The simulated response is compared to the one obtained in the previous study based on semi-empirical method. As an application of the method, the effective depth of ¹³⁷Cs artificial radionuclide is calculated in the same location where the experimental spectrum is measured.

2- Material and Methods

2-1 Method description

The transport of a gamma photon that starts from the source and propagates in the soil-air medium consists of a dynamic evolution in a 7-dimentional phase space where the position of the photon is described by the three coordinates x, y and z, its direction is described by the three components of its momentum direction p_x, p_y and p_z and its energy by one variable. In this description, the detector is represented by a geometric region in space. This region serves at scoring the useful events.

For the case of in-situ gamma ray field measurement, the dimensions of the detection region are typically small compared to the whole geometry of propagation. Monte Carlo methods are based in determining the detector response by simulating all relevant physical processes taking place along the path of a photon emitted by the source. The main disadvantage of this technique is that a large number of histories must be simulated to obtain a low statistical uncertainty. Therefore, Monte Carlo methods need long computing times. To circumvent this drawback, the detector response should be obtained indirectly.

The basic idea behind the actual work is that every event which corresponds to a photon track reaching the detector region is considered to be starting at a closed surface surrounding the detector. Let’s call this surface the detector phase space shell. The detector response to gamma ray photons emitted by radionuclides distributed in the ground could be obtained by tracking only the photons that emerge from the shell events. Every event is equivalent to a set of seven coordinates sampled from predetermined probability distributions.

It is possible to reconstruct the probability distribution of the phase space coordinates at the shell level by tracking every photon emitted by a source in the ground which has a high probability to reach the center of the detector region. This is achieved by simulating the soil-air medium by an optimized geometry allowing for the tracking of the photons most likely to be detected.

The optimization of the soil-air geometry is based on the fact that the probability for a radiation to reach the detector is strongly correlated to the thickness of matter separating the emitter from the detector. The contribution of the cross-sections of the physical interactions encountered by the photons is ensured throw the linear attenuation coefficient. Since the gamma ray is strongly attenuated in the soil within a nominal depth d equivalent to few tens of centimetres of soil thickness, a probability cut-off could be defined for an emitted radiation to reach the detector. If the detector is at distance h from the ground level, the attenuation factor for a reference point source at the nominal depth d on the vertical axe containing the detector is ${\mu }_{s}d+{\mu }_{a}h$. ${\mu }_s$ and ${\mu }_a$ are the linear attenuation coefficients in soil and air respectively. Thus, for an emitted radiation by a point source within the earth distant of $x_s+x_a$ from the detector the total attenuation factor verifies

${\mu }_s{x}_s+{\mu }_a{x}_a\le {\mu }_{s}d+{\mu }_{a}h$ (1)

It was demonstrated in ref. [8] that the set of the useful sources emitting this radiation are contained in a soil volume bordered by the surface defined by the equation

$r(z)={\left({(\widehat{\mu }h+d)}^2{\left(\frac{h-z}{\widehat{\mu }h-z}\right)}^2-{(h-z)}^2\right)}^{\frac{1}{2}}$ (2)

where r represents the radius at soil depth z, $\widehat{\mu }=\raisebox{1ex}{${\mu }_a$}\!\left/ \!\raisebox{-1ex}{${\mu }_s$}\right.$ is the ratio of linear attenuation coefficients. Following the same arguments, the air is supposed to be composed of virtual sources reemitting scattered radiation to the detector centre. The emitted radiation beyond a reference point in air distant of $x_0$ from the detector is supposed to have a probability to reach the detector lower than a given cut-off. This cut-off corresponds to an attenuation factor ${\mu }_a{x}_0$ in the air material considered. The corresponding air geometry is a spherical volume centred on the detector and of radius $x_0$. For convenience we choose ${\text{x}}_{\text{0}}=\sqrt{r^2(z=0)+h^2}$. This reference point corresponds to the furthest real source from the detector which lies on the ground. In this way, the suitable soil-air geometry depends on the nominal depth of soil as sole parameter. Detailed properties of the optimized geometry and its application to the gamma ray field measurement are found in [14] and [15].

Based on this methodology, the simulation of the response of a portable gamma detector is fulfilled throw two independent stages. In the first stage, the probability distributions of the phase space coordinates of the events at the shell level are reconstructed. In this stage, the soil-air medium is simulated by the optimized geometry. It is worth to mention that this stage is performed only once. The resulting distributions could be used to simulate the response of any kind of portable detector. In the second stage, the response of the studied portable detector to events generated according to the distributions reconstructed in the first stage is determined. A diagram explaining the methodology is shown in Fig. 1.

Fig. 1.Full-screen View

A sketch of the geometrical layout showing necessary parameters to calibrate a portable detector for in situ measurement by Monte Carlo simulation.

For computer simulation, we use the Geant4 toolkit. Geant4 is a Monte Carlo transport code that is based on object-oriented technology. It was developed using the C++ computer language and provides transparency of physical implementation [16, 17, 18]. The low-energy electromagnetic physics package of Geant4 was used in this simulation. This package includes improved models for the simulation of photon, electron, hadron and ion interactions. It also includes models for the simulation of atomic relaxation and other atomic shell effects as well as unique models describing the interactions of radiation with biological systems [19]. The package is capable of describing interactions at energies down to 100 eV and can perform very detailed simulations of particle transport in a medium. Geant4 geometry package allows a set of geometric shapes along with some Boolean operations, such as intersection and union, to build more complex geometries.

2-2 Probability distributions of phase space coordinates of the events at the detector phase space shell level

The soil-air medium was simulated by an optimized geometry of nominal depth equal to 50 cm. This choice is adequate since 90% of gamma radiation with the maximal energy of 3 MeV is attenuated within this soil thickness. The air radius is thus equal to 622 m. Soil and air compositions are shown in Table 1.

Table 1:Full-screen View

Compositions of soil and air in weight fraction considered in the simulation.

Elements	Soil (density=1.6 g/cm3) Weight fraction	Air(density=1.29 mg/cm3) Weight fraction
H	2.2	-
O	57.5	23.2
Al	8.5	-
Si	26.2	-
Fe	5.6	-
N	-	75.5
Ar	-	1.3

Show more

The detector phase space shell consists of a sphere of 80 cm radius, centered on the detector. This highly symmetric closed shape would enable an unbiased simulation. Four photon-interactions processes — photoelectric absorption, Compton scattering, Rayleigh scattering and pair production — were considered. The simulation was performed for primary energies corresponding to the natural radioactive series of ²³⁸U, ²³²Th and the ⁴⁰K radionuclide. In total 49 energies of the natural radioactive spectrum were simulated. The data was taken from [20]. These energies correspond to the most intense radionuclides.

A total number of 10⁸ of photons were generated uniformly in the soil per primary energy. For each inward crossing of the shell sphere by an incoming photon, its energy, position and momentum direction were recorded. The resulting distributions of the seven coordinates were normalized to the total number of records. In this manner we obtained the probability distribution of each coordinate. Some of the normalized distributions are shown in Fig. 2.

Fig. 2.Full-screen View

Probability distributions of the seven phase space coordinates at the phase space shell for ⁴⁰K radionuclide uniformly distributed in the ground.

As clearly seen, the distributions of position and momentum coordinates in the XY plane are almost flat. This indicates the rotational symmetry of the problem around the Z axis. It is worth to mention that for a soil volume simulated using the optimized volume; we achieved a variance reduction of two orders of magnitude compared with classical Monte Carlo simulation with a considerable reduction of computation time. An average time per event of 4×10^-5 s on a Core Duo processor was recorded.

3- Simulation of the response of p-type coaxial High Purity Ge (Li) portable detector

We simulated the response of a Canberra p-type coaxial High purity Ge (Li) detector (HpGe) of Φ 58 mm×56 mm with a relative photo-peak efficiency of 30% at 1332 keV. This is the same detector used in [13]. Its layout as supplied by the manufacturer is shown in Fig. 3. It has been shown in [13] that the detector dead layer value should be 1.84 mm instead of 0.84 mm provided by the manufacturer. We toke this detector dead layer value throughout the simulation. 10⁸ of events were generated at the shell level per energy of the natural radioactive spectrum. Each event corresponds to a photon with sampled energy, direction and position from the relevant distributions already prepared in the first stage. In GEANT4, the following electro-magnetic processes were included: Compton scattering, photo-electric effect, Rayleigh effect, pair production, multiple scattering, fluorescence and Auger effect, bremsstrahlung and ionization.

Fig. 3.Full-screen View

The portable HpGe detector geometry used in the present work. The structure is cylindrically symmetric.

The full energy peak efficiency is defined as the ratio of the number of recorded events corresponding to full absorption in the detector crystal of incident primary photons having the energy of the source, to the total number of the incident primary photons on the detector. In gamma ray spectroscopy, the full energy peak efficiency allows transforming total peak counts to concentration of the correspondent radionuclide in the ground. The detector full energy peak efficiency is function of the primary energy.

Fig. 4 shows the full energy peak efficiency for different radionuclides as deduced by our method, compared to the efficiency curve deduced from [13]. The values obtained by our method are in agreement with those calculated in [13]. The observed discrepancy at very low energy could be explained by the angular response of the detector considered to be uniform in the previous work, a hypothesis which is not true for the low energy gamma photons which experience severe scattering in the air medium. In our method, by generating events at the shell level according to the obtained distributions in the first stage, we take into account the real distribution of the gamma ray field in the detector proximity. Furthermore, in ref. [13], the efficiency curve has been extrapolated to higher energies. This could be the source of the observed discrepancy between the two methods at energy values higher than 2000 keV. For the energies between 600 keV and 2000 keV, efficiencies calculated by Monte Carlo method exceed the ones obtained from [13] by less than 5%. It is worth to mention that in ref. [13], the efficiencies were calculated with 5 % accuracy. We conclude that the observed discrepancy is not relevant in the energy range of interest.

Fig. 4.Full-screen View

The Full Energy Peak efficiency of the portable HpGe detector used in the present work as calculated by Monte carlo simulation for different natural radionuclides (markers) compared to the efficiency curve (line) deduced by semi-empirical method from ref. [13].

In Fig. 5, the spectra of partial absorption and full absorption are shown for different radionuclides of the natural decay series. It can be noted that a large amount of the partial absorption events occurred at the low energy part of the spectra. We also notice the appearance of some secondary peaks at this low part. These peaks are due to the escape and to the creation of X-ray and Auger electrons in the detector housing materials. Since the most important radioactive peaks are located above 20 keV, a cut-off applied at this energy permits to cancel the underlying undesirable background.

Fig. 5.Full-screen View

Full energy absorption (dark color) and partial energy absorption (clear color) spectra obtained by Monte Carlo simulation for the portable HpGe detector used in the present work. Simulations were conducted for different natural radionuclides uniformly distributed in the ground.

For radionuclides with higher energies, the partial absorption events became more and more frequent due to the escape of high energy photons. This photon escape is due to the small volume of the Germanium crystal. Thus, the detector used in this work is more sensitive to radionuclides emitting low energy gamma rays. In a real detector, if the incident gamma energy E_γ is above 1022 keV, pair production events result in the production of two 511 keV annihilation gamma-rays. If only one of these gamma-rays escapes while the other is completely absorbed in the detector, a separate peak in the spectrum appears at the energy E_γ-511 keV. This peak is called the single escape peak. If both annihilation gamma-rays escape this gives rise to the double escape peak at E_γ -1022 keV. The double escape peaks are shown in Fig. 5 for ⁴⁰K and ²⁰⁸ Tl radionuclide's spectra. These peaks are very pronounced because pair production is the most relevant interaction of high energy gamma photons with the detector crystal. Furthermore, in the ²⁰⁸ Tl spectrum, there is a second peak composed of full absorption events of 511 keV annihilation photons. The intensity of the single escape peak and double escape peak depends usually on the size of the Ge crystal. For a small detector compared with the mean free path of the secondary gamma radiation, all annihilation photons escape from the detector. The spectrum will show a double escape peak but no single escape peak. This could explain the absence of single escape peak in ⁴⁰K and ²⁰⁸ Tl radionuclide's spectra. The pronounced double escape peaks of high energy spectra represent a significant background proving the fact that the actual detector is suitable for identifying low energy gamma sources.

Fig. 6 shows the gamma ray simulated spectrum for the natural decay chain series after removing the partial absorption spectrum. This spectrum is normalized to its experimentally observed full-energy peak. In the same figure, the gamma ray experimental spectrum recorded in-situ is shown. A matrix stripping method explained in detail in [13] was applied to the experimental spectrum to remove the partial absorption events. The shapes of the two spectra are close. The observed discrepancy could be attributable to the accuracy of the stripping method described in [13]. The discrepancy could also be attributable to the background component due to the electronics and contamination of the detector material which is not taken into account in our model. This background component is dominant at the high energy part of the spectrum where the detector sensitivity to the emitted energies is degraded. As a result, the discrepancy at this part of the spectrum is the most relevant.

Fig. 6.Full-screen View

Total simulated spectrum compared to the experimental spectrum for the portable HpGe detector used in the present work. The markers indicate the experimental peaks identified by the peak finder utility of the Root toolkit system [21].

The discrepancy could also be attributable to how accurate is the assumption that the different radionuclides are homogeneously distributed in the soil and possibly to the different chemical composition and moisture of the soil where the measure was taken from the data used in the simulation (Table 1).

4- Application

As a consequence of nuclear probes and nuclear accidents a certain amount of the fission product ¹³⁷Cs was ejected into the atmosphere, from where it was deposited on the soil. ¹³⁷Cs is a gamma emitter with energy of 661.6 keV and with half-life of 30.17 years. Several models have been developed to describe ¹³⁷Cs migration in soils and to explain its vertical distribution [22, 23, 24, 25]. After the nuclear tests in the sixties and the Chernobyl accident it was initially possible to describe the depth profile of the ¹³⁷Cs isotope in the soil by an exponential distribution. It was also shown [25] that for any exponential distribution of ¹³⁷Cs in the soil, an equivalent planar distribution of ¹³⁷Cs radionuclides can be found at an effective depth, which yields the same photon fluence energy and angular distribution as the original one.

In the present work, the effective depth was determined in tow independent steps.

In the first step, we applied a cut-off at the bin energy containing the ¹³⁷Cs peak to the experimental and to the simulated in-situ spectra shown in Fig.6. We conserved only the low-energy part of the two spectra. It is worth to mention that the simulated spectrum is due to natural radionuclides uniformly distributed in the ground. The simulated spectrum was then stripped from the experimental one, leaving one with the residual spectrum belonging to the ¹³⁷Cs isotope only. This spectrum was utilized to obtain information on the ¹³⁷Cs depth distribution.

In the second step, we assumed a plane distribution of ¹³⁷Cs activity at a certain depth beneath the surface. For each location of the plane source, a Monte Carlo simulation was conducted. According to the optimization applied to the soil geometry, any plane source is reduced to a disc of radius r(z) given by Eq. (2). The effective depth was determined by comparing the simulated ¹³⁷Cs spectrum due to the plane source and normalized to the experimental full-energy peak value, with the residual experimental spectrum. The location of the ¹³⁷Cs plane beneath the soil surface is varied until a good match between the two spectra is obtained.

The optimal effective depth was found equal to 40 cm as inferred from Fig. 7. For an effective depth of 40 cm, the planar source is reduced to a disc of only 145 cm radius. It is worth to mention that the experimental in-situ measurement was taken in the soil region located in the vicinity of the site studied in ref. [26]. In this experimental study, the authors found a maximal activity of ¹³⁷Cs located between 30 cm and 40 cm which is in accordance with our prediction.

Fig. 7.Full-screen View

The Simulated spectrum for a disc source of ¹³⁷Cs radionuclide at a depth of 40 cm (full line) compared to the residual spectrum due to the same radionuclide (dash line).

5- Conclusion

A Monte Carlo based approach to study the response of a portable gamma detector to terrestrial gamma ray was proposed. This approach consists in two-stage Monte Carlo simulation. In the first stage, the probability distributions of the phase space coordinates at the detector phase space shell were reconstructed for natural radionuclides uniformly distributed in the ground. For this purpose, an optimized geometry of the soil-air medium was implemented to ensure the tracking of only the events most probable to be detected. These distributions serve at studying the response of any type of portable detector for in-situ measurement. The response to events originating at the phase space shell level according to the already prepared distributions of a p type HpGe detector was determined. This method allowed to separate the background signal, consisting of the partial absorption energy spectrum, from the desired signal, consisting of the full absorption energy spectrum, for any natural radionuclide uniformly distributed in the ground. Comparison with a previous study which uses a model based on semi-empirical method has shown a good agreement for the full energy peak efficiency curve and the total in-situ spectrum. As an application of this method, the effective depth of ¹³⁷Cs artificial radionuclide in the ground at the location where the in-situ spectrum was carried out was determined with good precision. This approach is among the first attempts to determine the entire response of portable HpGe detectors for in-situ measurement of gamma ray field by Monte Carlo simulation. We believe that we could improve this approach by taking into account other parameters such as soil and air real properties, the detector resolution [27] and the detector electronic background to enhance its predictive power.