Update of the Nuclear Criticality Slide Rule for the Emergency Response to a Nuclear Criticality Accident

AWE (UK), IRSN (France), LLNL (USA) and ORNL (USA) began a long term collaboration effort in 2015 to update the nuclear criticality Slide Rule for the emergency response to a nuclear criticality accident. This document, published almost 20 years ago, gives order of magnitude estimates of key parameters, such as number of fissions and doses (neutron and gamma), useful for emergency response teams and public authorities. This paper will present, firstly the motivation and the long term objectives for this update, then the overview of the initial configurations for updated calculations and preliminary results obtained with modern 3D codes.


Introduction
In 1997, Oak Ridge National Laboratory published the report "An Updated Nuclear Criticality Slide Rule" [1][2], as a tool for emergency response to a nuclear criticality accident. A similar document was produced by the Institut de Radioprotection et de Sûreté Nucléaire in 2000 [3].
This kind of document "permits continued updating of information during the evolution of emergency response, including exposure information about accident victims, estimates of potential exposures to emergency response re-entry personnel, estimates of future radiation field magnitudes, and number of fission ("fission yield") estimates, [1]" without precisely knowing the initial conditions leading to the criticality accident.
This document gives order of magnitude estimates of key parameters, useful for emergency response teams and public authorities. In practice, the "Slide Rule" provides estimates of the following information based upon variable times and distances from the accident: • The magnitude of the number of fissions based on personnel or field radiation measurements, • Neutron-and gamma-dose at variable unshielded distances from the accident, • The skyshine component of the dose, • Time-integrated radiation dose estimates, • One-minute gamma radiation dose, and • Dose-reduction factors for variable thicknesses of steel, concrete, and water. The Slide Rule provides estimates for five unreflected spherical uranium systems that give general characteristics of operations typical of facilities licensed by the US Nuclear Regulatory Commission. An example of this Slide Rule is presented in Figure 1. Several laboratories have determined the need to review, update, and expand the contents of the 1997 Slide Rule. In particular, the only conversion factors used to provide doses (Henderson flux-to-dose factors [4]) are outdated, so additional conversion factors will be used. In addition, the results from old radiation transport tools and nuclear data used in originally developing the Slide Rule will be updated with modern tools and nuclear data and some simplifying assumptions will be relaxed.

Long term objectives of this update
As a result, a long term collaboration effort between AWE (UK), IRSN (France), LLNL (USA) and ORNL (USA), performed in the framework of the US DOE Nuclear Criticality Safety Program [5], began in order to update the Slide Rule with modern tools and add new configurations, taking into account the experience of several laboratories in using the 1997 Slide Rule.
The complete work, envisioned to spread over many years, will be divided into five steps: 1. Redo with modern radiation transport tools and nuclear data libraries, for the same configurations and assumptions, the calculations performed initially for the 1997 estimation of the doses; 2. Perform additional calculations to improve the quality/quantity of information given to the user of the Slide Rule in order to not only give a value but also the possible variations and the area of applicability of this value. These additional calculations might include: a. New configurations (impact of the geometry and composition of the source, new fissile media including plutonium systems, multiple layers of shielding, etc.); b. New flux-to-dose conversion factors (for dosimetry, radiological protection and instrumentation purposes); c. Impact of parameters on the result (sensitivity/uncertainty studies, such as thickness and composition of the ground, humidity and density of the air, etc.); 3. Review and improve the section regarding the estimation of the number of fissions; 4. Add other sections to the document like a section regarding actions to stop an on-going criticality accident (for example, standards with neutron poison). 5. Based on the previous work, the final task will be the development of a Slide Rule "application" for a handheld device (e.g. smartphone). At the end, this work should improve the expertise for the real time response to a criticality accident in order to minimize the consequences of such an accident. In addition, this work will provide the opportunity to suggest experiments allowing the complete or partial validation of the tool results (benchmarking effort).
Another consequence of this collaborative effort might be the creation of "computer benchmarks" in order to test and validate the various variance reduction methods and to establish best practices when dealing with this kind of problem.
The present article is focused only on the first step with the presentation of preliminary results for the same configurations and assumptions used in the 1997 Slide Rule.

Critical systems (source)
The critical uranium systems considered for the initial configuration of the 1997 Slide Rule were: Neutron and gamma doses were calculated as a function of distance from the surface of the critical sphere from 0.3048 to 914.4 m (1 to 3000 feet).

Model Description
The geometry for the 1997 Slide Rule models consisted of a simple 2-D air-over-ground configuration with the source located at the radial center of a right-circular cylinder. The radius and the height of the air cylinder was 1530 m. The center of the critical assemblies (spheres) were all 1 m above the ground. The ground is modelled as a 30.48 cm (1 ft) layer of concrete. The dimensions of the critical spheres and the composition of all materials can be found in references [1], [2] and [6]. Figures 2 and 3 present the model for these initial calculations. For more clarity, all the information needed to calculate the initial configuration was written in a specific document using a "benchmark format" [6]. Furthermore, a common file naming convention for the various cases has been adopted. The same principle will be used, in the future, for additional new configurations. An example is the following: • SR-U-UN-G1-C1-d500-N.inp stands for: • SR: Slide Rule, • U: uranium, • UN: unreflected (no shielding), • G1: first case with ground, • C1: first case with uranium system (Uranyl fluoride (4.95%)), • d500: distance of 500 m from the critical sphere, • N: prompt neutron dose.

Initial configuration calculations
The development of the 1997 Slide Rule information utilized the 2-D, discrete-ordinates, coupled neutron gamma-ray radiation transport code, DORT [17]. The source is modelled as a point at the source location with an energy spectrum equivalent to the leakage spectrum generated via a 1-D simulation with XSDRNPM [8]. GRTUNCL [18] was used to calculate the uncollided flux so DORT could use a distributed first collision source in order to mitigate ray effects. The time-dependent sources were generated with the SCALE module ORIGEN-S which uses pointdepletion methods to solve for the radioisotopic compositions after arbitrary irradiation/decay periods. This work assumed an instantaneous event and then tabulates the expected dose rates for various periods after the event for all five critical systems. Specifically, ORIGEN-S first estimates the production of radioisotopes instantaneously during the event, decays the concentrations to the selected times, and then computes the gamma-ray radiation source based on discrete gamma-ray line data. These sources are then utilized in 1-D discrete-ordinates calculations using the same methods, geometry, and codes as those of the prompt dose calculations to determine the leakage spectrum and then the 2-D methods are applied to perform delayed dose calculations.

Results and discussion for the initial configuration
This section presents and discusses the preliminary simulation results for the cases described in the previous section as well as the comparison with the 1997 Slide Rule. As specified in reference [6], without taking into account the skyshine calculations, at least 900 results are needed to cover all the cases. The laboratories used various codes and methods, presented hereafter. Every laboratory used the Henderson flux-to-dose conversion factors [4] in order to compare the results with the 1997 Slide Rule results.

MCNP
MCNP6.1 [7] was used with the ENDF/B-VII.1 cross sections library with a continuous energy representation of the cross sections. An F4 tally (i.e. track length estimate of cell flux) in the shape of a cylindrical shell was used to take advantage of the azimuthal symmetry of the problem. This cylindrical shell was 20 cm tall in the axial direction and 20 cm thick in the radial direction. For prompt doses, the total nu-bar was used. Weight windows, in space and energy, were used for distances beyond 500 m, using the MCNP Weight Windows Generator (with an iterative process, i.e. the weight windows used will provide a result and help creating a new and more effective weight windows that can be used for a new calculation).
For prompt doses, a two steps method was used. The first step is a static calculation (KCODE mode) to determine the distribution of fission neutron production inside the uranium sphere. Many possibilities were considered (and will be discussed later) but all the MCNP results presented were obtained with 20 spherical meshes (SMESH), having the same volume, in which the fission reaction rate were tallied. The second step used the results of the first step to describe a fixed source (SDEF mode) of fission neutrons. A Watt spectrum was used for the energy distribution. The prompt gamma and neutron doses were determined in the same calculation, the gammas being produced by the neutron interactions inside the uranium sphere. Indeed, in the second step, the fission neutron production is turned off (treated as absorption) but the gammas are produced (NONU = 0).
In order to produce the delayed gamma source, the coupling KENO/MONACO/COUPLE/ORIGEN from the SCALE package was used to calculate spectra and intensities. The criticality accident was assumed to occur in 1 µs. KENO calculates the distribution of fission neutron production inside the sphere. MONACO uses the KENO result and calculates the neutron flux inside the sphere. COUPLE uses the resulting MONACO neutron flux and creates problem dependent flux weighted cross sections to produce reaction rates. ORIGEN uses these reaction rates and performs the depletion and decay of the uranium systems. This method was used because it is very close to the method First, KENO was run with a cartesian mesh tally of the fission neutron production, which captured the asymmetry due to the ground 1 m below the center of each fissile assembly. Then MAVRIC/MONACO used the KENO tally as a fixed source, generated variance reduction parameters, and simulated the prompt doses. Region tallies were used in the model to calculate doses at the desired distances by introducing cells that were cylindrical shells in the actual problem geometry, as was done with MCNP. For the prompt dose calculations total nu-bar was used like the MCNP calculations.
The delayed gamma source was produced with the same coupling of KENO/MONACO/COUPLE/ORIGEN from the SCALE package that the MCNP calculations used. Once the delayed gamma intensities were obtained from the ORIGEN calculations, the fission mesh source was updated to represent delayed gammas instead of fission neutrons by changing particle type and energy spectrum. The spatial distribution of the fission neutron mesh source and delayed gamma mesh source are the same.

COG
COG 11.1 [9] was used with ENDF/B-VII.1 crosssection library data. COG is a general purpose, multiparticle, high-fidelity Monte Carlo code developed by Lawrence Livermore National Laboratory (LLNL). It provides accurate simulation results for complex 3-D shielding, criticality safety, and activation problems. A newly developed feature in COG 11.1 can generate, track and score delayed fission gamma (DFG) rays born between two given times [10]. Point-wise continuous cross-sections are used in COG and a full range of biasing options are available for speeding up solutions for deep penetration problems.
A direct one-step criticality/detector calculation method was applied for all neutron, and prompt and delayed gamma ray dose calculations. Each neutron and gamma particle is tracked from its birth in fission within the spherical fissile volume to its absorption in the system or score at the detectors at various distances in one single, massively parallel, COG supercomputer run with no variance reduction biasing applied. To activate the DFG option, DELAYEDPHOTONS (and associated time interval values) and DGLIB, are input in the BASIC and MIX blocks, respectively. Boundary-crossing detector option was activated to score the dose calculations.

Comparison
Prompt doses comparisons were performed between the 1997 Slide Rule and the modern codes. Figures 4 to 6 show this comparison for three cases and Figures 7 and 8 are focused on the code to code comparison for all cases.    A good agreement is observed between the 1997 Slide Rule results and the modern codes results for both prompt neutron and gamma doses for all cases. The discrepancies between codes are generally below 5%. The corresponding relative error (2 V) on the ratio between two codes is generally lower than 13% for neutrons (whatever the distance) and for gammas for distance lower than 300 m. After 300 m, the corresponding relative error on the ratio for gamma increases and can reach 17% for the MCNP/SCALE comparison and 53% for the COG/SCALE comparison. This discrepancy between codes seems relatively independent of the distance for neutrons but increases for gamma doses with the distance. A possible explanation is that the prompt gamma contribution from the uranium sphere decreases with the distance increase. The gammas being produced by neutrons in the concrete and air, far from the uranium sphere, need to be correctly taken into account. This problem requires very good convergence not only for gammas created inside the uranium sphere, but also for neutrons that create gammas close to the detector. This problem seems more complicated than neutron dose determination and leads to higher discrepancies between codes.
One known difference between MCNP and SCALE gamma transport is MCNP's thick target bremsstrahlung model. This model accounts for the electromagnetic cascade of gammas and electrons that produce many low energy bremsstrahlung gammas. MCNP's thick target bremsstrahlung model accounts for these gammas, and allows users to not perform electron transport for geometries with thick shielding materials. All of SCALE's fixed-source radiation transport codes use gamma production data based on ENDF, which does not include this sort of bremsstrahlung. Regarding COG, whenever an electron-producing photon reaction occurs, COG checks whether the reaction occurred in a region enabled by the user for electron transport. If not (the case here), then the electron energy is immediately deposited.

Additional calculations
Many additional methods can be used to determine the doses, depending on the configurations. In particular, alternative ways to calculate the prompt doses were investigated with MCNP6.1 for some cases. These are divided into two categories.
The first category includes means to create the source and to use it (step one of the "two steps method"). One possibility is to use the Surface Source Write (SSW) card [11] to write the source points of a criticality calculation (KCODE mode), and then use it in a second calculation (via the Surface Source Read (SSR) card). It is an easy way to determine the source but, like the previous method (SMESH), one should pay attention to the convergence of the fission source location distribution, and run enough inactive batches. Besides, one thing to consider is that the continuum of the fission source is being approximated by a discrete set of points. It is then possible to run multiple instances of these same particles with different random number sequences. However it is not possible to generate any new points subsequently. Another limitation is that the sampling of the source variable can't be biased to improve the convergence rate of the problem.
The particle track (PTRAC) card can be used in order to write the starting energy of all energy neutrons and/or gammas produced by fission in order to be used in the second step. As an example, the Figure 9 illustrates the energy distribution of prompt neutron for uranyl nitrate. neutrons/gamma. In this case the resulting "method" looks like, in principle, the SSW/SSR method. This method, more accurate than taking the Watt spectrum, can have an impact on the results for some cases, in particular for highly "heterogeneous" systems (in terms of neutron spectrum and fissile material repartition). Another possibility is to consider, for the second step, the source as a point source, which will save some computer time by not simulating particles that never leave the uranium sphere. This method needs to know the intensity, direction and spectrum of particles leaving the sphere. This method is presented in particular in the reference [12], using MCNP or SCALE. This also can be done with a deterministic code like XSDRNPM, which was used in the 1997 Slide Rule.
The second category includes variance reduction methods used to improve the calculation of the second step.
For example, the use of the F5 Tally (for a point or a ring) may improve the calculation results. It may also help the determination of the first weight windows, with the MCNP weight windows generation (WWG), that can be used later with an F4 Tally. On the contrary, the use of an F4 tally with a small volume (for example small sphere) is sometime necessary but will complicate the calculation. In this case, the use of DXTRAN sphere may be used to improve the particle sampling in a given region around the small volume. This technique and others like forced collision are discussed in [11].
Another possibility to improve the calculation is to use a deterministic code in order to obtain the adjoint flux and to generate the weight windows and the biased source for the Monte Carlo fixed source calculation, like the MAVRIC module of the SCALE package. For MCNP, this step can be done, for example, with ADVANTG [13] or ATTILA [14]. An example of adjoint flux is illustrated in the Figure 10.

Delayed gamma doses
Delayed gamma doses comparisons were performed between the 1997 Slide Rule and the modern codes. Many decay times are considered in references [1] and [6] but only representative cases are presented hereafter. Figures 11 to 16 show these comparisons.    The previous figures present the results for the uranyl fluoride case but the following general observations are applicable for all cases: x the modern code results are higher than the 1997 Slide Rule for times less than 1 minute, but as time increases (to 1 minute) this difference decreases, x at 1 minute, the delayed modern code results agree very well with the 1997 Slide Rule, x at time greater than 1 minute, the modern code results are slightly smaller than the 1997 Slide Rule, but this difference is fairly constant and is not a function of time above 1 minute. These general observations might be slightly different for the different fissile media, because the characteristics of these media (in particular the enrichment and the moderation ratio) have an impact on the evolution of the isotopic composition, so on the delayed gamma source, which leads to some variations (see Figure 16 for example). The complete explanation of these discrepancies is difficult to find because some assumptions made for the 1997 Slide Rule are unknown. However, the following elements can be suggested. The first possibility to explain the discrepancies is the fact that ORIGEN data between 1997 and today has changed, especially with updated nuclear data. In particular for short decay times, the older gamma libraries were missing gamma spectral data for many of the short lived fission products [16]. This might be the best explanation for the discrepancies observed for results for times less than 1 minute. Use of 1997 nuclear data to run ORIGEN today was not tested here but could be done in the future in order to confirm and quantify the nuclear data effect.
It can also be noticed that, the modern code results are obtained by determining the instantaneous dose rate at a given time whereas the 1997 Slide Rule results considered an accumulated one minute dose, which take into account the intensity decrease of the source for short times. A part of the discrepancies between modern codes and 1997 Slide Rule may be explained because the intensity of delayed gamma greatly decreases for short decay time (in particular lower than 1 minute). As an illustration, the intensity of delayed gamma (given in gamma per second) as a function of the decay time after the accident is presented in Figure 17 for the uranyl fluoride case and shows that between 1 second and 1 minute, the intensity decreases by a factor~30. Finally, the 1997 Slide Rule states that "No delayed neutron contribution nor contributions from GHOD\HG JDPPDV EHWZHHQ ȝV DQG V ZHUH LQFOXGHG LQ the dose curves". This assumption should have an impact on the results for the short times.
The code to code comparison is illustrated in Figures 18 to 22.  It can be seen that the discrepancies between codes increase with the distance after 500 meters. It may be explained by large uncertainties in the doses calculation results and not necessarily by a difference in the delayed gamma source. At 10 meters from the criticality accident (with no convergence issues), Figure 22 shows the comparison for all cases. Two behaviors can be observed. The comparison between MCNP and SCALE results shows a slight decrease trend with decay time. A possible explanation is that the delayed gamma sources were determined with the same "method" but with different options (in particular the energy binning) and code versions, which have an impact in particular on the sources intensity for long decay time.
The comparison between COG and SCALE shows an "up and down" behavior. The COG/SCALE parameter presents a slightly important discrepancy (0.82 -0.95) for short time (less than 1 min), then a swing behavior (1.01 -1.15) to finally come back to lower values (0.88 -1.01). This behavior needs to be investigated.

Conclusions and perspectives
AWE (UK), IRSN (France), LLNL (USA) and ORNL (USA) began a long term collaboration effort in 2015 to update the nuclear criticality Slide Rule for the emergency response to a nuclear criticality accident. This document gives order of magnitude estimates of key parameters, useful for emergency response teams and public authorities. The first step of this update is to repeat with modern radiation transport tools, for the same configurations and assumptions, the calculations performed initially for the 1997 estimation of the doses.
Regarding the calculation of the prompt doses, the results from the modern tool used for this update are consistent. For the initial 2D configuration, the modern 3D tools results confirm the results obtained with the 2D tool used to create the 1997 Slide Rule. The impact of the use of new nuclear data is not visible on the update, with regard to the relative error (generally 5%). The interest of 3D capabilities will be developed for new configurations that will break the symmetry of the problem.
Regarding the determination of the delayed gamma source, some discrepancies occurred between codes. The detail analysis of these discrepancies was not the main purpose of this paper but additional effort should be done in order to understand the discrepancies. Other methods/codes should be investigated (like FISPACT or VESTA codes or the ACT card in MCNP) because only two methods were used for this paper.
In addition to the finalization of the calculations and the complete analysis of the results, effort will be devoted on the development and the calculation of new configurations (impact of the geometry and composition of the source, new fissile media including plutonium systems, impact of multiple layers of shielding, impact of the thickness and the composition of the ground, humidity of the air, etc.), using updated flux-to-dose conversion factors (for dosimetry, radiological protection and instrumentation purposes). The goal of this new work is to improve the quality/quantity of information given to the user of the "Slide Rule" in order to not only give a value but also the possible variations and the area of applicability of this value.
In parallel to this Slide Rule update effort (useful to limit the radiological consequences if a criticality accident should occur), this collaborative effort might be a good opportunity to create "computer benchmarks" in order to test and validate the various variance reduction methods and to establish best practices when dealing with this kind of problem (in particular fission source calculation).
Finally, this work will provide the opportunity to suggest experiments allowing the validation of the tool results (benchmarking effort). Indeed, the ICSBEP handbook for this subject is only composed of seven benchmarks, with only one case representative of a criticality accident [15].