the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
URANOS v1.0 – the Ultra Rapid Adaptable NeutronOnly Simulation for Environmental Research
Markus Köhli
Martin Schrön
Steffen Zacharias
Ulrich Schmidt
The understanding of neutron transport by Monte Carlo simulations led to major advancements towards precise interpretation of measurements. URANOS (Ultra Rapid NeutronOnly 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 cosmicray 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 layerwise, whereas in each layer a voxel geometry is extruded using a twodimensional map from pixel images representing predefined materials and allowing for the construction of objects on the basis of pixel graphics without a threedimensional editor. It furthermore features predefined cosmicray 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 cosmicray neutron footprint, signals in complex geometries like mobile applications on roads, urban environments and snow patterns.
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 socalled 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 NParticle 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 cosmicray neutron sensing, we developed the Monte Carlo code URANOS (Ultra Rapid Adaptable NeutronOnly 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) voxelbased 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 raycasting 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/BVII.1 (Chadwick et al., 2011), ENDF/BVIII.0 (Brown et al., 2018) and the JENDL/HE2007 (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 crossplatform framework.
The graphical user interface offers features specifically tailored to the needs of the field of cosmicray neutron sensing. This novel method retrieves subsurface soil moisture by measuring flux of cosmicrayinduced 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 bruteforce calculation technique, which is used for complex problems consisting of welldefined 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 generalpurpose 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 mid90s 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 (X5 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 HighEnergy 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 CascadeExciton Model (CEM) (Gudima et al., 1983) and the Los Alamos QuarkGluon 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 Xbranch into the main development branch. It provides an optional cosmicray 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 highenergy 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 SlowingDown Approximation) (Nelms, 1956) and intranuclear 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 followup developments is the PARMA (PHITSbased 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 highenergy 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 lowenergy 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 cosmicray neutron studies (Liu et al., 2021).
1.2 Why another code?
The choice for creating an own independently operating Monte Carlobased 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 nuclearrelated technology – whereas the underlying databases are free to access. Highprecision 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 timeconsuming code tuning. At best they received wrappers in C, rarely in C
++
. Today, facing multithreading, 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 boronbased hybrid detector requires two additional steps of chargedparticle 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^{5} due to the low thickness are computationally difficult to model.

All available codes propagate a takeoff 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 lowenergy 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 JEFF3.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, objectoriented language.
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 pseudorandom 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}\phantom{\rule{0.125em}{0ex}}\mathrm{937}}\mathrm{1}\approx \mathrm{4.3}\times {\mathrm{10}}^{\mathrm{6001}}$, low correlation between subsequent numbers (kdistributed 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
with the macroscopic cross section Σ_{t}, which in general is energydependent. Solutions of this type of differential equation are exponential functions. For the noninteraction probability one therefore can write
The probability distribution function for the distance to the next collision (2) assuming conditional probabilities transforms to
The free path length l is obtained by the cumulative probability distribution function of Eq. (3) by
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:
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 crosssection 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
one has to conserve the reaction rate, the integrand of Eq. (6),
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 centerofmass (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\mathrm{cos}\mathit{\vartheta}}$. Such a probability function can be constructed by
Defining the denominator of Eq. (8) as the normalization factor C and
as well as μ=cos ϑ, one obtains
In order to obtain a sampling scheme one can divide Eq. (10) into two parts such that
Here the reason for dividing and multiplying Eq. (10) by v+V is that g_{1} 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_{1}(x)g_{2}(x) can be sampled by sampling x^{′} from a normalized distribution q(x):
and accepting it with a probability of
with g_{1}(x) bounded. In order to determine q(V) it is necessary to integrate g_{2} into Eq. (11):
leading to sampling the probability distribution function
By substituting x=βV, likewise dx=βdV, and y=βv leads finally to
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 ξ_{1} is sampled from [0,1) and if
the function 2x^{3}exp (−x^{2}) is sampled, otherwise the function $\mathrm{4}/\sqrt{\mathit{\pi}}{x}^{\mathrm{2}}\mathrm{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 ξ_{2} in [0,1] by
and as the maximum of g_{1} is $\mathrm{4}\mathit{\sigma}\left({v}_{r}\right)/\sqrt{\mathit{\pi}}{C}^{\prime}$ another sampling random number ξ_{3} can be used to accept speed and angle by
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 socalled 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 highenergy branch JENDLHE 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).
The design of URANOS was motivated by the following general aspects:

The geometry is represented in a threedimensional coordinate space with dimensions from the centimeter to the kilometer scale.

In typical model runs the number of neutrons can easily reach 10^{9} 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 highenergy 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 substructured by twodimensional 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.
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 raycasting technique (Roth, 1982) refers to conducting a series of raysurface 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 Xray 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.
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.
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 twodimensional 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 ^{1}H, ^{10}B and ^{16}O, whereas the full list of available isotopes can be found in Appendix A.
Only MT numbers with significant contributions are taken into account, which translates to omitting processes with overall less than 10^{−2} % 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 nonskipped 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.
For calculating the total macroscopic cross section the individual contributions of elastic Σ^{e} and inelastic Σ^{in} scattering, as well as absorption Σ^{a}, are summed up
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 chargedparticle ejection by converters.
4.2 Geometry
URANOS uses analytical geometry definitions and voxels as introduced in Sect. 3. The following topdown structure is applied for describing the simulation environment:
Each layer of the stack is either entirely composed of a material or subdivided into several sections using a twodimensional 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).
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 tabseparated 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:
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.
Neutron tracks S are described by a mixed geometry definition of support vectors x in Cartesian coordinates and spherical direction vectors r:
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
and for determining the position on a layer at elevation z_{L}
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.
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:
For source definitions on a logarithmic support in [10^{a},10^{b}], like in Sect. 5.1, ξ is scaled to
5.1 The cosmic neutron source
The cosmic neutron source definition is specifically designed for the problem of soilmoisturedependent 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 cosmicray 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 downwardoriented component of the flux, which is used as the default cosmic neutron spectrum. Figure 5 exemplarily shows the URANOS cosmicray 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
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}}}}\mathrm{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\phantom{\rule{0.125em}{0ex}}\mathrm{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 semiempirical description for fission neutrons is the Watt spectrum (Watt, 1952), especially used for ^{235}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}\phantom{\rule{0.125em}{0ex}}\mathrm{sinh}\left(\sqrt{\mathrm{2}E}\right)\mathrm{exp}(E)\end{array}$$and for ^{252}Cf (Smith et al., 1957)
$$\begin{array}{}\text{(27)}& N\left(E\right)=\mathrm{sinh}\left(\sqrt{\mathrm{2}E}\right)\mathrm{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{\mathrm{exp}\left(ab/\mathrm{4}\right)}{\mathrm{0.5}\phantom{\rule{0.125em}{0ex}}\sqrt{\mathit{\pi}{a}^{\mathrm{3}}b}}}\phantom{\rule{0.125em}{0ex}}\mathrm{sinh}\left(\sqrt{bE}\right)\mathrm{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.
6.1 Loop nodes
The main calculation routine runs in two loops; for each neutron in the neutron stack, the geometry stack is traversed layerbylayer:
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 isotopebyisotope. 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 ^{10}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$\phantom{\rule{0.125em}{0ex}}\mathrm{keV}<E<\phantom{\rule{0.125em}{0ex}}$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
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=\mathrm{ln}\left(\mathit{\xi}\right)/{\mathrm{\Sigma}}_{\mathrm{t}}$. An example of the total crosssection 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_{0} and the layer z position z_{l} and height d_{l}.
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 zprojected length z_{m} of one lateral unit pixel s_{p} for the actual direction vector by ${z}_{\mathrm{m}}={s}_{\mathrm{p}}/\mathrm{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_{0}±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 voxelbased 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
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:
as the probability of selecting a reaction type i from all possible interaction channels is
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}:^{1}
For inelastic scattering additionally the excited state of the target isotope has to be determined. One vector
of numbers contains the cumulative crosssection 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
then the corresponding isotope is taken from v_{e}[i].^{2} Figure 8 schematically illustrates the calculation routine for all processes described above.
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):
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:
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
and added to the existing direction vectors
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 lowZ 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{\u203e}}{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 ^{235}U produces on average $\approx \mathrm{2.4}+E$ MeV^{−1} neutrons per fission.
6.3.3 Absorption
Neutrons are either absorbed by a nonradiating 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 highenergy 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^{2} g^{−1}. Experimental values for L_{prim} are in the range of (135–155) cm^{2} g^{−1}, depending on the latitude of the site. Here, the value from Sato (2016) is taken.
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.
Scoring options for CRNS
In the most simple case a uniform detection efficiency ϵ can be chosen within a specific range:
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 cosmicray 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
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.
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 cosmicray 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 nontrivial 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 lessdense 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 highenergy 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.
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 slowingdown 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 nonequilibrium 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.
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 energydependent comparison and on the interplay of moderator and absorber. Among the various existing technical realizations the heliumbased 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^{∘}.
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.
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 heightdependent scaling of the cosmicray 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 highenergy 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 energydependent 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 highenergy 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.
On the qualitative level the underlying physics model correctly calculates the response to the soil. Three peaks can be observed – the highenergy domain around 100 MeV with the incomingonly 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 highenergy cascade contributes to the production of neutrons by three different particle species with three different fractions and attenuation lengths (around 140 cm^{2} g^{−1} for neutrons, around 110 cm^{2} g^{−1} for protons and around 500 cm^{2} g^{−1} for muons). Future work will address these different attenuation lengths and model the highenergy 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 aboveground 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 realworld 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 boron10based 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 finetune the event reconstruction in the BODELAIRE (BOronbased 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 cosmicray 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 downwardshielded 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 directionalsensing 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”.
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 socalled “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.
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 rovertype 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.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.
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 theuranos.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^{−3}. These values are applied to the materials air and soil.i
andj
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 helium3 and boron10.
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 cosmicray neutron spectrum above the ground provides a quick overview of the scattering and absorption characteristics of the domain.
In the leftside 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 cosmicray 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 rightside 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 highresolution 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 modeldependent detection probabilities. Additionally, the full history, i.e., the complete neutron track list, can be exported. This feature allows us to backtrace 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.
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, threedimensional pixels, which can be stacked in layers and extruded from grayscale images or ASCII matrices. Neither a threedimensional modeling software for vector objects nor the writing of elaborate steering files is necessary in order to realize complex landscapes and domain structures. A builtin, validated cosmicray neutron spectrum leads to quick solutions, as it makes the extensive particle shower generation redundant. It allows us to record neutron flux by a domainwide scoring layer and a virtual detector with preconfigured instrument characteristics. These builtin 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.
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.
The URANOS source code is made available at the GitHub repository https://github.com/mkoehli/uranos/tree/URANOS (last access: 18 January 2023). URANOS v1.0 has been released under DOI https://doi.org/10.7910/DVN/THPNZW (Köhli, 2022a). Furthermore the code has been released, including a collection of examples and user guides under DOI https://doi.org/10.5281/zenodo.6578668 (Köhli, 2022b). This DOI represents all versions and will always resolve to the latest one.
Libraries and data used in this publication have been released in the abovementioned GitHub repository and are available from the Harvard Dataverse archive at DOI https://doi.org/10.7910/DVN/THPNZW (Köhli, 2022a).
The current model of URANOS is also available from the project website.
URANOS v1.0 is linked against ROOT 6.22.08, QT 5.15 and QCustomPlot 2.1.0.
MK and MS wrote the manuscript. SZ and US edited the manuscript. SZ and US supervised the development of URANOS. MK designed the structure of URANOS and wrote and implemented the physics core of URANOS. MK and MS designed the graphical user interface and tested the software.
Markus Köhli holds a CEO position at StyX Neutronica GmbH, Germany.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Markus Köhli acknowledges all URANOS users who have contributed to the development of URANOS: Jannis Weimar, Paul Schattan, Daniel Rasche, Jannis Jakobi, Cosimo Brogi, Sugan Kanagasenthinathan, Patrick Stowell and Tatsuhiko Sato. Martin Schrön and Markus Köhli further acknowledge Peter Dietrich for funding and continuous support of the model development. URANOS was developed within several projects. Initial funding was provided by the projects Neutron Detectors for the MIEZE method and Forschung und Entwicklung hochauflösender Neutronendetektoren, funded by the German Federal Ministry of Education and Research (BMBF) (grant identifiers: 05K10VHA and 05K16PD1) and by the German Research Foundation (DFG) research unit FOR 2694 Cosmic Sense via the project 357874777.
This research has been supported by the Bundesministerium für Bildung und Forschung (grant nos. 05K10VHA and 05K16PD1) and the Deutsche Forschungsgemeinschaft (grant no. 357874777).
This paper was edited by Wolfgang Kurtz and reviewed by two anonymous referees.
Agostinelli, S., Allison, J., Amako, K., et al.: GEANT4 – a simulation toolkit, Nucl. Instrum. Meth. A, 506, 250–303, https://doi.org/10.1016/S01689002(03)013688, 2003. a
Andreasen, M., Jensen, H. K., Zreda, M., Desilets, D., Bogena, H., and Looms, C.: Modeling cosmic ray neutron field measurements, Water Resour. Res., 52, 6451–6471, https://doi.org/10.1002/2015wr018236, 2016. a
Badiee, A., Wallbank, J., Pulido Fentanes, J., Trill, E., Scarlet, P., Zhu, Y., Cielniak, G., Cooper, H., Blake, J., Evans, J., Zreda, M., Köhli, M., and Pearson, S.: Using Additional Moderator to Control the Footprint of a COSMOS Rover for Soil Moisture Measurement, Water Resour. Res., 57, e2020WR028478, https://doi.org/10.1029/2020wr028478, 2021. a, b
Baroni, G., Scheiffele, L., Schrön, M., Ingwersen, J., and Oswald, S.: Uncertainty, sensitivity and improvements in soil moisture estimation with cosmicray neutron sensing, J. Hydrol., 564, 873–887, https://doi.org/10.1016/j.jhydrol.2018.07.053, 2018. a
Battistoni, G., Boehlen, T., Cerutti, F., Chin, P., Esposito, L., Fasso, A., Ferrari, A., Lechner, A., Empl, A., Mairani, A., Mereghetti, A., Ortega, P., Ranft, J., Roesler, S., Sala, P., Vlachoudis, V., and Smirnov, G.: Overview of the FLUKA code, Ann. Nucl. Energ., 82, 10–18, https://doi.org/10.1016/j.anucene.2014.11.007, 2015. a
Boudard, A., Cugnon, J., David, J.C., Leray, S., and Mancusi, D.: New potentialities of the Liège intranuclear cascade model for reactions induced by nucleons and light charged particles, Phys. Rev. C, 87, 014606, https://doi.org/10.1103/PhysRevC.87.014606, 2013. a
Bramblett, R. and Bonner, T.: Neutron evaporation spectra from (p, n) reactions, Nucl. Phys., 20, 395–407, https://doi.org/10.1016/00295582(60)901826, 1960. a
Bramblett, R., Ewing, R., and Bonner, T.: A new type of neutron spectrometer, Nucl. Instrum. Methods, 9, 1–12, https://doi.org/10.1016/0029554X(60)900434, 1960. a
Briesmeister, J.: MCNPA general Monte Carlo Nparticle transport code, Version 4C, LA13709M, https://permalink.lanl.gov/object/tr?what=info:lanlrepo/lareport/LA13709M, 2000. a
Brown, D., Chadwick, M., Capote, R., Kahler, A., Trkov, A., Herman, M., Sonzogni, A., Danon, Y., Carlson, A., Dunn, M., Smith, D., Hale, G., Arbanas, G., Arcilla, R., Bates, C., Beck, B., Becker, B., Brown, F., Casperson, R., Conlin, J., Cullen, D., Descalle, M.A., Firestone, R., Gaines, T., Guber, K., Hawari, A., Holmes, J., Johnson, T., Kawano, T., Kiedrowski, B., Koning, A., Kopecky, S., Leal, L., Lestone, J., Lubitz, C., Damiá¡n, J. M., Mattoon, C., McCutchan, E., Mughabghab, S., Navratil, P., Neudecker, D., Nobre, G., Noguere, G., Paris, M., Pigni, M., Plompen, A., Pritychenko, B., Pronyaev, V., Roubtsov, D., Rochman, D., Romano, P., Schillebeeckx, P., Simakov, S., Sin, M., Sirakov, I., Sleaford, B., Sobes, V., Soukhovitskii, E., Stetcu, I., Talou, P., Thompson, I., van der Marck, S., WelserSherrill, L., Wiarda, D., White, M., Wormald, J., Wright, R., Zerkle, M., Žerovnik, G., and Zhu, Y.: ENDF/BVIII.0: The 8th Major Release of the Nuclear Reaction Data Library with CIELOproject Cross Sections, New Standards and Thermal Scattering Data, Nucl. Data Sheets, 148, 1–142, https://doi.org/10.1016/j.nds.2018.02.001, 2018. a
Brun, R. and Rademakers, F.: ROOT – An Object Oriented Data Analysis Framework, in: Proceedings of AIHENP’96 Workshop, Lausanne, vol. 389, 81–86, https://doi.org/10.1016/S01689002(97)00048X, 1997. a, b
Caswell, R., Gabbard, R., Padgett, D., and Doering, W.: Attenuation of 14.1Mev Neutrons in Water, Nucl. Sci. Eng., 2, 143–159, https://doi.org/10.13182/NSE57A25383, 1957. a, b, c, d
Chadwick, M., Herman, M., Obložinský, P., Dunn, M., Danon, Y., Kahler, A., Smith, D., Pritychenko, B., Arbanas, G., Arcilla, R., Brewer, R., Brown, D., Capote, R., Carlson, A., Cho, Y., Derrien, H., Guber, K., Hale, G., Hoblit, S., Holloway, S., Johnson, T., Kawano, T., Kiedrowski, B., Kim, H., Kunieda, S., Larson, N., Leal, L., Lestone, J., Little, R., McCutchan, E., MacFarlane, R., MacInnes, M., Mattoon, C., McKnight, R., Mughabghab, S., Nobre, G., Palmiotti, G., Palumbo, A., Pigni, M., Pronyaev, V., Sayer, R., Sonzogni, A., Summers, N., Talou, P., Thompson, I., Trkov, A., Vogt, R., van der Marck, S., Wallner, A., White, M., Wiarda, D., and Young, P.: ENDF/BVII.1 Nuclear Data for Science and Technology: Cross Sections, Covariances, Fission Product Yields and Decay Data, Nucl. Data Sheets, 112, 2887–2996, https://doi.org/10.1016/j.nds.2011.11.002, 2011. a
Cnudde, V. and Boone, M.: Highresolution Xray computed tomography in geosciences: A review of the current technology and applications, EarthSci. Rev., 123, 1–17, https://doi.org/10.1016/j.earscirev.2013.04.003, 2013. a
Cugnon, J., Volant, C., and Vuillier, S.: Nucleon and deuteron induced spallation reactions, Nucl. Phys. A, 625, 729–757, https://doi.org/10.1016/S03759474(97)006027, 1997. a
Decker, A., McHale, S., Shannon, M., Clinton, J., and McClory, J.: Novel Bonner Sphere Spectrometer Response Functions Using MCNP6, IEEE T. Nucl. Sci., 62, 1689–1694, https://doi.org/10.1109/TNS.2015.2416652, 2015. a
Desilets, D. and Zreda, M.: Footprint diameter for a cosmicray soil moisture probe: Theory and Monte Carlo simulations, Water Resour. Res., 49, 3566–3575, https://doi.org/10.1002/wrcr.20187, 2013. a
Desilets, D., Zreda, M., and Prabu, T.: Extended scaling factors for in situ cosmogenic nuclides: new measurements at low latitude, Earth Planet. Sc. Lett., 246, 265–276, https://doi.org/10.1016/j.epsl.2006.03.051, 2006. a
Emmett, M.: The MORSE Monte Carlo Transport Code System, ORNL4972/R2, https://inis.iaea.org/collection/NCLCollectionStore/_Public/16/029/16029296.pdf, 1975. a
Federico, C., Gonçalez, O., Fonseca, E., Martin, I., and Caldas, L.: Neutron spectra measurements in the south Atlantic anomaly region, Radiat. Meas., 45, 1526–1528, https://doi.org/10.1016/j.radmeas.2010.06.038, 2010. a
Francke, T., Heistermann, M., Köhli, M., Budach, C., Schrön, M., and Oswald, S. E.: Assessing the feasibility of a directional cosmicray neutron sensing sensor for estimating soil moisture, Geosci. Instrum. Method. Data Syst., 11, 75–92, https://doi.org/10.5194/gi11752022, 2022. a, b
Franz, T. E., Zreda, M., Rosolem, R., and Ferre, T. P. A.: A universal calibration function for determination of soil moisture with cosmicray neutrons, Hydrol. Earth Syst. Sci., 17, 453–460, https://doi.org/10.5194/hess174532013, 2013. a
Galassi, M., Davies, J., Theiler, J., Gough, B., Jungman, G., Alken, P., Booth, M., Rossi, F., and Ulerich, R.: GNU Scientific Library: reference manual for GSL version 2.7, Free Software Foundation, https://www.gnu.org/software/gsl/manual/gslref.pdf (last access: 18 January 2023), 2016. a
Garny, S., Mares, V., and Rühm, W.: Response functions of a Bonner Sphere Spectrometer calculated with GEANT4, Nucl. Instrum. Meth. A, 604, 612–617, https://doi.org/10.1016/j.nima.2009.02.044, 2009. a
Gelbard, E.: Epithermal scattering in VIM, Tech. Rep. FRATM123, Argonne National Laboratory, IL, USA, technical Memorandum, 1979. a
Goldhagen, P., Reginatto, M., Kniss, T., Wilson, J., Singleterry, R., Jones, I., and Van Steveninck, W.: Measurement of the energy spectrum of cosmicray induced neutrons aboard an ER2 highaltitude airplane, Nucl. Instrum. Meth. A, 476, 42–51, https://doi.org/10.1016/S01689002(01)013869, 2002. a
Goldman, L.: Principles of CT and CT Technology, J. Nucl. Med. Tech., 35, 115–128, https://doi.org/10.2967/jnmt.107.042978, 2007. a
Goorley, T., James, M., Booth, T., Brown, F., Bull, J., Cox, L., Durkee, J., Elson, J., Fensin, M., Forster, R., Hendricks, J., Hughes, H., Johns, R., Kiedrowski, B., Martz, R., Mashnik, S., McKinney, G., Pelowitz, D., Prael, R., Sweezy, J., Waters, L., Wilcox, T., and Zukaitis, T.: Initial MCNP6 Release Overview, Nucl. Tech., 180, 298–315, https://doi.org/10.13182/NT11135, 2012. a, b
Gudima, K., Mashnik, S., and Toneev, V.: Cascadeexciton model of nuclear reactions, Nucl. Phys. A, 401, 329–361, https://doi.org/10.1016/03759474(83)905328, 1983. a
Gudima, K., Mashnik, S., and Sierk, A.: User Manual for the Code LAQGSM, LAUR016804, https://permalink.lanl.gov/object/tr?what=info:lanlrepo/lareport/LAUR016804, 2001. a
Hébert, A.: Multigroup Neutron Transport and Diffusion Computations, Springer US, Boston, 751–911, https://doi.org/10.1007/9780387981499_8, 2010. a
Heidbüchel, I., Güntner, A., and Blume, T.: Use of cosmicray neutron sensors for soil moisture monitoring in forests, Hydrol. Earth Syst. Sci., 20, 1269–1288, https://doi.org/10.5194/hess2012692016, 2016. a
Hertel, N. and Davidson, J.: The response of Bonner Spheres to neutrons from thermal energies to 17.3 MeV, Nucl. Instrum. Meth. A, 238, 509–516, https://doi.org/10.1016/01689002(85)904942, 1985. a
ISO group 85/SC 2 Radiological protection: Reference neutron radiations – Part 1, Standard, International Organization for Standardization, Geneva, CH, https://www.iso.org/standard/25666.html (last access: 18 January 2023), 2001. a
Iwamoto, O., Iwamoto, N., Kunieda, S., Minato, F., and Shibata, K.: The CCONE Code System and its Application to Nuclear Data Evaluation for Fission and Other Reactions, Nucl. Data Sheets, 131, 259–288, https://doi.org/10.1016/j.nds.2015.12.004, 2016. a
Iwase, H., Niita, K., and Nakamura, T.: Development of GeneralPurpose Particle and Heavy Ion Transport Monte Carlo Code, J. Nucl. Sci. Technol., 39, 1142–1151, https://doi.org/10.1080/18811248.2002.9715305, 2002. a, b
Iwema, J., Schrön, M., Da Silva, J. K., Schweiser De Paiva Lopes, R., and Rosolem, R.: Accuracy and precision of the cosmicray neutron sensor for soil moisture estimation at humid environments, Hydrol. Process., 35, e14419, https://doi.org/10.1002/hyp.14419, 2021. a
Iyer, M. and Ganguly, A.: Neutron Evaporation and Energy Distribution in Individual Fission Fragments, Phys. Rev. C, 5, 1410–1421, https://doi.org/10.1103/PhysRevC.5.1410, 1972. a
Jakobi, J., Huisman, J., Köhli, M., Rasche, D., Vereecken, H., and Bogena, H.: The Footprint Characteristics of Cosmic Ray Thermal Neutrons, Geophys. Res. Lett., 48, e2021GL094281, https://doi.org/10.1029/2021gl094281, 2021. a, b
Kawano, T., Talou, P., Stetcu, I., and Chadwick, M.: Statistical and evaporation models for the neutron emission energy spectrum in the centerofmass system from fission fragments, Nucl. Phys. A, 913, 51–70, https://doi.org/10.1016/j.nuclphysa.2013.05.020, 2013. a
Kittelmann, T. and Boin, M.: Polycrystalline neutron scattering for Geant4: NXSG4, Comput. Phys. Commun., 189, 114–118, https://doi.org/10.1016/j.cpc.2014.11.009, 2015. a
Knuth, D.: The Art of Computer Programming, Volume 2, 3rd Edn., Seminumerical Algorithms, AddisonWesley Longman Publishing Co., Inc., Boston, MA, USA, ISBN 9780321635778, 1997. a
Köhli, M.: URANOS v1.0, Harvard Dataverse V1 [code], https://doi.org/10.7910/DVN/THPNZW, 2022a. a, b
Köhli, M.: mkoehli/uranos: URANOS v1.0 (URANOS), Zenodo [data set], https://doi.org/10.5281/zenodo.6578668, 2022b. a
Köhli, M. and Schmoldt, J.P.: Feasibility of UXO detection via pulsed neutronneutron logging, Appl. Radiat. Isotopes, 188, 110403, https://doi.org/10.1016/j.apradiso.2022.110403, 2022. a, b
Köhli, M., Schrön, M., Zreda, M., Schmidt, U., Dietrich, P., and Zacharias, S.: Footprint characteristics revised for fieldscale soil moisture monitoring with cosmicray neutrons, Water Resour. Res., 51, 5772–5790, https://doi.org/10.1002/2015WR017169, 2015. a, b, c
Köhli, M., Allmendinger, F., Häußler, W., Schröder, T., Klein, M., Meven, M., and Schmidt, U.: Efficiency and spatial resolution of the CASCADE thermal neutron detector, Nucl. Instrum. Meth. A, 828, 242–249, https://doi.org/10.1016/j.nima.2016.05.014, 2016. a, b
Köhli, M., Desch, K., Gruber, M., Kaminski, J., Schmidt, F., and Wagner, T.: Novel neutron detectors based on the time projection method, Physica B, 551, 517–522, https://doi.org/10.1016/j.physb.2018.03.026, 2018a. a, b
Köhli, M., Schrön, M., and Schmidt, U.: Response functions for detectors in cosmic ray neutron sensing, Nucl. Instrum. Meth. A, 902, 184–189, https://doi.org/10.1016/j.nima.2018.06.052, 2018b. a, b, c
Köhli, M., Weimar, J., Schrön, M., Baatz, R., and Schmidt, U.: Soil Moisture and Air Humidity Dependence of the AboveGround CosmicRay Neutron Intensity, Front. Water, 2, 15, https://doi.org/10.3389/frwa.2020.544847, 2021. a, b, c
Koning, A., Bauge, E., Dean, C., Dupont, E., Nordborg, C., Rugama, Y., Fischer, U., Forrest, R., Kellett, M., Jacqmin, R., Leeb, H., Mills, R., Pescarini, M., and Rullhusen, P.: Status of the JEFF Nuclear Data Library, J. Korean Phys. Soc., 59, 1057–1062, https://doi.org/10.3938/jkps.59.1057, 2011. a
Lefmann, K. and Nielsen, K.: McStas, a general software package for neutron raytracing simulations, Neutron News, 10, 20–23, https://doi.org/10.1080/10448639908233684, 1999. a
Li, D., Schrön, M., Köhli, M., Bogena, H., Weimar, J., Jiménez Bello, M., Han, X., Martínez Gimeno, M., Zacharias, S., Vereecken, H., and Hendricks Franssen, H.: Can Drip Irrigation be Scheduled with CosmicRay Neutron Sensing?, Vadose Zone J., 18, 190053, https://doi.org/10.2136/vzj2019.05.0053, 2019. a, b
Liu, H., Hou, Y., Li, H., Song, Y., Hu, L., and Liang, M.: Cosmicray neutron fluxes and spectra at different altitudes based on Monte Carlo simulations, Appl. Radiat. Isotopes, 175, 109800, https://doi.org/10.1016/j.apradiso.2021.109800, 2021. a
Maire, E. and Withers, P.: Quantitative Xray tomography, Int. Mater. Rev., 59, 1–43, https://doi.org/10.1179/1743280413Y.0000000023, 2013. a
Mares, V. and Schraube, H.: Evaluation of the response matrix of a Bonner Sphere Spectrometer with LiI detector from thermal energy to 100 MeV, Nucl. Instrum. Meth. A, 337, 461–473, https://doi.org/10.1016/01689002(94)911169, 1994. a
Mares, V., Schraube, G., and Schraube, H.: Calculated neutron response of a Bonner Sphere Spectrometer with ^{3}He counter, Nucl. Instrum. Meth. A, 307, 398–412, https://doi.org/10.1016/01689002(91)90210H, 1991. a, b
Matsumoto, M. and Nishimura, T.: Mersenne Twister: A 623dimensionally Equidistributed Uniform Pseudorandom Number Generator, ACM T. Model. Comput. Sim., 8, 3–30, https://doi.org/10.1145/272991.272995, 1998. a
McConn Jr., R., Gesh, C., Pagh, R., Rucker, R., and Williams III, R.: Compendium of Material Composition Data for Radiation Transport Modeling, Tech. Rep. PNNL15870 Rev. 1, Pacific Northwest National Laboratory, Richland, Washington 99352, https://www.pnnl.gov/main/publications/external/technical_reports/pnnl15870rev1.pdf, 2011. a
McKinney, G.: MCNP6 Cosmic and Terrestrial Background Particle Fluxes, LAUR1324293, release 3, https://mcnp.lanl.gov/pdf_files/laur1324293.pdf, 2013. a
McKinney, G., Armstrong, H., James, M., Clem, J., and Goldhagen, P.: MCNP6 CosmicSource Option, LAUR1222318, https://permalink.lanl.gov/object/tr?what=info:lanlrepo/lareport/LAUR1222318, 2012. a
Metropolis, N. and Ulam, S.: The Monte Carlo Method, J. Am. Stat. Assoc., 44, 335–341, https://doi.org/10.2307/2280232, 1949. a
Nelms, A.: Energy loss and range of electrons and positrons, United States Department of Commerce, national Bureau of Standards Circular 577 and Suppl., https://catalog.hathitrust.org/Record/007291096 (last access: 18 January 2023), 1956. a
Niita, K.: QMD and JAM Calculations for High Energy NucleonNucleus Collisions, J. Nucl. Sci. Technol., 39, 714–719, https://doi.org/10.1080/00223131.2002.10875198, 2002. a
Niita, K., Takada, H., Meigo, S., and Ikeda, Y.: Highenergy particle transport code NMTC/JAM, Nucl. Instrum. Meth. B, 184, 406–420, https://doi.org/10.1016/S0168583X(01)007844, 2001. a
Otuka, N., Dupont, E., Semkova, V., Pritychenko, B., Blokhin, A., Aikawa, M., Babykina, S., Bossant, M., Chen, G., Dunaeva, S., Forrest, R., Fukahori, T., Furutachi, N., Ganesan, S., Ge, Z., Gritzay, O., Herman, M., Hlavač, S., Kato, K., Lalremruata, B., Lee, Y., Makinaga, A., Matsumoto, K., Mikhaylyukova, M., Pikulina, G., Pronyaev, V., Saxena, A., Schwerer, O., Simakov, S., Soppera, N., Suzuki, R., Takács, S., Tao, X., Taova, S., Tárkányi, F., Varlamov, V., Wang, J., Yang, S., Zerkin, V., and Zhuang, Y.: Towards a More Complete and Accurate Experimental Nuclear Reaction Data Library (EXFOR): International Collaboration Between Nuclear Reaction Data Centres (NRDC), Nucl. Data Sheets, 120, 272–276, https://doi.org/10.1016/j.nds.2014.07.065, 2014. a
Pazianotto, M., CortésGiraldo, M., Federico, C., Gonçalez, O., Quesada, J., and Carlson, B.: Determination of the cosmicrayinduced neutron flux and ambient dose equivalent at flight altitude, J. Phys. Conf. Ser., 630, 012022, https://doi.org/10.1088/17426596/630/1/012022, 2015. a
Peggs, S., et al.: European Spallation Source Technical Design Report, Technical Design Report ESS, ESS, Lundt, ESS2013001, http://docdb01.esss.lu.se/DocDB/0002/000274/006/TDR_final_130423_print_ch1.pdf, 2013. a
Prael, R. and Lichtenstein, H.: User Guide to LCS: The LAHET Code System, LAUR893014, https://permalink.lanl.gov/object/tr?what=info:lanlrepo/lareport/LAUR893014, 1989. a
Rasche, D., Köhli, M., Schrön, M., Blume, T., and Güntner, A.: Towards disentangling heterogeneous soil moisture patterns in cosmicray neutron sensor footprints, Hydrol. Earth Syst. Sci., 25, 6547–6566, https://doi.org/10.5194/hess2565472021, 2021. a
Romano, P. and Forget, B.: The OpenMC Monte Carlo particle transport code, Ann. Nucl. Energ., 51, 274–281, https://doi.org/10.1016/j.anucene.2012.06.040, 2013. a, b
Roth, S.: Ray casting for modeling solids, Comput. Vision Graph., 18, 109–144, https://doi.org/10.1016/0146664X(82)901691, 1982. a
Šaroun, J. and Kulda, J.: RESTRAX – a program for TAS resolution calculation and scan profile simulation, Physica B, 234236, 1102–1104, https://doi.org/10.1016/s09214526(97)000379, 1997. a
Sato, T.: Analytical Model for Estimating Terrestrial Cosmic Ray Fluxes Nearly Anytime and Anywhere in the World: Extension of PARMA/EXPACS, PLOS ONE, 10, 1–33, https://doi.org/10.1371/journal.pone.0144679, 2015. a
Sato, T.: Analytical Model for Estimating the Zenith Angle Dependence of Terrestrial Cosmic Ray Fluxes, PLoS ONE, 11, 1–22, https://doi.org/10.1371/journal.pone.0160390, 2016. a, b, c, d, e, f
Sato, T. and Niita, K.: Analytical Functions to Predict CosmicRay Neutron Spectra in the Atmosphere, Radiation Research, 166, 544–555, https://doi.org/10.1667/RR0610.1, 2006. a
Sato, T., Yasuda, H., Niita, K., Endo, A., and Sihver, L.: Development of PARMA: PHITSbased Analytical Radiation Model in the Atmosphere, Radiat. Res., 170, 244–259, https://doi.org/10.1667/RR1094.1, 2008. a, b, c
Schattan, P., Köhli, M., Schrön, M., Baroni, G., and Oswald, S.: Sensing AreaAverage Snow Water Equivalent with CosmicRay Neutrons: The Influence of Fractional Snow Cover, Water Resour. Res., 55, 10796–10812, https://doi.org/10.1029/2019WR025647, 2019. a, b
Schrön, M.: Cosmicray neutron sensing and its applications to soil and land surface hydrology, PhD thesis, University of Potsdam, https://nbnresolving.org/urn:nbn:de:kobv:517opus4395433, 2017. a, b
Schrön, M., Köhli, M., Scheiffele, L., Iwema, J., Bogena, H. R., Lv, L., Martini, E., Baroni, G., Rosolem, R., Weimar, J., Mai, J., Cuntz, M., Rebmann, C., Oswald, S. E., Dietrich, P., Schmidt, U., and Zacharias, S.: Improving calibration and validation of cosmicray neutron sensors in the light of spatial sensitivity, Hydrol. Earth Syst. Sci., 21, 5009–5030, https://doi.org/10.5194/hess2150092017, 2017. a, b
Schrön, M., Rosolem, R., Köhli, M., Piussi, L., Schröter, I., Iwema, J., Kögler, S., Oswald, S., Wollschläger, U., Samaniego, L., Dietrich, P., and Zacharias, S.: Cosmicray Neutron Rover Surveys of Field Soil Moisture and the Influence of Roads, Water Resour. Res., 54, 6441–6459, https://doi.org/10.1029/2017WR021719, 2018. a
Schrön, M., Zacharias, S., Womack, G., Köhli, M., Desilets, D., Oswald, S. E., Bumberger, J., Mollenhauer, H., Kögler, S., Remmler, P., Kasner, M., Denk, A., and Dietrich, P.: Intercomparison of cosmicray neutron sensors and water balance monitoring in an urban environment, Geosci. Instrum. Method. Data Syst., 7, 83–99, https://doi.org/10.5194/gi7832018, 2018. a, b, c
Schrön, M., Oswald, S., Zacharias, S., Kasner, M., Dietrich, P., and Attinger, S.: Neutrons on Rails: Transregional Monitoring of Soil Moisture and Snow Water Equivalent, Geophys. Res. Lett., 48, e2021GL093924, https://doi.org/10.1029/2021gl093924, 2021. a
Shibata, K., Iwamoto, O., Nakagawa, T., Iwamoto, N., Ichihara, A., Kunieda, S., Chiba, S., Furutaka, K., Otuka, N., Ohsawa, T., Murata, T., Matsunobu, H., Zukeran, A., Kamada, S., and Katakura, J.: JENDL4.0: A New Library for Nuclear Science and Engineering, J. Nucl. Sci. Technol., 48, 1–30, https://doi.org/10.1080/18811248.2011.9711675, 2011. a
Smith, A., Fields, P., and Roberts, J.: Spontaneous Fission Neutron Spectrum of ^{252}Cf, Phys. Rev., 108, 411–413, https://doi.org/10.1103/PhysRev.108.411, 1957. a
Solovyev, A., Fedorov, V., Kharlov, V., and Stepanova, U.: Comparative analysis of MCNPX and GEANT4 codes for fastneutron radiation treatment planning, Nucl. Energ. Technol., 1, 14–19, https://doi.org/10.1016/j.nucet.2015.11.004, 2015. a
Sun Programmers Group: Fortran 77 Reference Manual, Tech. rep., SunSoft, http://wwwcdf.pd.infn.it/localdoc/f77_sun.pdf, 1995a. a
Sun Programmers Group: Fortran 90 User's Guide, Tech. rep., SunSoft, http://smdc.sinp.msu.ru/doc/Fortran90UsersGuide.pdf, 1995b. a
Sweezy, J., Hertel, N., and Veinot, K.: BUMS – Bonner Sphere Unfolding Made Simple: an HTML based multisphere neutron spectrometer unfolding package, Nucl. Instrum. Meth. A, 476, 263–269, https://doi.org/10.1016/S01689002(01)014668, 2002. a
Terrell, J.: Fission Neutron Spectra and Nuclear Temperatures, Phys. Rev., 113, 527–541, https://doi.org/10.1103/PhysRev.113.527, 1959. a
Thomas, D. and Alevra, A.: Bonner Sphere Spectrometers – a critical review, Nucl. Instrum. Meth. A, 476, 12–20, https://doi.org/10.1016/S01689002(01)013791, 2002. a
Trkov, A., Herman, M., and Brown, D.: ENDF6 Formats Manual, National Nuclear Data Center, Brookhaven National Laboratory, data Formats and Procedures for the Evaluated Nuclear Data Files ENDF/BVI and ENDF/BVII, https://wwwnds.iaea.org/public/endf/endfmanual.pdf, 2012. a, b
van der Ende, B., Atanackovic, J., Erlandson, A., and Bentoumi, G.: Use of GEANT4 vs. MCNPX for the characterization of a boronlined neutron detector, Nucl. Instrum. Meth. A, 820, 40–47, https://doi.org/10.1016/j.nima.2016.02.082, 2016. a
Walsh, J., Forget, B., and Smith, K.: Accelerated sampling of the free gas resonance elastic scattering kernel, Ann. Nucl. Energ., 69, 116–124, https://doi.org/10.1016/j.anucene.2014.01.017, 2014. a
Waters, L., McKinney, G., Durkee, J., Fensin, M., Hendricks, J., James, M., Johns, R., and Pelowitz, D.: The MCNPX Monte Carlo Radiation Transport Code, AIP Conf. Proc., 896, 81–90, https://doi.org/10.1063/1.2720459, 2007. a, b
Watt, B.: Energy Spectrum of Neutrons from Thermal Fission of ^{235}U, Phys. Rev., 87, 1037–1041, https://doi.org/10.1103/PhysRev.87.1037, 1952. a
Wechsler, D., Zsigmond, G., Streffer, F., and Mezei, F.: VITESS: Virtual instrumentation tool for pulsed and continuous sources, Neutron News, 11, 25–28, https://doi.org/10.1080/10448630008233764, 2000. a
Weimar, J., Köhli, M., Budach, C., and Schmidt, U.: LargeScale BoronLined Neutron Detection Systems as a ^{3}He Alternative for Cosmic Ray Neutron Sensing, Front. Water, 2, 16, https://doi.org/10.3389/frwa.2020.00016, 2020. a, b
Weisskopf, V.: Statistics and Nuclear Reactions, Phys. Rev., 52, 295–303, https://doi.org/10.1103/PhysRev.52.295, 1937. a
X5 Monte Carlo Team: MCNPA general Monte Carlo Nparticle transport code, Version 5, LAUR031987, volume I: Overview and Theory, https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/laur031987.pdf, 2003. a
Yamashita, M., Stephens, L., and Patterson, H.: Cosmicrayproduced neutrons at ground level: Neutron production rate and flux distribution, J. Geophys. Res., 71, 3817–3834, https://doi.org/10.1029/JZ071i016p03817, 1966. a
Ziegler, J.: Terrestrial cosmic ray intensities, IBM J. Res. Dev., 42, 117–140, https://doi.org/10.1147/rd.421.0117, 1998. a
Zreda, M., Desilets, D., Ferré, T., and Scott, R.: Measuring soil moisture content noninvasively at intermediate spatial scale using cosmicray neutrons, Geophys. Res. Lett., 35, L21402, https://doi.org/10.1029/2008GL035655, 2008. a
Zreda, M., Shuttleworth, W. J., Zeng, X., Zweck, C., Desilets, D., Franz, T., and Rosolem, R.: COSMOS: the COsmicray Soil Moisture Observing System, Hydrol. Earth Syst. Sci., 16, 4079–4099, https://doi.org/10.5194/hess1640792012, 2012. a
Zweck, C., Zreda, M., and Desilets, D.: Snow shielding factors for cosmogenic nuclide dating inferred from Monte Carlo neutron transport simulations, Earth Planet. Sc. Lett., 379, 64–71, https://doi.org/10.1016/j.epsl.2013.07.023, 2013. a
For example natural boron; see also Table 2, contains ≈80 % ^{11}B and ≈20 % ^{10}B. For a reaction the cross section ${\mathrm{\Sigma}}_{i}={\mathrm{\Sigma}}_{i}\left({}^{\mathrm{10}}\mathrm{B}\right)+{\mathrm{\Sigma}}_{i}\left({}^{\mathrm{11}}\mathrm{B}\right)$ is accompanied by a vector
of individual contributions $\left[{\mathrm{\Sigma}}_{i}\right(\mathrm{0}),{\mathrm{\Sigma}}_{i}(\mathrm{0})+{\mathrm{\Sigma}}_{i}(\mathrm{1}\left)\right]$ and a vector
of isotopes [^{10}B,^{11}B].
If i=0, then ξ≤v_{s}[1].
 Abstract
 Introduction
 Calculation routines
 Model design
 Computational structure
 Sources and energy
 Calculation scheme
 Detector configurations
 Evaluation of the URANOS model
 Performance benchmarks
 User interface
 Conclusions
 Appendix A: Elements, isotopes and reaction types
 Appendix B: Material codepages
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Calculation routines
 Model design
 Computational structure
 Sources and energy
 Calculation scheme
 Detector configurations
 Evaluation of the URANOS model
 Performance benchmarks
 User interface
 Conclusions
 Appendix A: Elements, isotopes and reaction types
 Appendix B: Material codepages
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References