URANOS v1.0 – the Ultra Rapid Adaptable Neutron-Only Simulation for Environmental Research

The understanding of neutron transport by Monte Carlo simulations led to major advancements towards precise interpretation of measurements. URANOS (Ultra Rapid Neutron-Only Simulation) is a free software package which has been developed in the last few years in cooperation with particle physics and environmental sciences, specifically for the purposes of cosmic-ray neutron sensing (CRNS). Its versatile user interface and input/output scheme tailored for CRNS applications offers hydrologists straightforward access to model individual scenarios and to directly perform advanced neutron transport calculations. The geometry can be modeled layer-wise, whereas in each layer a voxel geometry is extruded using a two-dimensional map from pixel images representing predefined materials and allowing for the construction of objects on the basis of pixel graphics without a three-dimensional editor. It furthermore features predefined cosmic-ray neutron spectra and detector configurations and also allows for a replication of important site characteristics of study areas – from a small pond to the catchment scale. The simulation thereby gives precise answers to questions like from which location do neutrons originate? How do they propagate to the sensor? What is the neutron's response to certain environmental changes? In recent years, URANOS has been successfully employed by a number of studies, for example, to calculate the cosmic-ray neutron footprint, signals in complex geometries like mobile applications on roads, urban environments and snow patterns.

1 Introduction

The physical processes of neutron transport depend on the atomic composition of materials and the individual neutron energy and act across many orders of spatial scales. It is therefore not feasible to find generalized, analytical solutions under realistic conditions. Statistical and computational approaches are the only way to take all relevant physical interactions into account. In so-called Monte Carlo codes millions of particles can be summoned with randomly sampled initial conditions. Their paths can be tracked, and their interactions with nuclei obey the laws of physics. Finally, the summary statistics of those neutrons can reveal insights into their collective behavior. In the last few decades, the Monte Carlo code MCNP6 (Monte Carlo N-Particle 6) (Goorley et al., 2012) and its predecessor MCNPX (MCNP eXtended) (Waters et al., 2007) were often consulted to study the behavior of neutrons near the surface (Desilets et al., 2006; Zreda et al., 2008; Franz et al., 2013; Zreda et al., 2012; Zweck et al., 2013; Desilets and Zreda, 2013; Andreasen et al., 2016). The conventional model accounts for all kinds of particles and various interactions, decreasing the computational efficiency and resulting in complex model structures (and interfaces). As such it hampers flexible use and is particularly difficult for new users to access. As an alternative for the growing user community of cosmic-ray neutron sensing, we developed the Monte Carlo code URANOS (Ultra Rapid Adaptable Neutron-Only Simulation), which was specifically tailored to address open and recurring questions of the CRNS community (Köhli et al., 2015; Schrön et al., 2017, 2018; Köhli et al., 2018b; Li et al., 2019; Schattan et al., 2019; Weimar et al., 2020; Köhli et al., 2021; Badiee et al., 2021; Jakobi et al., 2021; Francke et al., 2022). As the model has evolved, it has been proven to be useful for neutron spin echo detectors as in other research fields as well (Köhli et al., 2016, 2018a; Köhli and Schmoldt, 2022). URANOS is computationally very efficient as it only accounts for the most relevant neutron interaction processes, namely elastic collisions, inelastic collisions, absorption and evaporation.

The main model features are (1) tracking of particle histories from creation to detection, (2) detector representation as layers or geometric shapes, and (3) voxel-based model extrusion and material setup based on color codes in ASCII matrices or bitmap images.

URANOS is designed as a Monte Carlo tool which exclusively simulates contributions in an environment of neutron interactions. The standard calculation routine features a ray-casting algorithm for single neutron propagation and a voxel engine. The physics model follows the routines declared by the Evaluated Nuclear Data File (ENDF) database standard and descriptions of implementations by OpenMC (Romano and Forget, 2013). It features the treatment of elastic collisions in the thermal and epithermal regime, as well as inelastic collisions, absorption and emission processes such as evaporation, the delayed emission of MeV neutrons from excited nuclei. Cross sections, energy distributions and angular distributions were taken from the databases ENDF/B-VII.1 (Chadwick et al., 2011), ENDF/B-VIII.0 (Brown et al., 2018) and the JENDL/HE-2007 (Japanese Evaluated Nuclear Data Library High Energy 2007) database (Shibata et al., 2011). The entire software is developed in C++ and linked against CERN's analysis toolbox ROOT (Brun and Rademakers, 1997), whereas the GUI uses the QT cross-platform framework.

The graphical user interface offers features specifically tailored to the needs of the field of cosmic-ray neutron sensing. This novel method retrieves subsurface soil moisture by measuring flux of cosmic-ray-induced neutrons that scatter at the soil interface. With typical footprint ranges of hundreds of meters for stationary and beyond 1 km for mobile sensors, it specifically addresses research questions in complex environments.

The paper is divided into 10 sections. After an overview about neutron Monte Carlo codes the physics concepts and the mathematical computation routines in URANOS are discussed in Sect. 2. It is followed by a description of the design ideas of URANOS with the geometrical layer concept and its computational flow in Sect. 4. The individual steps of the calculation are then presented in Sect. 5 with the source configuration and in Sect. 6 with the computation of neutron interactions. Finally Sect. 7 discusses the scoring options. Sections 2 to 7 therefore provide a description of the core of URANOS in terms of computation and design. In Sect. 8 a variety of examples illustrate the capabilities of URANOS and the precision for typical use cases and provide links to previous research questions. Section 9 addresses the computational performance, and Sect. 10 presents the graphical user interface.

1.1 Neutron Monte Carlo codes

The Monte Carlo (Metropolis and Ulam, 1949) method is a brute-force calculation technique, which is used for complex problems consisting of well-defined and/or independent subtasks. It solves a problem by repeated random sampling from a set of initial conditions and interactions.

MCNP was developed as a general-purpose software to treat neutrons, photons, electrons and the coupled transport thereof, excluding magnetic field effects. Versions until MCNP4 (Briesmeister, 2000) were written in FORTRAN 77 (Sun Programmers Group, 1995a), which until the mid-90s was considered the standard in scientific computing. MCNP4 is capable of simulating neutrons up to 20 MeV, which is the maximum of most of the cross sections available in the evaluated databases. With version 5 (X-5 Monte Carlo Team, 2003) the development was forked to the MCNPX (Waters et al., 2007) branch, which converted the code to Fortran 90 (Sun Programmers Group, 1995b) and included the LAHET (Los Alamos High-Energy Transport) (Prael and Lichtenstein, 1989) framework. This revision especially introduced the extension of the energy range for many isotopes up to 150 MeV and some to several GeV by using the Cascade-Exciton Model (CEM) (Gudima et al., 1983) and the Los Alamos Quark-Gluon String Model (LAQGSM) (Gudima et al., 2001). It can also treat (heavy) ion transport for charged particles with energies larger than 1 MeV/nucleon by tabulated interaction ranges. The actual version 6 (Goorley et al., 2012) merged the X-branch into the main development branch. It provides an optional cosmic-ray source (McKinney et al., 2012) which can be used to produce a cosmic neutron spectrum (McKinney, 2013).

A more recent general purpose tool is the PHITS (Particle and Heavy Ion Transport System) code (Iwase et al., 2002), as an extension of the high-energy particle transport code NMTC/JAM (Nucleon Meson Transport code/Jet AA Microscopic) (Niita et al., 2001), which, besides the features mentioned above, also supports charged particles in magnetic fields, d$E/$dx calculations in the CSDA (Continuous Slowing-Down Approximation) (Nelms, 1956) and intra-nuclear cascade JAM transport (Niita, 2002) models up to 1 TeV. PHITS is typically linked against the JENDL databases, consisting of files evaluated by CCONE (Comprehensive COde for Nuclear data Evaluation) (Iwamoto et al., 2016), which is a more sophisticated model compared to INCL (Boudard et al., 2013) and JAM. It comes along with many adjustable parameters for each nucleus, which often leads to a better accuracy compared to other physics models. One of the recent follow-up developments is the PARMA (PHITS-based Analytical Radiation Model in the Atmosphere) (Sato et al., 2008). It calculates the spectra of leptons and hadrons, providing effective models for fluxes of particles of different species, especially with the aim of dose estimations.

The FLUKA (FLUktuierende KAskade) (Battistoni et al., 2015) code is mostly oriented towards charged hadronic transport and nuclear and particle physics experiments. For neutron calculations, the full spectrum is divided into 260 energy groups. FLUKA is not directly linked against an evaluated database, but they operate on their own set of reprocessed and simplified mean values. Especially for neutrons and geometrical representations, it contains reimplementations from the MORSE (Emmett, 1975) neutron and gamma ray transport code.

GEANT4 (GEometry ANd Tracking) (Agostinelli, et al., 2003) can be regarded as FLUKA's successor, based on multithreaded C++ code and OpenGL visualizations. It is designed specifically for the needs of high-energy and accelerator physics. GEANT4 especially excels in describing complex geometries. Since 2011, also driven by requests from the European Spallation Source (Peggs et al., 2013), an increasing number of low-energy neutron calculation features were introduced. Meanwhile the software has advanced to a level where there is a good agreement with other codes like MCNP for fast neutrons (Solovyev et al., 2015), as well as slow neutrons (van der Ende et al., 2016), including cosmic-ray neutron studies (Liu et al., 2021).

1.2 Why another code?

The choice for creating an own independently operating Monte Carlo-based program apart from the mentioned codes was based on evaluating the specific demands of understanding the physics of neutron detectors. The key ideas are as follows:

• Most of the existing codes are not publicly available and fall under the export control law for nuclear-related technology – whereas the underlying databases are free to access. High-precision detector development or environmental sciences are not use cases which are envisaged by the authorities.

• Most of the existing codes were developed in the 1970s or 1980s. Written in the procedural programming language Fortran, which has been proven useful in the ages of limited execution orders and memory, these tools nowadays suffer the drawback of requiring sophisticated and time-consuming code tuning. At best they received wrappers in C, rarely in C++. Today, facing multi-threading, distributed network services and distributed memory in abundance, the changes of computing technology also have a strong impact on the code design and coding strategies.

• Meanwhile complex mathematical operations are readily available from standard packages like the GSL (GNU Scientific Library) (Galassi et al., 2016) and frameworks such as ROOT (Brun and Rademakers, 1997).

• The majority of codes focuses on the evaluation of radiation sources, including gamma emissions. For example, signal generation in a boron-based hybrid detector requires two additional steps of charged-particle transport mechanisms – within the conversion layer itself and subsequently in the gas. In most cases it is not possible to integrate such a calculation path directly, but it would have to be added on top of the simulation. Furthermore, typical codes expect for the geometry objects of roughly equal size – boron layers having an aspect ratio of 10⁵ due to the low thickness are computationally difficult to model.

• All available codes propagate a take-off amount of neutrons in time, as in typical applications like criticality calculations the neutrons themselves change the state of the environment, for example by generating a significant amount of heat. Therefore, the whole ensemble has to be propagated in time, especially until an equilibrium state is reached. Due to limited computing resources this required simplifications like the multigroup method.

• The multigroup method is a technique which allows for the significant improvement of the calculation speed by not treating every neutron track individually but assigning an effective weight to a propagating particle. This weight gets increased for production processes and reduced, if a neutron is absorbed or loses enough energy to drop out of a specific interval. The method is derived from solving Fermi age diffusion equations (Hébert, 2010) and is applied in many codes. However, it requires many interactions to generate enough randomness, and thus it leads to a significant bias in situations when neutrons will likely undergo only 1–2 collisions. For the study of background contributions in detectors or albedo neutrons, such a systematic error should be avoided.

• Restricting the calculation to neutrons and ignoring other nuclear reactions has been proven useful to increase computational speed in programs dedicated to neutron instrumentation and their representation in virtual experiments, like McStas (Lefmann and Nielsen, 1999), VITESS (Wechsler et al., 2000) and RESTRAX (Šaroun and Kulda, 1997).

The only software package which does not suffer from the mentioned drawbacks was GEANT4. Yet when the work on URANOS started in 2014, the GEANT4 code did not at all feature any accurate low-energy neutron calculation. Materials in GEANT4 are usually described under a free gas assumption with unbound cross sections with no information about interatomic chemical bindings. This especially comes into play when treating hydrogen collisions – GEANT4 though can be coupled to the constantly developed models for evaluating the JEFF-3.X (Joint Evaluated Fission and Fusion) (Koning et al., 2011) ACE (A Compact ENDF) formatted thermal scattering law files. For scattering in crystal structures, meanwhile, the NXSG4 (Neutron cross Section library for GEANT4) extension (Kittelmann and Boin, 2015) has been released, which reduced the amount of relevant physics necessary to be integrated. Still, the computational speed of GEANT4 in typical scenarios is significantly lower than those of other codes. In conclusion it has been decided to focus on a design from scratch in a modular, object-oriented language.

2 Calculation routines

2.1 Sampling

The Monte Carlo approach is a stochastic method, in which properties of generated neutrons are randomly chosen from a predefined probability distribution. Examples for sampled neutron transport properties are (a) the path length l, sampled from the probability of an interaction on a distance dx in a homogeneous material of cross section Σ using the random variable r, $l=-\text{ln}\left(r\right)/\mathrm{\Sigma }$, or (b) the thermal neutron velocity distribution.

2.2 Random number generation

The pseudo-random number generator TRandom3 uses the Mersenne Twister algorithm MT 19937 (Matsumoto and Nishimura, 1998) based on the Mersenne prime number 19937. It has a very long period of $p={\mathrm{2}}^{\mathrm{19}\mathrm{937}}-\mathrm{1}\approx \mathrm{4.3}\times {\mathrm{10}}^{\mathrm{6001}}$, low correlation between subsequent numbers (k-distributed for the output sequence) and is relatively fast, as it generates the output sequence of 624 32 bit integers at once – TRandom3 takes approximately 10 ns for each random number on a modern architecture CPU (< 5 years). The generator is seeded at the initialization of the program by the system time in milliseconds. This timestamp is taken as the first integer of the seed sequence, the remaining 623 numbers are generated by the multipliers from Knuth (1997).

2.3 Sampling free path length

The probability p of an interaction along a distance dx in a homogeneous material can be described as

$$\begin{array}{}\text{(1)}& \mathrm{d}p={\mathrm{\Sigma }}_{\mathrm{t}}\mathrm{d}x\end{array}$$

with the macroscopic cross section Σ_t, which in general is energy-dependent. Solutions of this type of differential equation are exponential functions. For the non-interaction probability one therefore can write

$$\begin{array}{}\text{(2)}& p\left(x\right)=\exp \left(-x{\mathrm{\Sigma }}_{\mathrm{t}}\right).\end{array}$$

The probability distribution function for the distance to the next collision (2) assuming conditional probabilities transforms to

$$\begin{array}{}\text{(3)}& p\left(x\right)\mathrm{d}x={\mathrm{\Sigma }}_{\mathrm{t}}\exp \left(-x{\mathrm{\Sigma }}_{\mathrm{t}}\right)\mathrm{d}x.\end{array}$$

The free path length l is obtained by the cumulative probability distribution function of Eq. (3) by

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\underset{\mathrm{0}}{\overset{l}{\int }}p\left(x\right)\mathrm{d}x& =\underset{\mathrm{0}}{\overset{l}{\int }}{\mathrm{\Sigma }}_{\mathrm{t}}\exp \left(-x{\mathrm{\Sigma }}_{\mathrm{t}}\right)\mathrm{d}x\\ & =\mathrm{1}-\exp \left(-{\mathrm{\Sigma }}_{\mathrm{t}}l\right)=P\left(l\right).\end{array}\end{array}$$

In order to retrieve a path length, Eq. (4) can be sampled using the inversion method. This means that the normalized cumulative function is set equal to a random number ξ on a unit interval:

$$\begin{array}{}\text{(5)}& l=-{\displaystyle \frac{\ln (\mathrm{1}-\mathit{\xi })}{{\mathrm{\Sigma }}_{\mathrm{t}}}}\widehat{=}-{\displaystyle \frac{\ln \left(\mathit{\xi }\right)}{{\mathrm{\Sigma }}_{\mathrm{t}}}}.\end{array}$$

As ξ is uniformly distributed in [0,1) the same holds true for 1−ξ, justifying the latter transformation.

It is assumed in Eq. (5) that the material is homogeneous and the cross section along with the kinetic energy stay constant. In case of inhomogeneous material it is possible that the integral cannot be resolved in a closed form. The solution is to split the domain into entities of homogeneous materials and only evaluate the path to the respective border. This procedure is equal to the prerequisite already stated in Eq. (3) that the probability at any point x does not depend on the individual path history.

2.4 Sampling thermal velocity distributions

Scattering thermal processes require sampling a Maxwell–Boltzmann velocity distribution and the relative velocities of neutrons with respect to the target isotope. This also has an influence on the cross section and therefore on the interaction probability. An algorithm has to be applied which preserves the thermally averaged reaction rate, i.e., taking into account the Doppler effect on the cross-section evaluation. Such has been introduced by Gelbard (1979), whereas this modified version follows the implementation by Romano and Forget (2013) and Walsh et al. (2014). By using the effect of thermal motion on the interaction probability

$$\begin{array}{}\text{(6)}& v\stackrel{\mathrm{‾}}{\mathit{\sigma }}(v,T)=\int v_r\mathit{\sigma }\left(v_r\right)f_{\mathrm{M}}^{\left(\mathrm{1}\right)}\left(V\right)\text{d}\mathit{V},\end{array}$$

one has to conserve the reaction rate, the integrand of Eq. (6),

$$\begin{array}{}\text{(7)}& R\left(V\right)=\parallel \mathit{v}-\mathit{V}\parallel \mathit{\sigma }\left(\parallel \mathit{v}-\mathit{V}\parallel \right)f_{\mathrm{M}}^{\left(\mathrm{1}\right)}\left(V\right),\end{array}$$

whereas $f_{\mathrm{M}}^{\left(\mathrm{1}\right)}\left(V\right)$ denotes the velocity distribution for target nuclei of temperature T, velocity V and magnitude of velocity V. The center-of-mass (CM) system of the collision of a neutron with velocity v moves at $v_r=\parallel {\mathit{v}}_{\mathit{r}}\parallel =\parallel \mathit{v}-\mathit{V}\parallel =\sqrt{v^{\mathrm{2}}+V^{\mathrm{2}}-\mathrm{2}vV\cos \mathit{\vartheta }}$. Such a probability function can be constructed by

$$\begin{array}{}\text{(8)}& p\left(V\right)\text{d}V={\displaystyle \frac{R\left(V\right)\text{d}V}{\int R\left(V\right)\text{d}V}}.\end{array}$$

Defining the denominator of Eq. (8) as the normalization factor C and

$$\begin{array}{}\text{(9)}& \mathit{\beta }=\sqrt{{\displaystyle \frac{m}{\mathrm{2}{k}_{\mathrm{B}}T}}},\end{array}$$

as well as μ=cos ϑ, one obtains

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }p(V,\mathit{\mu })\mathrm{d}V\mathrm{d}\mathit{\mu }=\\ \text{(10)}& {\displaystyle }& {\displaystyle \frac{\mathrm{4}\mathit{\sigma }\left(v_r\right)}{\sqrt{\mathit{\pi }}{C}^{\prime }}}\sqrt{v^{\mathrm{2}}+V^{\mathrm{2}}-\mathrm{2}vV\mathit{\mu }}{\mathit{\beta }}^{\mathrm{3}}{V}^{\mathrm{2}}\exp \left(-{\mathit{\beta }}^{\mathrm{2}}{V}^{\mathrm{2}}\right)\mathrm{d}V\mathrm{d}\mathit{\mu }.\end{array}$$

In order to obtain a sampling scheme one can divide Eq. (10) into two parts such that

$$\begin{array}{ll}{\displaystyle }p(V,\mathit{\mu })& {\displaystyle }=g_{\mathrm{1}}(V,\mathit{\mu })g_{\mathrm{2}}\left(V\right)\\ {\displaystyle }{g}_{\mathrm{1}}(V,\mathit{\mu })& {\displaystyle }={\displaystyle \frac{\mathrm{4}\mathit{\sigma }\left(v_r\right)}{\sqrt{\mathit{\pi }}{C}^{\prime }}}{\displaystyle \frac{\sqrt{v^{\mathrm{2}}+V^{\mathrm{2}}-\mathrm{2}vV\mathit{\mu }}}{v+V}}\\ \text{(11)}& {\displaystyle }{g}_{\mathrm{2}}\left(V\right)& {\displaystyle }=(v+V){\mathit{\beta }}^{\mathrm{3}}{V}^{\mathrm{2}}\exp \left(-{\mathit{\beta }}^{\mathrm{2}}{V}^{\mathrm{2}}\right).\end{array}$$

Here the reason for dividing and multiplying Eq. (10) by v+V is that g₁ is bounded. As $\parallel \mathit{v}-\mathit{V}\parallel$ can take on arbitrarily large values, dividing by the sum of the speeds as the maximum value ensures it to be bounded. In general a probability distribution function q(x)=g₁(x)g₂(x) can be sampled by sampling x^′ from a normalized distribution q(x):

$$\begin{array}{}\text{(12)}& q\left(x\right)\mathrm{d}x={\displaystyle \frac{g_{\mathrm{2}}\left(x\right)}{\int g_{\mathrm{2}}\left(x\right)}}\end{array}$$

and accepting it with a probability of

$$\begin{array}{}\text{(13)}& p_{\text{accept}}={\displaystyle \frac{g_{\mathrm{1}}\left(x^{\prime }\right)}{\text{max}\left[g_{\mathrm{1}}\right(x\left)\right]}},\end{array}$$

with g₁(x) bounded. In order to determine q(V) it is necessary to integrate g₂ into Eq. (11):

$$\begin{array}{}\text{(14)}& \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathrm{d}V(v+V){\mathit{\beta }}^{\mathrm{3}}{V}^{\mathrm{2}}\exp \left(-{\mathit{\beta }}^{\mathrm{2}}{V}^{\mathrm{2}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\beta }}}\left(\sqrt{\mathit{\pi }}\mathit{\beta }v+\mathrm{2}\right),\end{array}$$

leading to sampling the probability distribution function

$$\begin{array}{}\text{(15)}& q\left(V\right)\mathrm{d}V=\left({\displaystyle \frac{\mathrm{4}{\mathit{\beta }}^{\mathrm{4}}v{V}^{\mathrm{2}}}{\sqrt{\mathit{\pi }}\mathit{\beta }v+\mathrm{2}}}+{\displaystyle \frac{\mathrm{4}{\mathit{\beta }}^{\mathrm{4}}{V}^{\mathrm{3}}}{\sqrt{\mathit{\pi }}\mathit{\beta }v+\mathrm{2}}}\right)\exp \left(-{\mathit{\beta }}^{\mathrm{2}}{V}^{\mathrm{2}}\right).\end{array}$$

By substituting x=βV, likewise dx=βdV, and y=βv leads finally to

$$\begin{array}{ll}{\displaystyle }q\left(x\right)\mathrm{d}x& {\displaystyle }=[\left({\displaystyle \frac{\sqrt{\mathit{\pi }}y}{\sqrt{\mathit{\pi }}y+\mathrm{2}}}\right){\displaystyle \frac{\mathrm{4}}{\sqrt{\mathit{\pi }}}}{x}^{\mathrm{2}}\exp \left(-x^{\mathrm{2}}\right)\\ \text{(16)}& {\displaystyle }& {\displaystyle }+\left({\displaystyle \frac{\mathrm{2}}{\sqrt{\mathit{\pi }}y+\mathrm{2}}}\right)\mathrm{2}{x}^{\mathrm{3}}\exp \left(-x^{\mathrm{2}}\right)]\mathrm{d}x.\end{array}$$

The terms outside the parentheses are normalized probability distribution functions which allow to be sampled directly, and the expressions inside the parentheses are always <1.

The thermal neutron scattering sampling scheme therefore is the following: a random number ξ₁ is sampled from [0,1) and if

$$\begin{array}{}\text{(17)}& {\mathit{\xi }}_{\mathrm{1}}<{\displaystyle \frac{\mathrm{2}}{\sqrt{\mathit{\pi }}y+\mathrm{2}}},\end{array}$$

the function 2x³exp (−x²) is sampled, otherwise the function $\mathrm{4}/\sqrt{\mathit{\pi }}{x}^{\mathrm{2}}\exp \left(-x^{\mathrm{2}}\right)$. The retrieved x gives the value for V by dividing by β. The decision to accept this velocity is based on Eq. (13). The cosine of the angle can be sampled by another random number ξ₂ in [0,1] by

$$\begin{array}{}\text{(18)}& \mathit{\mu }=\mathrm{2}{\mathit{\xi }}_{\mathrm{2}}-\mathrm{1},\end{array}$$

and as the maximum of g₁ is $\mathrm{4}\mathit{\sigma }\left(v_r\right)/\sqrt{\mathit{\pi }}{C}^{\prime }$ another sampling random number ξ₃ can be used to accept speed and angle by

$$\begin{array}{}\text{(19)}& {\mathit{\xi }}_{\mathrm{3}}<{\displaystyle \frac{\sqrt{v^{\mathrm{2}}+V^{\mathrm{2}}-\mathrm{2}vV\mathit{\mu }}}{v+V}}.\end{array}$$

If this condition is not met, speed and cosine of the angle have to be resampled.

2.5 Nuclear data

Experimental and theoretical results on neutron–nuclear interactions and their subsequent products are collected in libraries. The main database is the Experimental Nuclear Reaction Data Library (EXFOR) (Otuka et al., 2014), which stores most of the accepted published results. Results of neutron interaction measurements are sometimes contradictory and often not comprehensive, therefore so-called evaluated databases have been created. They are the result of literature assessment and intercomparison of different results, yielding standardized values. The databases which are used for this work are the United States Evaluated Nuclear Data File (ENDF/B) and the Japanese Evaluated Nuclear Data Library (JENDL). The high-energy branch JENDL-HE especially provided the largest neutron interaction data set relevant for environmental studies. The ENDF format uses the “MT numbers” to identify neutron reaction types and “MF numbers” to classify the data type of the respective file set (Trkov et al., 2012).

3 Model design

The design of URANOS was motivated by the following general aspects:

• The geometry is represented in a three-dimensional coordinate space with dimensions from the centimeter to the kilometer scale.

• In typical model runs the number of neutrons can easily reach 10⁹ with only one in a million neutron contributing to an observable. For this reason ensemble statistics would not be applicable. Such a calculation method requires a very large number of particles to derive laws without necessarily taking into account each individual state vector.

• The relevant interactions are typically not deterministic but of statistical (random) nature.

• Important parameters like cross sections cannot be derived analytically but have to be extracted from tabulated databases.

• Often neutron interactions can be reduced to a subset of relevant interactions (among them absorption, elastic and inelastic scattering, as well as evaporation), predominately of not more than two different types.

• Most particles other than neutrons are typically not contributing. URANOS, however, is still capable of modeling consecutive conversion ions necessary for signal generation.

• For the creation and propagation of high-energy neutrons within particle cascades effective models can be applied.

3.1 Geometrical layer concept

One specific feature of URANOS is its layer geometry, which takes advantage of the lateral symmetry of typical modeling problems, may it be an air–ground interface or the buildup of a neutron detector. The concept is presented in Fig. 1. While along the horizontal and vertical axes the geometric scales vary significantly, the mean free path lengths in both directions are comparable with respect to the spatial extension. For example the absorption probability for a neutron in a 500 nm film of boron might be around 3 %, the scattering probability in a polymer foil of 100 times the thickness is approximately the same number and in an air volume 100 times larger it may be 0.3 %. This also means that these significant differences in the spatial dimensions are a challenge in terms of defining the geometries for simulation. The solution of URANOS is using layers. This allows us to easily build a geometry of homogeneous materials with the main parameter being the position and height of such a layer. Each layer furthermore can be sub-structured by two-dimensional matrices into voxels. Applying periodic or reflecting boundary conditions to the domain can improve the statistics or reduce the effort for building the simulation model.

Figure 1 Two-dimensional representation of the three-dimensional model setup (not to scale) and processes. Neutrons are generated in a source layer (black-filled circle) following a source spectrum that resembles the energy distribution of particles propagated through the air column above. High-energy neutrons create fast neutrons by nuclear evaporation in the soil (gray-filled circle). They further scatter within (empty red circle) or through the material layers (empty blue circle) until they eventually reach the detector layer (empty gray circle) which absorbs the neutrons depending on the chosen detector response function. While the soil is modeled as a combination of silicate, air and water, various other components of the environment can be modeled with the available material and density options.

3.2 Ray casting

In contrast to other Monte Carlo models focusing on fuel calculations, URANOS uses the method of ray casting in order to keep track of the particles. This improves the accuracy in cases where only a specific subset of conditions will meet the criteria for scoring. Ray casting follows tracks from the source to the point of detection, contrary to ray tracing, which follows tracks backwards from the point of detection but requires mostly deterministic interactions. The ray-casting technique (Roth, 1982) refers to conducting a series of ray-surface intersection tests in order to determine the first object crossed by tracks from a source. These intersections are either defined by analytical surfaces, like the layer structure, or computed from extruded voxels, which do not at all consist of surfaces. Similar types of geometry definitions with mixed volume and surface data were for example used in early computer games when no powerful hardware acceleration was available and nowadays for X-ray tomography image reconstruction in material research (Maire and Withers, 2013), geosciences (Cnudde and Boone, 2013) and especially medical imaging (Goldman, 2007). The method of ray casting also allows us to only record and store the variables necessary for each run. The neutron is propagated forward in time through the domain and flags are used as Boolean operators for each possible output. If for example the recording observable is the flux density above the surface, not the whole track, but only the tracklet within the layer above the ground is kept in the memory.

4 Computational structure

The basic concept of URANOS relies on looping over a set of neutrons, which features initial conditions, predefined or randomized, and consecutively on loop tracking of the path of each neutron through the geometry. Both entities are referred to as “stacks”. In each step the geometrical boundaries are determined and handed over to the physics computation unit. For specific cases actual variables of the neutron or its track history are recorded, emulating a real or a virtual detector. This process is called “scoring”. It can be invoked when passing a specific volume or the track is terminated. A track is defined as the shortest path between two points of interaction. As will be seen later, it can be cut by layer or material boundaries, which dissects it into tracklets. The workflow shown in Fig. 2 illustrates the entire simulation process, which will be described in the following.

Figure 2 Internal workflow of URANOS. Each calculation step is represented by a block describing the structural function in orange and the corresponding physics variables.

4.1 Startup

Before the main calculation routine three steps are carried out:

• assigning memory to objects, which will be used throughout the calculation, by creating empty containers. These are at least 50 one- and two-dimensional root histograms;

• reading the configuration files, creating the geometry and, if available, reading the voxel extrusion matrices; and

• reading the necessary tables from the ENDF library and loading them into the system memory.

The configuration is split into two files, one containing the basic settings for URANOS, like the number of neutrons to calculate and furthermore import and export folders for the data, and one containing information on how to geometrically structure the layers; see here also the next Sect. 4.2. Cross sections and angular distributions are read from tabulated ENDF files, exemplarily shown in Fig. 3, and grouped into absorption, elastic and inelastic scattering. Table 1 exemplarily shows the selected cross sections to be loaded for ¹H, ¹⁰B and ¹⁶O, whereas the full list of available isotopes can be found in Appendix A.

Table 1 Example cross sections according to the ENDF standard.

Only MT numbers with significant contributions are taken into account, which translates to omitting processes with overall less than 10⁻² % of the total cross section. Furthermore, the cross section tables are compressed before it is loaded into the memory. Except for hydrogen, the algorithm skips every consecutive value with a relative difference of less than 1 % to its non-skipped predecessor. This removes 0 % (rare elements) to 98 % (iron) of data, which saves a significant amount of iteration steps of determining the cross section. The smallest error listed on cross sections can be found for elastic scattering of hydrogen with 0.3 %; other isotopes exhibit standard deviations of 1 % and larger, which justifies the compression method.

Figure 3 Examples of cross sections for the light isotopes hydrogen, an efficient moderator; boron as an efficient absorber; and silicon, which can be considered nearly transparent, from the ENDF library from thermal energies ranging to several MeV.

For calculating the total macroscopic cross section the individual contributions of elastic Σ^e and inelastic Σⁱⁿ scattering, as well as absorption Σ^a, are summed up

$$\begin{array}{}\text{(20)}& {\mathrm{\Sigma }}_{\mathrm{t}}={\mathrm{\Sigma }}^{\mathrm{a}}+{\mathrm{\Sigma }}^{\mathrm{e}}+{\mathrm{\Sigma }}^{\mathrm{in}},\end{array}$$

whereas for inelastic cross sections only the main contributors are summed up; see Table A1, and absorption itself is understood as a sum of MT numbers stated in Table 1, which can either lead to its capture without consecutive particles or the creation of new neutrons by for example evaporation or charged-particle ejection by converters.

4.2 Geometry

URANOS uses analytical geometry definitions and voxels as introduced in Sect. 3. The following top-down structure is applied for describing the simulation environment:

$$\mathrm{geometry}\to \mathrm{layer}\to \mathrm{voxelmesh}\to \mathrm{material}\to \mathrm{isotope}.$$

Each layer of the stack is either entirely composed of a material or subdivided into several sections using a two-dimensional matrix from which voxels are extruded. The entities are filled with predefined materials. A material is a specific composition of isotopes with atomic weight and density. Table 2 provides an example of such a definition, whereas all materials available in URANOS can be found in Appendix B. Most compounds are taken from McConn Jr. et al. (2011).

Table 2 Examples of the composition of the material “dry air” (at normal temperature and pressure (NTP)) and a neutron converter.

The voxel mesh is automatically loaded and generated if a file with a name corresponding to a layer number is found. It can be either a tab-separated ASCII matrix of equal row and column rank or a quadratic portable network graphics (PNG) image. The integer values w or grayscale values denote the material numbers which primarily override the global layer definition. Typically solids are directly extruded from these values, yet there are four further declaration modes:

• the material is soil and w defines the amount of water in volume percent;

• the material is soil and w defines the porosity;

• the material is defined globally by the layer or by voxels and w scales the density; and

• the material is defined globally by the layer, w scales the height of this material and the remaining volume extended to the full layer height is filled with air to represent the remaining soil porosity.

The layers can be stacked on top of each other with individual definitions to realize complex geometries. Figure 4 provides examples to illustrate the scope of applications and the scales which can be targeted. The images of one single layer hereby act as a sectional view. Especially landscapes can be modeled using the last declaration mode; an example is provided in Fig. 4.

The geometry of each layer is simply defined by an array of eight elements:

$$\begin{array}{ll}{\displaystyle }g=& {\displaystyle }[x\text{lower bound},x\text{upper bound},\\ {\displaystyle }& {\displaystyle }y\text{lower bound},y\text{upper bound},\\ \text{(21)}& {\displaystyle }& {\displaystyle }\text{upper}z\text{position},\text{height},\text{material},\text{layer number}],\end{array}$$

whereas the lateral lower and upper bounds are defined globally, and the layer number acts as an additional identifier to create subgroups within the stack. Furthermore, the forward and backward propagation direction are defined according to if the layer number along the path increases or decreases, respectively.

Figure 4 Examples of layers for voxel geometry definitions (all in top view) with grayscale values defining preconfigured materials: g cm⁻³ a 2 in. proportional counter with 1 in. of moderator, (b) the rooftop of the institute of physics in Heidelberg, (c) a part of a lake (Schrön, 2017), (d) voxel geometry for a digital environmental model (Kaunertal Glacier at $\mathrm{46}{}^{\circ }{\mathrm{52.2}}^{\prime }$ N, ${\mathrm{10}}^{\circ }{\mathrm{42.6}}^{\prime }$ E) with 500 × 500 pixels at a lateral resolution of 1 m and 0.5 m in height, and (e) shaded illustration of the resulting layered voxel structure from (d).

Neutron tracks S are described by a mixed geometry definition of support vectors x in Cartesian coordinates and spherical direction vectors r:

$$\begin{array}{}\text{(22)}& \mathit{x}=\left(\begin{array}{c}x\\ y\\ z\end{array}\right)\mathrm{and}\mathit{r}=\left(\begin{array}{c}r\\ \mathit{\vartheta }\\ \mathit{\varphi }\end{array}\right),\end{array}$$

denoting the three spatial coordinates $x,y,z$ and the angles ϑ,ϕ with the range r. The choice for this system is due to the fact that this characterization provides direct access to the necessary observables. Examples are point sources which are randomly distributed in both angles or detector planes for which the beam inclination is an important parameter considering sensitivity. Hence new coordinates x^′ are calculated by

$${\mathit{x}}^{\prime }=\left(\begin{array}{c}x\\ y\\ z\end{array}\right)+\left(\begin{array}{c}r\cos \left(\mathit{\varphi }\right)\sin \left(\mathit{\vartheta }\right)\\ r\sin \left(\mathit{\varphi }\right)\sin \left(\mathit{\vartheta }\right)\\ r\cos \left(\mathit{\vartheta }\right)\end{array}\right)$$

and for determining the position on a layer at elevation z_L

$${\mathit{x}}_{\mathrm{L}}=\left(\begin{array}{c}y\\ z\end{array}\right)+\left(z-z_{\mathrm{L}}\right)\left(\begin{array}{c}r\cos \left(\mathit{\varphi }\right)\tan \left(\mathit{\vartheta }\right)\\ r\sin \left(\mathit{\varphi }\right)\tan \left(\mathit{\vartheta }\right)\\ \mathrm{1}\end{array}\right).$$

Time is an indirect quantity. It is derived from the geometrical position of the neutron calculated from energy and initial conditions.

In URANOS three layers can each be assigned specific functions. These are source layer, detector layer and ground layer. The source layer defines the origin for all neutron histories. Especially all height values for starting positions, see also Sect. 5, are restricted to be initiated here. This layer may neither be the upper- nor lowermost as otherwise neutrons would escape the computational domain. The ground layer is used in cosmic neutron simulations to record the spectra at the air–ground interface and monitor the penetration depth. In the detector layer either single real or virtual detectors can be placed, or the layer itself acts as a virtual detector and records every neutron passing; see also Sect. 7.

5 Sources and energy

URANOS provides a variety of sources. A source is defined by a spatial distribution and an energy spectrum from which random values are sampled. As explained in Sect. 4.2 the source has to be placed in the source layer, which defines its height (z position). Available source definitions are the following:

• point sources with all neutrons starting from the same coordinate vector;

• a plane source with all neutrons sharing the same z coordinate within lateral boundaries; and

• a volume source, which randomly distributes neutrons in the source layer within lateral boundaries, and alternatively extends the volume source downwards to the ground layer with exponentially distributed height values.

The cosmic source used for studies of environmental neutrons is typically represented by a plane source. For the coordinates $(x,y)\in A$ in the source area A in case of plane or volume sources the options are the following:

• rectangular boundaries with either equal aspect ratio (square) or any other, sampling the origins uniformly from possible positions in (x,y), and

• circular boundaries, sampling the origins either uniformly in radius r from the center or in (x,y).

Furthermore, the starting angle ϑ can be set to either one of the following:

• full or half sphere, sampling ϑ in $[\mathrm{0},\mathrm{\dots },\mathit{\pi }]$ or $[\mathrm{0},\mathrm{\dots },\mathit{\pi }/\mathrm{2}]$, and

• unidirectional beam, which allows us to set ϑ to a specific inclination. Additionally, a divergence s_ϑ can be chosen. Then, angles are sampled from a Gaussian function centered around ϑ with a width of s_ϑ.

The starting energies are derived from normalized distributions, which are described in the following sections. For source definitions on a linear support in [a,b], like in Sect. 5.2, the random variable $\mathit{\xi }\in [\mathrm{0},\mathrm{1}]$ is scaled to the abscissa test quantity:

$${\mathit{\xi }}_{\mathrm{t}}=a+(b-a)\mathit{\xi }.$$

For source definitions on a logarithmic support in [10^a,10^b], like in Sect. 5.1, ξ is scaled to

$${\mathit{\xi }}_{\mathrm{t}}={\mathrm{10}}^{a+(b-a)\mathit{\xi }}.$$

5.1 The cosmic neutron source

The cosmic neutron source definition is specifically designed for the problem of soil-moisture-dependent neutron transport in the vicinity of the atmosphere–soil interface. Instead of propagating primary particles through several kilometers of atmosphere, a source definition near the ground level is chosen. Recent works, especially from Sato and Niita (2006), Sato et al. (2008) and later Sato (2015), have provided analytical functions modeling cosmic-ray spectra for various conditions like atmospheric depth and cutoff rigidity. In their latest version Sato (2016) introduced a “black hole” mode which enables us to exclusively model the downward-oriented component of the flux, which is used as the default cosmic neutron spectrum. Figure 5 exemplarily shows the URANOS cosmic-ray neutron spectrum (black) and the total spectrum above ground for 10 % soil moisture. The energy of neutrons can range over more than 12 orders of magnitude. The plot here, as well as the following, will be presented logarithmically in units of lethargy. The intensity I or flux density per logarithmic unit of energy is given in units of

$$\begin{array}{}\text{(23)}& I=\mathrm{d}\mathrm{\Phi }/\mathrm{d}\left(\log \right(E\left)\right)=E\mathrm{d}\mathrm{\Phi }/\mathrm{d}E.\end{array}$$

Figure 5 The URANOS cosmic neutron source spectrum (downward only, black) and total angular integrated flux after interaction with the soil 50 m below the source (blue).

5.2 General sources

Besides the cosmic neutron source definition, energy distributions for various specific source types have been implemented. These additional spectra can be used to investigate special applications with regard to sensor development, among them the sensitivity to thermal or epithermal neutrons. These available source configurations allow for sampling from thermal and fission spectra. Whereas the cosmic neutron source would typically be released from an extended layer, usually the geometry for other types of sources may be limited to a considerably small plane, scaling down to even a point source. Exemplarily some are shown in Fig. 6.

• Monoenergetic. Neutrons of energy E or wavelength λ.

• Thermal. Neutrons at a temperature T described by a Maxwellian distribution

  $$\begin{array}{}\text{(24)}& N\left(E\right)={\displaystyle \frac{E}{(k_{\mathrm{B}}T{)}^{\mathrm{2}}}}\exp \left(-{\displaystyle \frac{E}{k_{\mathrm{B}}T}}\right).\end{array}$$

• Predefined. Americium–beryllium spectrum from ISO group 85/SC 2 Radiological protection (2001).

• Evaporation. Assuming the nucleus to form a degenerate Fermi gas (Weisskopf, 1937) one can derive various forms of density distributions in the form

  $$\begin{array}{}\text{(25)}& N\left(E\right)\propto E\exp \left(-{\displaystyle \frac{E}{k_{\mathrm{B}}T}}\right),\end{array}$$

  which are simply described by a temperature parameter (Terrell, 1959). The energy distribution of the neutrons released by fission are commonly represented either by a Maxwellian distribution or the following Watt spectrum (Iyer and Ganguly, 1972).

• Fission. A semi-empirical description for fission neutrons is the Watt spectrum (Watt, 1952), especially used for ²³⁵U, which can be selected as a source although the isotope itself is not implemented. The following spectra can be selected:

  $$\begin{array}{}\text{(26)}& N\left(E\right)=\mathrm{0.4865}\sinh \left(\sqrt{\mathrm{2}E}\right)\exp (-E)\end{array}$$

  and for ²⁵²Cf (Smith et al., 1957)

  $$\begin{array}{}\text{(27)}& N\left(E\right)=\sinh \left(\sqrt{\mathrm{2}E}\right)\exp (-\mathrm{0.88}E),\end{array}$$

  which are both specific cases of the more general form of a Maxwellian distribution. A more accurate modeling can be performed by specifying the Watt parameters a and b, taking into account the mean neutron kinetic energy of and those of the fission fragments:

  $$\begin{array}{}\text{(28)}& N\left(E\right)={\displaystyle \frac{\exp \left(-ab/\mathrm{4}\right)}{\mathrm{0.5}\sqrt{\mathit{\pi }{a}^{\mathrm{3}}b}}}\sinh \left(\sqrt{bE}\right)\exp \left(-{\displaystyle \frac{E}{a}}\right).\end{array}$$

  The parameters a,b are usually tabulated as a function of energy, element and isotope in ENDF libraries.

Figure 6 Different preconfigured source distribution functions in URANOS in the MeV range covering different use cases: laboratory test sources like americium–beryllium (yellow), spontaneous fission (green), evaporation (cyan) from nuclear de-excitation and fusion (blue).

6 Calculation scheme

6.1 Loop nodes

The main calculation routine runs in two loops; for each neutron in the neutron stack, the geometry stack is traversed layer-by-layer:

$$\textcolor{red}{}.$$

Each onset neutron is a placeholder and only initialized at runtime with the particle number not being conserved due to physical processes generating further neutrons. The layer stack is created at startup and consists of a fixed amount of elements which are traversed by an iterator either forwards or backwards, depending on the spatial direction vector.

The possible initial conditions for neutrons are the following:

• Energy. Available source definitions from Sect. 5, which can be either real values, normalized functions to be sampled from or lookup tables.

• Geometry. Definition from Sect. 4.2, which can be either a fixed vector from a source or a distribution function to be sampled from or lookup tables, which are normalized at startup.

Using these initial conditions the loop over the layer stack commences. Each layer, which is geometrically described in Sect. 4.2, can either consist of a homogeneous material defined by its isotope composition, a material defined by an analytical function, or an input matrix from which voxels are extruded. A comprehensive material list is provided in Appendix B. The neutron iterates to the following layer if it geometrically leaves the boundaries without absorption, and no change of materials can be found in the collision detection.

6.2 Tracking in finite geometry regions

For each layer the material composition is loaded according to the actual position of the neutron. The definition either accounts for the whole layer or for regions, which can be described by analytic functions or voxels. For the selected material the total macroscopic cross section Σ_t is set up isotope-by-isotope. The amount and type of reactions (MT identifiers) depends on the element; see also the description in Sect. 4.1 and the isotope list in Appendix A. Elemental hydrogen for example cannot undergo inelastic scattering and ¹⁰B exhibits a negligible radiative capture probability, so only charged reaction paths are relevant. The selection criteria in detail are as follows:

• elastic and absorption cross sections are always calculated if available, and

• inelastic cross sections are loaded on demand for energies 750$\mathrm{keV}<E<$50 MeV.

The macroscopic cross section of a compound with weight fractions w_i of n elements and energy dependent absorption σ^a and scattering σ^s cross sections is defined as

$$\begin{array}{}\text{(29)}& {\mathrm{\Sigma }}_{\mathrm{t}}=\mathit{\rho }{N}_A\sum _{i=\mathrm{1}}^n{w}_i{\displaystyle \frac{{\mathit{\sigma }}_i}{M_i}}\text{with}{\mathit{\sigma }}_i\left(E\right)={\mathit{\sigma }}^{\mathrm{a}}\left(E\right)+{\mathit{\sigma }}^{\mathrm{s}}\left(E\right).\end{array}$$

Other cross sections can also contribute. The free path length l is sampled from a random number ξ as described in Sect. 2.3 from Eq. (5): $l=-\ln \left(\mathit{\xi }\right)/{\mathrm{\Sigma }}_{\mathrm{t}}$. An example of the total cross-section evaluation is shown in Fig. 7 for the two typical materials air and water. In case the material definition contains a density multiplication factor, it is applied to Σ_t before evaluating Eq. (5). The distance to the border l_trj is calculated by the z coordinate of the last interaction of the neutron z₀ and the layer z position z_l and height d_l.

Figure 7 Mean free path $\mathrm{1}/{\mathrm{\Sigma }}_{\mathrm{t}}$ for neutrons in the MeV range. The dominant peaks originate from the contribution of the elastic scattering cross section in dry air (NTP) mainly by nitrogen and in water by oxygen; see also Fig. 8.

In case the material is defined by voxels, additionally a procedure is applied which samples the trajectory according to the underlying pixel matrix.

• Determination of the z-projected length z_m of one lateral unit pixel s_p for the actual direction vector by $z_{\mathrm{m}}=s_{\mathrm{p}}/\tan \left(\mathit{\vartheta }\right)$. The unit pixel size is determined by the spatial extension of the domain divided by the number of pixels.

• If the material of the voxel at x^′ for z₀±z_m is different from the actual, stop and repeat the range calculation for the actual composition and geometry. If the material does not change, iterate ±z_m until the layer border is reached. The propagation direction, forward or backward, determines sgn(z_m).

If l_trj>l no interaction takes place, and the neutron can proceed to the following layer. It has to be noted that for a voxel-based geometry definition the resolution is the size of a pixel. A neutron in that case will not interact with the side faces of a voxel but with its volume. If l_trj<l the spatial coordinates of the interaction x_i are calculated by

$$\begin{array}{}\text{(30)}& {\mathit{x}}_i=\left(\begin{array}{c}{x}_{\mathrm{0}}\\ y_{\mathrm{0}}\\ z_i\end{array}\right)+\left(\begin{array}{c}\cos \left(\mathit{\varphi }\right)\left|\tan \right(\mathit{\vartheta }\left)\right(z_i-z\mathrm{0}\left)\right|\\ \sin \left(\mathit{\varphi }\right)\left|\tan \right(\mathit{\vartheta }\left)\right(z_i-z\mathrm{0}\left)\right|\\ \mathrm{0}\end{array}\right),\end{array}$$

with afterwards updating the new z coordinate.

Consecutively, the type of reaction is determined by another random number ξ. The relative fractions of the respective macroscopic cross section define the probability:

$$\begin{array}{ll}{\displaystyle }\mathit{\xi }<{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{el}}}{\mathit{\sigma }}}& {\displaystyle }⟶\mathrm{\text{scattered elastically}},\\ {\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{el}}}{\mathit{\sigma }}}<\mathit{\xi }<{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{el}}+{\mathit{\sigma }}_a}{\mathit{\sigma }}}& {\displaystyle }⟶\mathrm{\text{absorbed}},\\ \text{(31)}& {\displaystyle }\mathit{\xi }>{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{el}}+{\mathit{\sigma }}_a}{\mathit{\sigma }}}& {\displaystyle }⟶\mathrm{\text{scattered inelastically}},\end{array}$$

as the probability of selecting a reaction type i from all possible interaction channels is

$$\begin{array}{}\text{(32)}& p_i={\displaystyle \frac{{\mathrm{\Sigma }}_i}{\sum _{j=\mathrm{0}}^n{\mathrm{\Sigma }}_j}}={\displaystyle \frac{{\mathrm{\Sigma }}_i}{{\mathrm{\Sigma }}_{\mathrm{t}}}}.\end{array}$$

The target interaction element is determined by randomly choosing from a proportional lookup table. Each reaction type is accompanied by two vectors – one number represents the cumulative probability distribution v_s and one number a list of corresponding elements v_e:¹

$$\begin{array}{}\text{(33)}& {\displaystyle }{v}_{\mathrm{s}}\left[n\right]& {\displaystyle }=\left\{i|\sum _{j=\mathrm{0}}^i{\mathrm{\Sigma }}_j\right\},\text{(34)}& {\displaystyle }{v}_{\mathrm{e}}\left[n\right]& {\displaystyle }=\left\{i|\text{isotope of }i\right\}.\end{array}$$

For inelastic scattering additionally the excited state of the target isotope has to be determined. One vector of numbers contains the cumulative cross-section distribution and two support vectors contain the q values representing the energy loss in MeV and the inelastic angular distributions in a <TMatrixF> list. The individual contributor, which will be used to determine the reaction target, is chosen by a random number ξ. If

$$\begin{array}{}\text{(35)}& v_{\mathrm{s}}\left[i\right]\le \mathit{\xi }\le v_{\mathrm{s}}[i+\mathrm{1}],\forall i>\mathrm{0},\end{array}$$

then the corresponding isotope is taken from v_e[i].² Figure 8 schematically illustrates the calculation routine for all processes described above.

Figure 8 Range calculation in URANOS is as follows: for a given neutron energy, here in the MeV range, the cross sections from the isotope list are evaluated according to elastic, inelastic and absorption processes. Only possibly relevant contributions are evaluated. The left panel shows such a list of reaction probabilities for water. Inelastic levels are only displayed up to MT56 and are linked additionally to energy loss and angular distribution. The cross section multiplied by the atom number density in Eq. (29) yields the macroscopic cross section. By sampling a random number ξ as in Eq. (5), a free path length value for the range (4) is obtained.

6.3 Interaction channels

For each interaction the following quantities are updated:

• the position vector x, including time, by adding the path length l to the last position;

• the direction vector r; and

• energy, including velocity v and wavelength λ.

6.3.1 Elastic and inelastic scattering

Scattering is described by the collision of a neutron with a nucleus of mass A assuming energy and momentum conservation. The problem has a radial symmetry regarding the impact parameter, therefore only one angle ϑ_CMS needs to be calculated. The second angle can be determined by a random number ξ in [0,1):

$$\begin{array}{}\text{(36)}& {\mathit{\varphi }}_{\mathrm{cm}}=\mathit{\pi }\left(\mathrm{2}\mathit{\xi }-\mathrm{1}\right).\end{array}$$

For inelastic scattering the energy loss is substituted by the q value obtained from Eqs. (35) and (34). The target velocity V can be neglected for kinetic energies E of the neutron:

$$V\approx \mathrm{0}\mathrm{if}\left\{\begin{array}{ll}\mathrm{0.11}\mathrm{eV}<E<\mathrm{1}\mathrm{MeV}& \text{in case of hydrogen,}\\ \mathrm{0.15}\mathrm{eV}<E<\mathrm{0.01}\mathrm{MeV}& \text{otherwise. }\end{array}\right.$$

For lower energies the interaction result has to be calculated by laws of thermal scattering, taking into account the velocity distribution of the target material. In the case of amorphous material or fluids there is no analytical form to describe such a state; therefore, only sampling from an effective thermal spectrum like Eq. (24) is carried out. For solids with a crystal lattice Bragg scattering is the dominant channel. The kinetic theory of gases allows for a cohesive description of the scattering process. For such calculation the energy and angle are sampled according to Eqs. (17), (18) and (19).

In case of higher energies than stated in the above limits the angular distribution in the center of the mass frame can be found in ENDF cards either tabulated or described by Legendre polynomials. With increasing energy the forward direction is preferred, except for hydrogen – here the asymmetry is much weaker than for heavy elements and only for very high energies a significant deviation from an even distribution can be observed.

For inelastic scattering with an energy transfer E^* the evaluation of the angular distributions is carried out likewise, whereas the lowest possible energy is given by the first q value. Hence, the reaction kinematics of inelastic processes share some similarities with elastic processes of corresponding kinetic energies $E^{\prime }=E-E^{*}$.

As the scattering kinematics have been calculated in the center of mass system, a transformation to the laboratory system is carried out via

$$\begin{array}{}\text{(37)}& {\mathit{\vartheta }}_{\mathrm{l}}=\arccos \left({\displaystyle \frac{\mathrm{1}+A\cos \left({\mathit{\vartheta }}_{\mathrm{cm}}\right)}{\sqrt{A^{\mathrm{2}}+\mathrm{1}+\mathrm{2}A\cos \left({\mathit{\vartheta }}_{\mathrm{cm}}\right)}}}\right)\end{array}$$

and added to the existing direction vectors

$$\begin{array}{ll}{\displaystyle }{\mathit{\vartheta }}_u& {\displaystyle }=\cos \left({\mathit{\vartheta }}^{\mathrm{o}}\right)\cos \left({\mathit{\vartheta }}_{\mathrm{cm}}\right)\\ \text{(38)}& {\displaystyle }& {\displaystyle }+\sin \left({\mathit{\vartheta }}^{\mathrm{o}}\right)\sin \left({\mathit{\vartheta }}_{\mathrm{cm}}\right)\cos \left(\mathit{\pi }+{\mathit{\varphi }}_{\mathrm{cm}}\right),\text{(39)}& {\displaystyle }{\mathit{\vartheta }}_{\mathrm{l}}^{\mathrm{new}}& {\displaystyle }=\arccos \left({\mathit{\vartheta }}_u\right),\text{(40)}& {\displaystyle }{\mathit{\varphi }}^{\mathrm{new}}& {\displaystyle }={\mathit{\varphi }}^{\mathrm{o}}\pm \arccos \left({\displaystyle \frac{\cos \left({\mathit{\vartheta }}_{\mathrm{cm}}\right)-\cos \left({\mathit{\vartheta }}^{\mathrm{o}}\right){\mathit{\vartheta }}_u}{\sin \left({\mathit{\vartheta }}^{\mathrm{o}}\right)\sin \left({\mathit{\vartheta }}_{\mathrm{cm}}\right)}}\right).\end{array}$$

Due to the choice of the coordinate system, see also the geometry definition (22), adding direction vectors is less convenient than the otherwise usual declaration. The method presented here equals to an Euler rotation in θ and ϕ around the direction axis given by the trajectory of the particle.

6.3.2 Evaporation

URANOS simplifies the calculation of the evaporation process, as in the low-Z and intermediate energy range most quantities relevant for fissionable elements are approximately constant. The mean number of evaporated neutrons can be considered constant ${\stackrel{\mathrm{‾}}{n}}_{\mathrm{evap}}\approx \mathrm{1}$ for projectile energies below several hundred MeV and mass numbers of A<100 (Cugnon et al., 1997). Furthermore, for the emission energy a Maxwellian spectrum according to Eq. (25) with a mean neutron energy of 1.8 MeV (Kawano et al., 2013) and a flat angular distribution (Bramblett and Bonner, 1960) is assumed. In order to provide upper limits in comparison ²³⁵U produces on average $\approx \mathrm{2.4}+E$ MeV⁻¹ neutrons per fission.

6.3.3 Absorption

Neutrons are either absorbed by a non-radiating process and consequently the calculation is terminated, or the material is a specific absorber, which leads to a scoring by the detection unit. A specific case is the high-energy cascade: URANOS mainly carries out neutron interactions. For the generation of high energetic radiation in the atmosphere charged particles are also largely contributing to the production of the neutron component (Sato, 2016; Ziegler, 1998). As far as for low energetic and albedo neutrons such can be neglected; in order to simulate more than 100 m of atmosphere the generation of the primary spectrum is emulated by an effective model: for any interaction occurring above 16 MeV with the possible release of secondary neutrons the primary neutron is not eliminated if a random number ξ is below a specific value k_HE, receiving only a fractional energy loss and angular deviation. This value k_HE is tuned to emulate an effective atmospheric attenuation length L_prim of the primary spectrum component of 145 cm² g⁻¹. Experimental values for L_prim are in the range of (135–155) cm² g⁻¹, depending on the latitude of the site. Here, the value from Sato (2016) is taken.

7 Detector configurations

Neutrons can be scored in three different ways: one layer can be defined as the “detector layer”, a virtual entity which can record any particle of chosen characteristics to pass through. Secondly, a virtual detector with limited spatial extension can be placed within the detector layer. Specific materials, if voxel definitions are used, can additionally score neutrons like the virtual detector. For a thermal neutron detector this material would be the converter. Output options are either only hits, partial tracks within a material or full tracks. The detector layer stretches out to the full domain dimension and is placed at a fixed height. It is most useful for mapping the spatial distribution of the neutron flux. The virtual detector can emulate an instrument at a specific position and can be used for analyzing the origin distribution of neutrons. It can be set to transparent for a generalized recording or to absorbing for mimicking a detector more realistically.

Figure 9 Detection efficiencies for the standard CRNS detector model. (a) Energy-dependent absorption probability for perpendicular irradiation, here simulation of a monoenergetic beam with results (red markers) averaged over the surface. (b) Energy independent, averaged, angular dependence relative to (a).

Scoring options for CRNS

In the most simple case a uniform detection efficiency ϵ can be chosen within a specific range:

$$\begin{array}{}\text{(41)}& \mathit{ϵ}=\left\{\begin{array}{ll}\mathrm{1}& \text{for }{E}_{\min }<E<E_{\max },\\ \mathrm{0}& \mathrm{otherwise}.\end{array}\right.\end{array}$$

If thermal neutrons are not considered, the flux in the epithermal/fast region can be considered a plateau region, partially justifying the established choice of Eq. (41). In order to not model an entire cosmic-ray neutron detection system in an environment larger by orders of magnitude, a set of “standard” detectors has been simulated independently and integrated as an effective model. In many possible scenarios the “simple” scoring option (41) leads to incorrect results. This approach therefore allows us to directly analyze the neutron flux measured by an actual device and at the same time allowing us to scale up the scoring volume beyond the physical dimensions of the instrument. Figure 9 shows the implemented functions, which represent averaged values for different side faces of a detector.

In URANOS cubic spline interpolation is used for describing the absolute efficiency and the angular dependence is modeled by

$$\begin{array}{}\text{(42)}& {\mathit{ϵ}}_{\mathit{\vartheta }}=\mathrm{1.24}-\mathrm{0.254}\exp \left({\displaystyle \frac{\mathit{\theta }}{\mathrm{0.92}}}\right).\end{array}$$

The options above can be applied to the whole detector layer. If for example angular resolution or single neutron tracking information are required, one can place two types of scoring units within the detector layer, which is either a plane in z direction, a sphere or a vertical cylinder, whereas in both cases the radius can be specified. The cylinder height corresponds to the detector layer. If due to positioning or the choice of the radius of the sphere an intersection with the layer boundary occurs, only the volume inside the detector layer is taken into account.

8 Evaluation of the URANOS model

As the implemented physics is similar, model results are expected to be comparable to other codes. In preparatory studies, we explored the performance of the URANOS model in reproducing results from standard software like MCNP(X). The tests successfully agreed in many different setups as the one presented by Köhli et al. (2021). However, more rigorous validation experiments may be conducted in future projects using MCNP and GEANT4. Apart from that, URANOS model results have found remarkable agreement with observations, e.g., in terms of footprint experiments (Köhli et al., 2015) and applications to many studies in cosmic-ray neutron sensing (Heidbüchel et al., 2016), while it should be noted that the major challenges with these comparisons are the experimental uncertainties of soil profiles and heterogeneity, among others (Baroni et al., 2018; Iwema et al., 2021). Furthermore, the physics of neutron attenuation in water can be validated with measurements from Caswell et al. (1957), conducting measurements in water tanks.

8.1 Basic performance examples

In order to visualize the tracking capabilities of URANOS Fig. 10 shows two non-trivial neutron paths from generation until absorption, exemplarily in air (top) and in the ground (bottom). It acts as a demonstrator for the interactions at this specific interface. In air the main scattering partners are nitrogen and oxygen, which leads to a large amount of scatterings with small energy decrements. By the long path lengths in the less-dense medium the neutron can also acquire hundreds of meters of integrated travel distance. Inside the soil typical scattering lengths are far below 1 m. For high-energy neutrons, the main scattering partners can be silicon, aluminum and oxygen. However, due to the presence of water a few interactions with light nuclei can thermalize a neutron (blue lines). Then it will carry out a random walk which will be dominated by hydrogen scattering.

Figure 10 Projection of track calculations in an air–ground interface. The simulated neutrons, which are artificially released from 1 m above the soil, are rainbow-colored according to the logarithm of the corresponding energy scaling from 10 MeV (red) to thermal (blue). (a) A neutron which mainly scatters in the air. (b) A neutron thermalizing inside the soil. To be noted: both x and y axes are not scaled equally.

8.2 Diffusion length in water

The attenuation of fast neutrons by efficient moderators is a basic example of neutron physics and the main source of thermal neutrons. Modeling the slowing-down process properly requires the correct description of interaction lengths, energy loss and geometric transport. Therefore, it can be regarded as a validation test of the Monte Carlo code. Only a few examples of well controlled and simple measurements can be found in available literature. Caswell et al. (1957) describe an experiment of determining the radial distribution of neutrons in a water tank from 14.1 MeV to thermal energies and 1.46 eV. A deuterium beam is delivered by an aluminum tube onto a tritium target, inducing fusion. The tank measures 2.4 m in length and 1.2 m in height, whereas the particle injector is located at a distance of 0.6 m from one wall and vertically centered. The flux is measured pointwise by indium foil activation, which provides data for the non-equilibrium state above 1 eV, and thermal neutron detectors with cadmium shielding. Although both energy regimes are supposed to exhibit similar range distributions, they have to be treated by different methods of neutron transport. Until reaching the indium resonance a maximum mean energy loss by elastic collisions, including a few inelastic reactions, can be attributed to hydrogen interactions. Below this limit the kinetics of neutrons are dominated by thermal scattering, leading to a constant average energy. Figure 11 shows the measured fluxes from Caswell et al. (1957) in comparison to the simulation results. Both attenuation distributions are in good agreement. The particle density in both cases peaks at around 15 cm followed by a nearly exponential decay with similar attenuation lengths.

Figure 11 Comparison of the attenuation length from Caswell et al. (1957) for deuterium–tritium fusion neutrons emitted into water. The spherical surface flux for thermal and indium resonance neutrons as a function of distance from the source is compared to the simulation results from URANOS.

8.3 Simulation of the response function of Bonner spheres

A similar case are Bonner spheres (Bramblett et al., 1960), proportional counters surrounded by shells of polyethylene. As this spectrometer type of array is used to monitor environmental fluxes, various studies were carried out for the modeling of such (Hertel and Davidson, 1985; Mares and Schraube, 1994; Garny et al., 2009; Decker et al., 2015). Whereas the neutron range distribution in water in the previous example demonstrated geometric transport and collision treatment, the Bonner sphere offers the possibility to focus on an energy-dependent comparison and on the interplay of moderator and absorber. Among the various existing technical realizations the helium-based version was chosen, equipped with a 3.2 cm spherical counter. For reasons of convenience, the whole model has been discretized in 17 layers, which are symmetrically arranged around the center. Laterally the resolution by the pixel matrix was set to 1 mm, therefore the voxel size of an X in. sphere is 1 mm × 1 mm × (X/17) in. For the simulation the model was irradiated by a neutron beam of the same diameter as the sphere under an angle of 0^∘.

Figure 12 Flux calculation of Bonner spheres of 2 in. diameter. The simulated neutron tracks (E_kin=10 keV) of 10⁶ histories are displayed in a central cross section of 3 mm height (top row) and the full domain of 13 cm × 13 cm × 5.4 cm (bottom row).

Furthermore, hydrogen atoms in polyethylene have been emulated by the scattering kernel derived from the (oxygen-)bound cross section in water. This can yield exclusively for thermal energies a systematic uncertainty of around 10 %. Due to the statistical nature of neutron transport the actual geometrical shape of a body has a minor influence compared to other parameters like overall volume or thickness. Exemplarily for the calculation routine some track views are shown in Fig. 12 as a central cut through the model and for the whole domain.

Figure 13 Comparison of simulations of the energy-dependent response function of Bonner spheres of URANOS and MCNP calculations by Mares et al. (1991). The detectors are HDPE (high-density polyethylene) spheres with diameters in the range of (2–5) in, equipped with a 3.2 cm ³He counter.

Comparing to calculations from Mares et al. (1991), see Fig. 13, there is a good agreement in the energy sensitivity between response curves from literature and URANOS results. This successfully validates the simulation for basic scattering calculations.

8.4 Cosmic spectrum evaluation

Although for more than 50 years the general shape and height-dependent scaling of the cosmic-ray neutron spectrum at ground level has been known (Yamashita et al., 1966), there is a perpetual discussion about the precise features of the intensity distribution, especially at the soil interface. The reasons are that high-energy neutron interaction cross sections above 20 MeV were originally not seriously investigated nor integrated into transport codes. Their evaluation and corresponding measurements are recent developments, mainly of the 21st century. Furthermore, the invention of the Bonner sphere could standardize dosimetric flux evaluations, yet, the neutron spectrum is measured indirectly. The flux distributions are the result of unfolding algorithms (Sweezy et al., 2002), which rely only on a few absolute values and energy-dependent response functions from Monte Carlo models of the detectors (Thomas and Alevra, 2002). This means that different simulations can produce slightly different weightings for different parts of the spectrum. In Fig. 14 an overview of different results from the most widely used codes is presented along with experimental results. The main differences appear in the high-energy regime. These uncertainties were partly compensated by effective nuclear interaction models. The peak structures at around 1 MeV, which are spectral lines of (in)elastic resonances, mostly oxygen, cannot be resolved experimentally by spectrometers and are displayed only at times.

In order to minimize this general problem URANOS uses a validated neutron spectrum near the surface as a source. It is released directly onto the ground to reduce typical uncertainties of atmospheric propagation. The implementation of the works presented by Sato (2016), which are based on Iwase et al. (2002) and Sato et al. (2008), is discussed in Sect. 5.1. Fig. 15 presents the result from URANOS for the calculated neutron flux (green) above the surface in an infinite domain. Note that, although atmospheric propagation of cosmic rays are not necessary within the URANOS concept, there is still the option to model any variation of atmospheric conditions using the layer, material and density definition system.

Figure 14 Energy-dependent neutron flux at a 1 m altitude calculated by PARMA, MCNPX, and GEANT4 and determined experimentally from Goldhagen et al. (2002) and from Federico et al. (2010). Environmental conditions are not the same (Pazianotto et al., 2015).

Figure 15 Cosmic-ray neutron spectrum comparison for 3 % (light colors) and 50 % (dark colors) soil moisture. The spectrum, which is generated from the analytical functions from Sato (2016), is shown in blue. The black hole mode spectra from Sato (2016), see also Fig. 5, are used as input distributions at a height of 500 m and released onto the soil. The resulting intensity distribution is shown for MCNP (black) and URANOS (green) in comparison.

On the qualitative level the underlying physics model correctly calculates the response to the soil. Three peaks can be observed – the high-energy domain around 100 MeV with the incoming-only flux, the region around 1 MeV with the evaporation peak and the thermal peak at 25 meV. As there are no significant sources in the range of 1 eV to 0.1 MeV, there is a flat plateau between neutron generation and thermalization. This plateau can feature a slant angle in case there are significant absorption processes involved. An incision between high energy and evaporation peak can occur from not using cascades but an effective model for neutrons above 20 MeV. The high-energy cascade contributes to the production of neutrons by three different particle species with three different fractions and attenuation lengths (around 140 cm² g⁻¹ for neutrons, around 110 cm² g⁻¹ for protons and around 500 cm² g⁻¹ for muons). Future work will address these different attenuation lengths and model the high-energy cascade accordingly. Due to the lack of a generally accepted standardized spectrum or a consensus in the literature, the evaluation of the URANOS code focuses on the capability to reproduce the above-ground cosmic neutron spectrum for typical conditions. This implies that the input spectrum released on the ground should reproduce the same densities as the input formulae (blue) or other simulation toolkits like MCNP (black).

8.5 Previous studies using URANOS

The URANOS model has been in active development since 2014 while it has been employed by many international studies to support the CRNS data analysis and interpretation. We here provide a short summary of the published studies which indirectly confirm the proper modeling capabilities of URANOS by their comparisons to real-world field observations.

In the first application, Köhli et al. (2015) simulated the radial CRNS footprint and confirmed the results with measurements at the shoreline between land and a lake. Köhli et al. (2016) used a URANOS model of the boron-10-based CASCADE detector including signal generation in the gas and its projection onto the readout structure in order to understand the signals generated specifically during spin echo measurements. The results agreed with reference efficiency measurements and helped to improve the event selection algorithm. The model was also employed by Schrön (2017) to find the optimal location of a buoy detector on a lake with minimal influence from the land. By challenging six different CRNS sites, Schrön et al. (2017) consistently improved their soil moisture measurement and calibration performance with the URANOS predictions of distance and depth contribution. This weighting approach has later been confirmed by many other authors worldwide. Schrön et al. (2018) simulated the spatial neutron density in an urban environment and confirmed agreement to corresponding observations. In Köhli et al. (2018a) URANOS was employed to generate reference track and energy deposition models in order to fine-tune the event reconstruction in the BODELAIRE (BOron-based DEtector with Light And Ion REconstruction) detector with its Timepix readout. Schrön et al. (2018) used URANOS to describe the local road effect for mobile CRNS applications and were able to experimentally verify and correct these effects using the predicted relationships. Köhli et al. (2018b) simulated detector response functions of various neutron detectors used in the context of cosmic-ray neutron sensing and found very good agreement with previously published data and models. Li et al. (2019) were able to explain the observations in an irrigated orchard with the help of detailed URANOS simulations. Schattan et al. (2019) discovered the hysteresis effect in CRNS signals of complex snow patterns and explained this effect accurately with dedicated URANOS simulations. Weimar et al. (2020) employed URANOS to develop new CRNS detector and shielding models which eventually improved the sensor performance. Köhli et al. (2021) revised the neutron–water relationship N(θ,h) and found substantially improved performance for CRNS observations particularly in dry regions. Badiee et al. (2021) employed URANOS simulations to support new developments of downward-shielded CRNS robots. Jakobi et al. (2021) investigated and evaluated the footprint characteristics of thermal neutron transport using URANOS. Francke et al. (2022) used URANOS to explain the angular distribution of measurements with a directional-sensing CRNS variant. Rasche et al. (2021) were able to explain unconventional behavior of neutron transport at a heterogeneous wetland site with the help of URANOS. The model also supported theoretical considerations by Schrön et al. (2021) about the influence of train wagons, rail tracks and shielding material on a CRNS detector in moving trains, which then have been confirmed by experiments. Köhli and Schmoldt (2022) investigated the possibilities of detecting unexploded ordnance using the pulsed neutron–neutron technique with the operation mode “URANOS Downhole”.

Table 3 The standard setup for a layer composition in cosmic neutron sensing.

9 Performance benchmarks

The performance of the code heavily depends on the simulation setup. The most significant contributor is the mean lifetime of neutrons in terms of scatterings within the domain. In simple configurations like for the analysis of detectors most neutrons undergo only a few interactions before either being absorbed or leaving the domain. Atmospheric neutrons can have up to hundreds of scatterings before ending up thermalized. Additionally, the performance depends on the chosen materials and energy range, as for elastic scattering or absorption only a few cross sections are evaluated. In the MeV domain, however, adding up all inelastic channels scales up to dozens of address requests. In order to provide practical estimations a standard setup can be defined as in Table 3. The domain measures 900 m × 900 m with a source dimension of 840 m × 840 m. It contains the minimum configuration of six layers.

A similar setup has also been used for simulating the so-called “UFZ site” in Schrön et al. (2018) – an urban environment with many concrete buildings, streets, green spaces, railroad lines, a lake and trees. Of the 12 layers, a total of 8 contained a matrix of 1800 × 1800 pixels. In the provided minimum configuration, the source is positioned relatively close to the ground. In case atmospheric effects like the influence of air humidity needs to be studied, the source should be located at an altitude of 500–600 m.

Table 4 Benchmark results for the single core performance of URANOS (v0.99ρ) and MCNP for a number of practically relevant scenarios. URANOS was tested on a 4 GHz i7-6700K Skylake CPU in a Windows 10 environment, MCNP6 on a 2.7 GHz E5-2680 Sandy Bridge CPU in Windows Server 2016 and GEANT4 (10.7) on a 3.4 GHz E-2124G Coffee Lake CPU in Debian 10.

In order to provide references for the computational speed in realistic scenarios the following Table 4 summarizes the single core performance of the code in terms of neutrons per second per GHz per core. The standard setup and UFZ site are described above and are simulated in combination with the cosmic neutron spectrum like presented in Fig. 15; the detector is a rover-type instrument (Köhli et al., 2018b), which is a setup similar to the Bonner sphere models, and the other benchmarks are synthetic. Without additional voxel geometry descriptions URANOS requires approximately 230 MB of memory, mainly for storing ENDF data. The benchmark results show that the performance of URANOS scales proportionally to the number of underlying calculations, i.e., scatterings. Applying a voxel geometry definition instead of a homogeneous layer has only a minor influence on the computational speed as long as the whole model can be kept in the system memory. Voxel models do not require collision tests for every surface within the domain; the number of intersection tests scales only by the voxel density. Compared to URANOS MCNP6 is less susceptible to the environmental topology or soil water content. From dry to moist conditions there is roughly a difference of 15 %, and another 15 % can be gained by deactivating thermal neutron transport.

URANOS can be compiled under Windows and Linux/OSX, with the latter not specifically considered in the performance evaluation. In general, with respect to implementations on clusters, running several instances of URANOS scales linearly with the number of physical CPU cores.

10 User interface

10.1 Creation of input files

The layered geometry of URANOS allows for building the geometry as a stack of materials. Especially some layers can be described by a single homogeneous material, and others can contain complex patters constructed on the basis of images or ASCII matrices. This allows us to easily realize for example a model domain composed of a soil with structures (e.g., buildings) erected on top. The procedure is explained in Sect. 4.2, and some examples for input image files are shown in Fig. 4. Beyond the material input definition each layer can contain a density multiplicator matrix, which scales the density of the respective voxel by a certain factor and a porosity matrix which allows for adapting the soil bulk density. Images or matrices are stored in the active working folder and have to be labeled as the layer they refer to, for example 6.dat or 6.png define the geometry in layer 6. 6d.dat would be the density matrix and 6p.dat the porosity matrix. In order to create input image files, it is recommended to use a software which works on a pixel level. Contrasts must not be smeared out as usually done in photography software. Grayscale values between two colors represent completely different materials.

Figure 16 URANOS main user interface with (1) main simulation control bar, (2) configuration tabs, (3) Live View tabs, (4) global environmental settings, (5) geometry stack with two layers defined by voxels, (6) Bird's eye View of the neutron flux within the domain in the detector layer and (7) cosmic-ray neutron spectrum above the ground.

10.2 Command line

URANOS can be run with different startup parameters. One of the options is to run the simulation without the GUI using a config file. This enables us to create a batch file to start several jobs at once. The arguments that can be used in combination with UranosGUI.exe are the following:

• silent. Disables console output.

• noGUI location/of/uranos.cfg. Starts the simulation without the GUI and takes all configuration options necessary for the model run from the uranos.cfg. The file also contains input and output directories. If no file is specified the configuration from the active directory is taken.

• batchrun i j. Conducts a batch run with 11 × 11 = 121 values for soil moisture (1, 3, 6, 10, 12, 15, 20, 25, 30, 40, 50) in volumetric percent and air humidity (1, 2, 4, 7, 10, 12, 15, 18, 21, 27, 35) in g m⁻³. These values are applied to the materials air and soil. i and j specify numbers for the setting to start and end with.

• detectorBatchRun. Conducts a batchrun with 40 different energies from 0 eV to 20 MeV released orthogonally onto the domain (ϑ=0). The domain is supposed to contain a detector model. An additional output file is created for the number of neutrons which are absorbed by helium-3 and boron-10.

Figure 17 URANOS user interface showing the detector tabs. (left) Settings for changing energy and angular sensitivity, as well as geometric shape and extension. (right) Detector scoring output showing the origins of detected neutrons.

10.3 GUI

The design of the graphical user interface provides all settings for adjusting the simulation alongside with the data output and information which are relevant to ensure a proper configuration at runtime. The main window is displayed in Figs. 16 and 17. The control bar (1) with its functions to start and stop the model run shows the simulation progress and the expected time to complete the simulation. The main user area is split into the settings pages (2) and the Live View tabs (3). The sliders (4) control general features of air and soil, which can be used as materials in the geometry stack (5). After starting the simulation the neutron distribution within the domain can be viewed in the visualization tab (6). The cosmic-ray neutron spectrum above the ground provides a quick overview of the scattering and absorption characteristics of the domain.

In the left-side settings pages the available tabs contain a variety of configuration options, such as:

• Physical parameters. In the main control tab the central table is used to build the geometry stack by defining the height and thickness, as well as the base material, for each layer. Matrix definitions for materials, densities and porosity are automatically loaded and assigned if the respective option is used. The sliders on the left side control the general settings for the moisture content of the materials air and soil. For the incoming spectrum the vertical cutoff rigidity can be chosen.

• Computational parameters. The lateral size of the domain and source can be adjusted in this tab alongside with specific options like boundary conditions or the activation of calculation models like thermal neutron transport. These options largely influence the computational speed and therefore need to be adjusted to fit the research question.

• Detector. In URANOS the settings for the detector and detector layer can be chosen independently. For both entities a flat detection efficiency within upper and lower energy bounds can be chosen or the physics model, which emulates the response function of an actual CRNS instrument. By using an enlarged virtual detector the simulation can gather much more statistics than using a detector model with its actual dimensions, as this significantly enlarged representation can score neutrons using the same characteristics as a realistic model. Users can create their own response functions and load them into the simulation.

• Showcase. The cosmic-ray input spectrum as well as the footprint range can be viewed using this tab.

• Folders. This tab allows one to set input and output folders for URANOS.

• Export. The type of individual data sets to be exported from the detector and detector layer recordings can be selected.

• Display. This tab contains the controls for the Live View graphs, energy range, color scheme and scaling. For the Bird's eye View either hit or track distribution are available.

The right-side visualization elements are the following:

• Bird's eye View and Spectra. The main tab shows the intersections of neutrons with the detector layer as a representation of track or hit density. The whole layer can act as a virtual detector by applying the respective response function. The + button opens a window with an enlarged high-resolution output of the currently shown detector layer display. The lower graph shows the energy distribution of neutrons above the ground. The full spectrum can be divided into an incoming part, which has been used as the initial spectrum, and an albedo part, which represents neutrons having at least one soil contact.

• Range view. In this tab the range distributions for detector and detector layer are presented which are also called horizontal weighting functions. The range is defined as the distance from the first collision in the soil to detection.

• Spatial view. The top graph shows a horizontal cut view through the innermost 8 % of the detector layer for the chosen energy range. Furthermore in the depth distribution three plots are shown: the interaction of all neutrons within the soil; the maximum probe depth, which is the lowest point a neutron reached when being registered by the detector layer; and the last probe depth, which shows the depth distribution of interaction points before leaving the soil and being detected. Both distributions can be regarded as approximations for the minimum and maximum depths for probing the soil.

• Detector. This tab shows the points of first soil contact for neutrons scored by the virtual detector. By setting the range of interest in the display tab to a value different from none, the plot changes from showing hits to the spatial density of origins.

10.4 Output formats

URANOS provides three output channels for the simulation data. The neutron flux distribution within the detector layer can be exported as an ASCII matrix or image. This data do only provide spatial information. The virtual detector can export the full information of neutrons passing through it: direction, energy, last interaction, first soil contact or model-dependent detection probabilities. Additionally, the full history, i.e., the complete neutron track list, can be exported. This feature allows us to back-trace only those neutrons which match specific detection criteria. All simulation data are furthermore stored in two compressed ROOT files; one contains all graphs and distributions seen in the GUI itself, and one contains support information mainly addressing mechanisms inside the domain like element- and materialwise scattering and absorption distributions.

11 Conclusions

URANOS is a new neutron Monte Carlo tool based on C++ with a graphical user interface specifically adapted to the needs of environmental physics. It uniquely features a modeling geometry using the concept of voxels, three-dimensional pixels, which can be stacked in layers and extruded from grayscale images or ASCII matrices. Neither a three-dimensional modeling software for vector objects nor the writing of elaborate steering files is necessary in order to realize complex landscapes and domain structures. A built-in, validated cosmic-ray neutron spectrum leads to quick solutions, as it makes the extensive particle shower generation redundant. It allows us to record neutron flux by a domain-wide scoring layer and a virtual detector with preconfigured instrument characteristics. These built-in tools shorten the data analysis for modeling results without requiring programming skills. URANOS features the computation of all relevant neutron interactions below 20 MeV. As for higher energies particles other than neutrons contribute to their production; URANOS uses an effective cascade model which reproduces the empirically know flux by adjustable parameters. Comparisons show a close agreement of the URANOS model with other simulation toolkits like MCNP. URANOS is now available for all users and regularly maintained to address the growing needs of the community.

Appendix A: Elements, isotopes and reaction types

The database of URANOS materials relies on a library of predefined elements. Such are described by ENDF cards, which are extracted from the existing sources mentioned in Sect. 2.5 and stored individually. The following Table A1 is a comprehensive list of isotopes, which have been selected and implemented.

Table A1 Available isotopes in URANOS and cross sections used, identifiers according to Trkov et al. (2012).

Appendix B: Material codepages

URANOS provides a list of predefined materials, which are combinations of elements described in Appendix A. Table B1 summarizes all available compositions which are implemented as materials.

Table B1 List of preconfigured materials available in URANOS with their composition and density.