Skin dose contamination conversion coefficients. Benchmark with three simulation codes

Handling of radioactive material by operators can lead to contamination at the surface of the skin in case of an accident. The quantification of the dose received by the skin due to a contamination scenario is performed by means of dedicated dose coefficients as it is the case for other radiation protection dose quantities described in the literature. However, most available coefficients do not match realistic scenarios according to state-of-the-art of science and technology. Therefore, this work deals with dedicated dose conversion factors for skin contamination. Since there is an increasing demand on dose coefficients in general, these specific coefficients can be used for various calculations in radiation protection. In this work a method to evaluate such coefficients for the skin contamination dose related to photons, electrons, positrons, alpha and neutron particles is proposed. The coefficients are generated using Monte-Carlo simulations with three well established calculation codes (FLUKA, MCNP, and GEANT4). The results of the various codes are compared against each other for benchmarking purposes. The new dose coefficients allow the computation of the skin received dose, in the case of skin contamination scenario of an individual, taking into account the decay radiation of the radionuclides of interest. To benchmark the quantity derived here, comparisons of radionuclide contamination doses to the skin using the VARSKIN code available in the literature are performed with the results of this work.


Introduction
Handling of radioactive material by operators can lead to radioactive contamination at the surface of the skin in case of an accident which goes beyond the risk of exposure. When skin is contaminated, it can be required to quantify the dose received by an individual due to a known surface-related activity of radionuclides. Hence, applicable dose conversion factors are required. Since there is an increasing demand for dose coefficients in general, specific coefficients for skin contamination doses, as described in this work, can be used for various radiation protection calculations [1,2] such as for an accident scenario or for the determination of activity limits, e.g. in International Atomic Energy Agency (IAEA) regulations [3][4][5].
The majority of the available data in the literature for the dose contamination quantity have been created more than 25 years ago before the development era of modern calculation tools such as intensive Monte-Carlo calculations. Therefore, the current derivation of skin contamination dose data relies on calculations that were performed with outdated tools and data available that have evolved consequently for more accuracy today. Moreover, only beta skin contamination dose tabulated data is available in the literature. For instance, Cross et al [6,7] produced beta dose data for skin contamination considering a water volume phantom. In [8], the authors computed beta dose rate for skin contamination for 21 nuclides at different depth.
Some tools allow the assessment of skin contamination dose, due to electrons, positron and gamma particles, for various radionuclides. The code VARSKIN [9], initiated in 1987, has evolved among the years to maintain accurate and efficient physics models. It is widely used to compute the skin dose contamination originating from beta and gamma radiation, for specific nuclides with known activities. Also, Covens et al [10] computed absorbed skin dose rate conversion factors for 39 radiopharmaceutical nuclides of interest.
In this paper, a method to construct data for skin contamination dose is proposed. Various particles from radionuclides decay channels are considered: electron, positron, photon, alpha and neutron particles. The data produced are cross compared among three different Monte-Carlo calculation codes that are extensively used worldwide: FLUKA [11,12], GEANT4 [13][14][15] and MCNP [16,17]. The geometry model is derived from the protection quantity standard described in the ICRP116 [18] publication.
Section 2 focuses on the geometry model and a description of the Monte-Carlo codes used to derive the skin contamination dose coefficients (SCDC). Then, section 3 details the results and focuses on the discrepancies obtained with the three Monte-Carlo codes. A numerical validation is finally presented in section 4 for a set of radionuclides by comparison with VARSKIN.

Calculation models
The geometry model used to assess the SCDC is based on the one provided by ICRP116 [18] for local skin dose except for two differences. A cube of 10 × 10 × 10 cm 3 of skin material is located in air whereas the ICRP116 model for local skin dose, the skin cube is in vacuum (see figure 1). The skin material is described in table 1. The scoring volume, in which the deposited energy is scored, is a cylinder, located between 50 and 100 µm depth in the skin cube and with a surface of 1 cm 2 . The source is a disk of 3.5 cm radius. The source is uniformly distributed within the disk and emits isotropically at contact with the skin cube, contrary to the ICRP116 setup in which the source is a beam normally impinging on the skin cube. The two differences with ICRP116, detailed above, originate from the SCDC quantity definition created here. While ICRP116 provides the geometry to build Local Skin Dose Conversion Coefficients, skin dose coefficients from contamination are constructed in this work. Hence, a particle source that is in contact with the skin and emits isotropically is simulated. In summary, this approach is optimized for the calculation of SCDC.
For each primary beam type, the total averaged energy deposited over the scoring volume by all particle types is obtained in GeV/g or MeV/g per incident particle. This quantity is converted into an absorbed dose in pGy. Finally, the result is multiplied by the surface of the source in contact with the skin, to get the SDSC conversion coefficient in pGy.cm 2 . Three Monte-Carlo calculation codes are used to simulate five different primary particles of interest that are produced during the decay emissions of radionuclides. The calculation codes are FLUKA 2011.3 [11,12], MCNP version 6.2 [16,17] and GEANT4.10.06 (or higher) toolkit [13][14][15]. The primary particles of interest are photons, electrons, positrons, as well as alpha and neutron particles.

FLUKA simulations
FLUKA Monte-Carlo code is widely used for radiation protection purposes in the particle physics and accelerators community. The code is based on original physics models that cannot be changed by the standard user. Therefore, the user can only change precision parameters such as energy thresholds and activation/deactivation of physics processes.
FLUKA is used in this work for the simulation of photons, electrons, positrons, alphas and neutrons as primary source particle. All the secondary particles created during the transport process are also taken into account.
For photons, electrons and positrons, transport and production energy thresholds are lowered to the minimum possible value of 1 keV. Electron and positron interactions are based on different stopping power and bremsstrahlung data. Positron annihilation in flight or at rest is simulated within the FLUKA code. Continuous and discrete energy loss depending on particle energy is modelled exactly. For photons, pair production, Compton and photoelectric effects are modelled. All secondary particles are transported.
For neutrons with energy lower than 20 MeV, only multigroup formalism is available in FLUKA with 260 groups discretization. The cross-section data are derived from different evaluation databases that are available in the literature and the most commonly used: ENDF/BVII, ENDF/BVI [19] and JEFF-3.1 [20]. The results obtained with FLUKA are subsequently rebinned considering the ICRP107 [21] neutron binning. The energy deposition is assessed in FLUKA by means of kerma approximation factors (see FLUKA's manual 4 ). A detailed description of kerma factors impact on deposited energy is available in [22]. More precisely, among  • Recoil protons; • Protons from (n, p) reactions.
All secondary particles originating from neutron primaries and contributing to the dose are taken into account: electrons, positrons, photons, heavy and light ions, alpha and protons. Secondary photons from neutron interactions are generated according to a multigroup treatment as well.
The thermal scattering correction is activated for H bound with water, for the skin material.

GEANT4 simulations
GEANT4 toolkit is a fully integrated Monte-Carlo code for interaction and transport of particles developed by the GEANT4 collaboration at CERN. In this work GEANT4 is used for the simulation of photons, electrons, positrons, alpha and neutron primary particles. All the secondary particles created during the transport process are also taken into account. Within the Geant4 toolkit many reference physics lists are provided to investigate specific issues. For this work the QGSP_BIC_HP physics list is used which includes the electromagnetic physics option 4. This uses the most accurate standard and low-energy electromagnetic physics models. With its concentration on the best possible physics, electromagnetic option 4 is slower than the standard EM package but best suited for the intended use of this work. Stopping powers for electrons and positrons are different. All secondary particles are transported.
For neutrons, the mentioned physics list adopts the high-precision (HP) neutron package to describe neutron interactions from thermal energies to 20 MeV. The G4NDL4.6 nuclear databases is used. This allows for considering point wise cross sections. The thermal scattering correction is activated for H bound with water, for the skin material. This choice of physics for neutron simulation is driven by [23] in which the authors provide recommendations by comparisons with data from the literature. Simulations are done with the same energy binning than FLUKA for comparison purpose. Then, results are rebinned for compatibility with ICRP107 [21] binning. No energy threshold is applied. The two patches regarding exit direction of hadronic elastic collisions 5 and inelastic excitation energy at the exit of inelastic collisions 6 are applied.

MCNP simulations
MCNP is the simulation code developed by Los Alamos National Laboratory. Like the previous codes described in this paper, it is a general-purpose Monte-Carlo code designed to track many different particles. In this work, it is used to simulate the transport of photons, electrons, positrons, alphas and neutrons as primary source particle. All the secondary particles created during the transport process are also taken into account.
MCNP version 6.2 [17,24] is used. The considered nuclear libraries are derived from ENDF/B-VII [19] for neutrons and protons and from ENDF/B-VI.8 for photons. The electron transport algorithms and nuclear libraries are derived from the ITS code version 3 [25]. Positron physics in MCNP6 is identical to electron physics, except for tracking directions in magnetic fields and consideration of positron annihilation. Whereas electrons below the energy cutoff are terminated, positrons below the energy cutoff produce annihilation photons. All secondary particles are transported.
The thermal scattering correction in light water is taken into account for the skin material as recommended in [26] for the neutron's simulations.
For all particles apart from neutrons, the transport energy threshold is lowered to 1 keV. The tally '+ F6' is used to record the deposited energy and relies on kerma approximation factor. This tally has no particle designation and automatically applies to all charged particles.

Particularity of alpha models
For simulations of alpha primary particles with GEANT4 and FLUKA, the standard physics model is considered. For MCNP [27], the Vavilov model is used. Tests are performed with MCNP turning off the Vavilov model and using instead the Continuous Slowing Down Approximation and with FLUKA turning on the IONFLUCT to activate ionization fluctuations. Slight differences were observed when changing these models. However, it has to be pointed out that alpha simulations rely on stopping powers that are computed using the Bethe theory. However, the analytical formulation of Bethe theory is not valid at our energies of interest. Therefore, at low energies, corrections are done to complement the physics model. As described in [28], FLUKA, MCNP and GEANT4 have all unique ways of computing the stopping powers, which could lead to differences in alpha transport results.

Radiation weighting factors
The absorbed dose is a measurable quantity. ICRP103 [29] defines radiation weighting factors to account for the effectiveness of the radiation type and express the equivalent dose. In the next, our results are expressed in picogray (pGy) for the absorbed dose and picoSievert (pSv) for the equivalent dose.
For photon, electron and positron primary particles, the radiation weighting factor (W R ) value is 1. Hence, the results for such primary particles can be expressed interchangeably as absorbed or equivalent dose.
The considered radiation weighting factors (w R ) for neutrons are those listed in [29]. W R is applied at the energy of incident neutrons before the interaction with the skin 7 and is described as a continuous function of the incident neutron energy. As in our simulations there is only air around the neutron source and the skin, we can just assume that the neutrons enter the skin with the neutron energy of your source. The equation below summarizes the values of W R for incident neutron at energy E in MeV: Finally, the radiation weighting factor (W R ) for alpha particles is set to 20, as suggested in [29].

Results
The simulations are performed for different primary particle as a function of energy. Results are presented in tables 3-7 (appendix). Corresponding figures 2-4 illustrate the results and ratios obtained between the different Monte-Carlo calculation codes.  Table 3 provides the SCDC for primary photons, electrons and positrons. The energy ranges of the primary particles vary from 1 keV to 10 MeV and comparisons are provided between FLUKA, GEANT4 and MCNP. In figure 2, the corresponding curves are shown and the ratio of each code with the average of the three code results is represented by crosses. The statistical uncertainties are below 2% at 1 σ for all calculated points apart from one exception with primary electrons as discussed below.
For photons, the three curves are almost in total agreement as the ratios are comprised between 0.95 and 1.05, i.e. in the range of the statistical uncertainties. The plot in figure 2 shows an increase from 1 keV up to 4 keV which is the photon energy necessary to record events in the scoring volume. As the energy increases, photons escape the scoring volume explaining the decrease until around 100 keV. Above 100 keV, the Compton effect outside the scoring volume causes photons to generate secondary particles depositing energy in the scoring volume.
For electrons, at low energies the ratios between each code and the average results reach a value of 1.10 and 0.9 for the lowest primary energies of 3 and 4 keV, respectively. The probability that a radiation reaches  the scoring volume from electrons emitted at those energies is low, hence, the convergence of the calculation is complex and requires larger calculation time. For those two energies, the statistical uncertainty is 50% and 20%, respectively, explaining the large discrepancy between the two calculation codes. Doses from electrons below 55 keV are caused by bremsstrahlung and fluorescence and only secondary photons reach the scoring volume. For the electrons, the curves present a rise at around 55 keV primary energy. This is explained by the electron energy needed to reach the scoring volume. Above that energy, they start reaching the scoring volume and deposit part of their energy through collisions. As the electron energy increases, we observe a flat behavior from 100 keV to 10 MeV in figure 2. This can be explained by the isotropic electron emission from the source. At high energies, electrons can deposit energy entering the scoring volume from the sides as they have larger track length compared to the minimum distance between the source and the scoring volume. From 5 to 50 keV the statistical uncertainties are below 10% also explaining the discrepancies between the calculation codes. Finally, for electrons emitted between 50 and 90 keV, the ratios vary from 0.95-1.05. This difference cannot be explained by the statistical uncertainties which are for these energy points below 2%. However, from 50 to 90 keV, there exists a transition in the particle contributions to the dose. At this transition, electrons start reaching the scoring volume while, at low energy, only secondary photons from electron energy conversion are hitting the scoring volume. This effect depends on the stopping power that vary between FLUKA, GEANT4 and MCNP. Outside of these energy ranges, the SCDC are similar between the three calculation codes with a ratio in the range of the statistical uncertainty.
Finally, for positrons, the agreement between the three codes is within the range of the statistical uncertainty except for primary positrons emitted between 8 and 120 keV in which discrepancies are observed between MCNP, FLUKA and GEANT4. From 8 to about 80 keV, the dose is mainly deposited by photons from positron-electron annihilation. The differences in the physics models (annihilation) between the codes could be responsible of these discrepancies. Regarding the range from 80 to 120 keV, the discrepancies could be explained by the same reason as for electrons. In this range, a transition is happening and positrons are starting to reach the scoring volume, explaining the dose increase. The effects of stopping power are very important in this region and all codes compute different stopping powers. As for the electrons, positron curves present a rise at around 55 keV primary energy. This is explained by the positron energy needed to reach the scoring volume. Above that energy, they start reaching the scoring volume and deposit part of their energy through collisions. As the positron energy increases, we observe a flat behavior from 100 keV to 10 MeV in figure 2 with a slight slope. This can be explained by the isotropic positron emission from the source. At high energies, positrons can deposit energy entering the scoring volume from the sides as they have larger track length compared to the minimum distance between the source and the scoring volume. Figure 3, tables 4 and 5 provide the SCDC for primary neutrons emitted between 4.14.10 −07 -15 MeV. The statistical uncertainty is below 2% at 1σ. The results are given with the multigroup formalism corresponding to the ICRP107 binning of neutron spectra [21]. This choice has been oriented to convolute SCDC with ICRP107 decay data and directly obtain the skin contamination dose for a radionuclide of interest. The agreement between the three calculation codes is within the range of the statistical uncertainties from 0.1 to 15 MeV. Slightly larger discrepancies are observed below 1 keV where GEANT4 seems to underestimate the SCDC values compared to the other codes. Even though the discrepancies are negligible, the differences in the physics can be pointed out: • The neutron nuclear databases used by the simulation codes. GEANT4 simulations rely on G4NDL4.6 which is mainly extracted from JEFF-3.3. MCNP calculations are done using ENDF/BVII for neutrons and END-F/BVI.8 for photons. FLUKA neutron nuclear data takes into account different evaluation files but is mainly base on ENDF/BVII which could explain why FLUKA and MCNP results have a better agreement than with GEANT4 results. • The physics models for the transport of charged particles are different for the three codes, as explained in section 2. FLUKA and MCNP use kerma factors to estimate energy deposition by neutrons as it is detailed in [22]. On the other hand, GEANT4 considers a real neutron transport physics model. Figure 4, tables 6 and 7 below present the results obtained with alpha primary particles and comparisons between the three Monte-Carlo calculation codes FLUKA, MCNP and GEANT4. For primary alpha particles below 5 MeV, no dose was recorded in the scoring volume due to the high stopping power related to alpha particles. The statistical uncertainties are below 3% at 1 σ for all calculated data points apart from 5 to 6 MeV where the uncertainties are respectively below 20%, below 12% and below 7% for 5, 5.5 and 6 MeV primary alpha, due to the high calculation time required since alpha's probability to go through 50 µm skin is low at this energy. The results are in agreement between the three calculation codes in the range from 9 to 20 MeV. However, from 5 to 9 MeV, large discrepancies are observed in particular at 5 and 6 MeV where the ratio can reach two order of magnitudes especially by comparing MCNP with the two other codes. These discrepancies could be explained by the differences in the stopping power used by the three calculation codes, as detailed for alpha particles in [28]. In figure 4, we notice that Geant4 presents an SCDC rise between 6.5 and 7 MeV of primary alpha energies while FLUKA and MCNP present the rise between 6 and 6.5 MeV of primary alpha energies. This feature could be an indication of higher stopping power values for GEANT4 at this energy range. At 6.5 MeV, the dose ratio between FLUKA and GEANT4 is 280 and the dose ratio between FLUKA and MCNP is 400 indicating a lower stopping power for FLUKA.

Discussions on comparisons with values available in the literature
In order to ensure the correctness of our calculation model, comparisons with VARSKIN6.1 [9] are performed. In VARSKIN, a skin thickness of 7.5 µm is considered with a disk source area of 38.5 cm 2 (diameter of 7 cm). ICRP-107 decay database is used [21]. The VARSKIN calculations are compared to calculations based on the SCDC produced in this work with the three codes. Skin Contamination Dose (SCD) for five radionuclides are provided using the following equation with the SCDC produced in this work. where SCDC (E) is the value of the SCDC computed in this work at the energy E; θ rn (E) is the continuous spectrum of the radiation (for neutrons and beta spectra) for the radionuclide rn of interest and I rn (E) is the discrete intensity emission for the radionuclide rn of interest. θ rn (E) and I rn (E) are extracted from ICRP107 [21]. 'Discrete radiations' refers to particle emissions with defined energies such as alpha, photons, Auger and Internal Conversion Electrons. We refer to the Auger and Internal Conversion Electrons as Discrete Electrons.
The results show very strong agreement (see table 2) between the three calculation codes despite the differences pointed out in the stopping power. In particular, we observe the highest discrepancy from MCNP to the average of the three codes for electron dose in the case of Sr90 (ratio value to the average of 0.98). The results are in qualitative and quantitative agreement with regard to the statistical uncertainty between the three calculation codes.
Concerning the comparison with VARSKIN6 (table 2), the results globally show good agreement except for skin contamination dose from photon in the case of Na22 (ratio of 0.5) and skin contamination dose from electrons and photons in the case of U238 (ratio of 10). The difference in the case of Na22 is coming from annihilation photons. In the case of VARSKIN, the dose from annihilation photons is included in the photon dose using the decay data from ICRP107. In the case of this work, the photon dose from annihilation is included in the positron dose as, annihilation photons are simulated in the simulation codes and transported through the scoring volume. Regarding U238 (ratio of 10 with discrete electron dose), a part of the dose is coming from electrons emissions below 10 keV that are not considered in VARSKIN, leading to a lower dose compared to the one obtained in this work. Finally, the discrepancies in the photon dose for U238 between VARSKIN and our codes can also be explained. In VARSKIN6, photons emissions are considered only when the photon yield is above 1% while we considered all the emissions from ICRP107.

Conclusions
New tabulated data for SCDC related to electron, positron, photon, alpha and neutron particles are produced. Those results are compared between three Monte-Carlo simulation codes providing qualitative and quantitative agreement considering the statistical uncertainties. Some discrepancies are observed due to differences in the physics models and nuclear data. In particular, the three codes do not agree for primary alpha emitted below 7 MeV. The main contributor to these discrepancies could be the stopping power for charged particles that can vary consequently from the different calculation codes as detailed in articles from the literature [27,28]. The new values can be used to predict doses from skin contamination for known activities of radionuclide of interest. A validation is performed with VARSKIN6 by comparing electron and photon skin contamination doses for a set of radionuclides. The present work is not in a position to conclude the best Monte Carlo code to use for assessing skin dose contamination without experimental data.
Appendix. Tables of skin contamination dose   Table 3. SCDC in pGy.cm 2 or pSv.cm 2 or pSv.Bq −1 .s −1 .cm 2 of source area for photons, electrons and positrons. Statistical uncertainties are all below 2% except for electrons between 3 and 4 keV where uncertainties are below 50% and electrons between 5 and 50 keV where uncertainties are below 10%.