Issue 
EPJ Nuclear Sci. Technol.
Volume 7, 2021



Article Number  12  
Number of page(s)  15  
DOI  https://doi.org/10.1051/epjn/2021010  
Published online  15 June 2021 
https://doi.org/10.1051/epjn/2021010
Regular Article
Muon radiography to visualise individual fuel rods in sealed casks
^{1}
GRS gGmbH,
Schwertnergasse 1,
50667
Cologne,
Germany
^{2}
GRS gGmbH,
Boltzmannstr. 14,
85748
Garching,
Germany
^{3}
GRS gGmbH,
Kurfürstendamm 200,
10719
Berlin,
Germany
^{4}
BGZ mbH,
Dammstraße,
84051
Essenbach,
Germany
^{*} email: thomas.braunroth@grs.de
Received:
26
February
2021
Received in final form:
4
May
2021
Accepted:
11
May
2021
Published online: 15 June 2021
Cosmicray muons can be used for the nondestructive imaging of spent nuclear fuel in sealed dry storage casks. The scattering data of the muons after traversing provides information on the thereby penetrated materials. Based on these properties, we investigate and discuss the theoretical feasibility of detecting single missing fuel rods in a sealed cask for the first time. We perform simulations of a vertically standing generic cask model loaded with fuel assemblies from a pressurized water reactor and muon detectors placed above and below the cask. By analysing the scattering angles and applying a significance ratio based on the KolmogorovSmirnov test statistic we conclude that missing rods can be reliably identified in a reasonable measuring time period depending on their position in the assembly and cask, and on the angular acceptance criterion of the primary, incoming muons.
© T. Braunroth et al., Published by EDP Sciences, 2021
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
The operation of nuclear power plants generates highlevel radioactive wastes which need to be stored and disposed of. A wellestablished part of nuclear waste management is the dry storage of used fuel assemblies in designated casks. Depending on the availability of a final repository, the fuel assemblies might remain inside the casks placed in interim storage facilities for decades. The casks are designed to enclose the highlevel radioactive waste and separate it from the biosphere.
Inspecting the interior of a storage cask directly would require an opening of the cask, a difficult task due to the radiation. So far, conventional radiography using neutrons or photons could not be applied successfully, partly due to the rich scattering history of traversing particles as a direct result of the dimensions of the storage cask [1]. Other methods such as threedimensional temperature field measurements or antineutrino monitoring [2] require a detailed knowledge of the fuel history or are not suitable for the assessment of individual storage casks. Cosmic muons, created directly or indirectly in the atmosphere by the interaction of cosmic radiation with particles, have been also used for imaging purposes (muography) [3]. These muons and antimuons (muons in the following) are characterized by a broad energy spectrum spanning several orders of magnitude and have a mean momentum of approximately 4 GeV/c [4]. For muons with an absolute momentum above 1 GeV/c the integral vertical intensity I_{V} is approximately 70 m^{−2} s^{−1} sr^{−1} [5,6] and the decrease of the flux intensity scales approximately with cos(θ), with θ being the angle with respect to the vertical [4]. In experimental physics, a value of I_{V} ≈ 1 cm^{−2} min^{−1} has been established for working with horizontal detectors.
Alvarez and colleagues were one of the first to use muons for noninvasive imaging and published their landmark paper on the search for hidden chambers in the pyramids of Giza in 1970 [7]. These first studies were based on the measurement of the attenuation of the cosmic muon flux and provided twodimensional projection images (muon radiography). Since then, muon radiography has been used in various fields such as the study of volcanoes [8], geological applications [9], the identification of cavities in archaeology [10] as well as in industrial applications [11].
In 2003, a Los Alamos research group proposed using the scattering angle of the outgoing muons as the basic information for imaging [12]. This approach requires the measurement of the incoming as well as the outgoing trajectories of the muons and allows to obtain threedimensional images (muon tomography) of volumes not exceeding tens of meters. This technique has already been used in various fields, such as nuclear control [13], transport control [14] and the monitoring of historical buildings [15]. In addition to experimental studies, MonteCarlo simulations play an important role in muon imaging, e.g. in terms of feasibility studies or with respect to the detector design.
The application of muography for the purpose of noninvasive control and monitoring of the interior of dry storage casks has gained an increasing interest and fostered experimental and simulation studies. Besides the fundamental suitability of muography for this purpose, these efforts also addressed methodological and time requirements. Durham et al. [16] applied muonscattering radiography to a Westinghouse MC10 cask and showed experimentally that cosmic muons can indeed be used to determine if spent fuel assemblies are missing without the need to open the cask. A number of simulation studies were performed using two planar tracking detectors placed on opposite sides of an object, each focusing on different aspects, for example: Jonkmans et al. [17] investigated the capabilities of muons to image the contents of shielded containers to detect enclosed nuclear materials with highZ. Clarkson et al. [18] performed Geant4 simulations of a scintillatingfibre tracker for tomographic scans of legacy nuclear waste containers. Chatzidakis [19] applied a Bayesian approach to monitor sealed dry casks to infer on the amount of spent nuclear fuel and to investigate the limitiations of this approach. Using the attenuation and scattering characteristics of the muons derived from Geant4 simulations, Ambrosino et al. [20] found that a 10 cm^{3} Uranium block inside a concrete structure could be identified after a onemonth period of measurement. Poulsen et al. [21] were the first to apply filtered backprojection algorithms to muon tomography imaging of dry storage casks using simulated data and could show that this technique can be applied to the detection of missing fuel assemblies. In a more recent work,Poulsen et al. [22] used the experimental data of reference [16] for a numerical study using Geant4 to distinguish different loads of a cask. The study indicates that a oneweek muon measurement for the given experimental setup is sufficient to detect a missing fuel assembly or to identify a dummy assembly made out of iron or lead.
With respect to the resolving power, both experimental and simulation studies have been focused on the level of fuel assemblies so far. In addition, the majority of studies are based on a transversal configuration, where the detectors are placed on the sides of the cask.
In this study, we use MonteCarlo simulations to investigate a longitudinal configuration, with the detectors placed above and below the cask. We will investigate if muography allows detecting individual missing fuel rods. To unravel insights independent of reconstruction algorithms, this work will focus on radiographic images based on transmission ratios as well as scatteringangle information.
The guiding questions of this work are as follows: is it possible to even detect individual missing fuel rods with muon radiography? If so, are there any constraints or requirements with respect to the experimental setup and how much time does a measurement require? What can be used as a significance measure to detect a missing fuel rod? Does the significance depend on the relative position of the considered fuel rod within the fuel assembly and is the significance dependent on the number of considered events?
This contribution is structured as follows: in Section 2 we provide a detailed description of the simulation tool as well as the investigated geometry. Moreover, we describe the data processing and aspects related to the validation of the simulation. Section 3 features the analysis and discussion. We address two levels of detail concentrating on the recognition of (missing) fuel assemblies and individual fuel rods. This is followed by a summary and conclusion in Section 4. We end with a short outlook in Section 5.
2 Simulation and data processing
Simulations were performed with a dedicated tool based on the MonteCarlo toolkit Geant4 [23–25]. Geant4 allows simulating the passage of particles through matter and has been used for numerous applications, including high energy physics, nuclear physics, accelerator physics and others.
In this section we describe the key aspects of the tool, i.e. the geometry (Sect. 2.1), the treatment of primary particle properties (Sect. 2.2), aspects related to physics and tracking (Sect. 2.3) as well as optimization strategies (Sect. 2.4) to reduce the computational time.
The tool was compiled against v10.06.p2 of Geant4 and allows using multithreading. The results are written eventbyevent into ROOT [26,27] container files, which allows performing postprocessing in a flexible manner. Finally, Section 2.5 discusses the validation of the code by comparing simulated and tabulated (or empirically established) energy losses and angular straggling for different target materials and projectile energies.
2.1 Geometry
2.1.1 Generic cask
The key component of the geometry is a generic cask model (referred to as generic model in the following) which mimics the features of the CASTOR^{®} V/19 cask [28], e.g. in terms of major components, dimensions, materials and masses. The CASTOR^{®} V/19 cask is used for transport as well as storage purposes and is designed to carry up to 19 fuel assemblies from pressurized water reactors (PWR). All information on geometries and materials specifying the generic model was derived from public available sources such as references [28,29].
A visualisation of the generic model was generated with the Geant4 OpenGL interface and is depicted in Figure 1. The individual components of the generic model can be identified in the exploded view shown in Figure 2.
A comparison of some key properties of the generic model on the one hand and the CASTOR^{®} V/19 cask on the other hand can be found in Table 1 and highlights the mutual similarities.
Each of the 19 fuel compartments can be occupied with one (modelled) fuel assembly. Each modelled fuel assembly comprises top and bottomnozzle, fuel rods as well as control rods. The fuel rods consist of nuclear fuel (UO_{2}) and cladding tubes (Zirconium alloy). A complete 18x1824 fuel assembly consists of 300 fuel rods and 24 control rods. The arrangement of fuel and control rods within a complete fuel assembly is shown in Figure 3. The modelled top and base components of the assembly are simplified. They are assumed to be boxlike, with heights chosen to comply with real masses. Basic properties related to the modelled fuel assemblies are summarized in Table 2.
All parts of the Geant4 model were derived from basic Geant4 solid objects and refined with Boolean operations. Each component can be switched on and off by a commandbased user interface, which easily allows performing simulations for different target geometries. In addition, it is possible to remove arbitrary components from a fuel assembly, e.g. specific fuel rods at specific slot positions. This allows, among others, for investigating the contribution of specific fuel rods to radiographic (or tomographic) images in more detail.
Fig. 1
Visualisation of the generic model and its upward orientation within the present work. The coordinate system as it is used in this study is provided in the lower right corner. The origin of the coordinate system coincides with the center position of the bottom part of the model. The red and green areas above and below the model indicate the incoming and outgoing detectors. 
Fig. 2
Exploded view of the generic model showing its individual components: monolithic body (1), basket with 19 fuel compartments (2), trunnions (3), primary lid (4), polyethylene plate (5), secondary lid (6), protection plate (7), inner moderator rods and plugs (8), outer moderator rods and plugs (9), polyethylene plate (10), base plate (11) and a representative fuel assembly (12). See text for details. 
Comparison of key dimensions between the generic model and the CASTOR^{®} V/19 cask (in storage configuration).
Fig. 3
Left: topview on the arrangement of fuel rods (yellow), cladding tubes (green) as well as control rods (blue) in case of a complete fuel assembly. Each element can be specified unambiguously based on its slot position given by x_{id}, y_{id}. Right: same as left figure, but several elements along the diagonal are missing (x_{id} ∕y_{id}= 1, 5, 9). This configuration is investigated in more detail in Section 3.2. 
Properties of the modelled fuel assemblies.
2.1.2 Boxlike object
Instead of the generic model, it is possible to generate a boxlike object, whose basic properties – i.e. dimensions, placement and material – can be specified individually for each simulation run. All interactions within this box are recorded within ROOT container files, which is particularly useful for validation purposes.
2.1.3 Coordinate system
The zaxis coincides with the symmetry axis of the generic model and is oriented upwards, i.e. its orientation is selected so that the zcoordinate of the model’s top is larger than the zcoordinate of its bottom, z_{top} > z_{bottom}. The x and y axis are orientated in such a way that the (x, y, z)coordinate system generates a righthanded euclidean space. The orientation of the individual axes is indicated in Figure 1. The angle θ reflects the angle of a given vector r⃗ and the inverted zaxis. Using this convention, a muon from the zenith is characterized by θ = 0°. The angle φ is the angle between the projection of a vector r⃗ onto the (x, y)plane and ê_{x}.
2.1.4 Detectors
Detector systems are mimicked by two rectangular detector planes (vanishing thickness, area of (3 × 3) m^{2}), whose normal vectors are parallel to ê_{z}. The detector plane placed above the generic model is called incoming detector whereas the detector plane placed below the generic model is called outgoing detector. The gap between the detector surface and the top (or the bottom) of the generic cask is ≈ 10 cm. Muons traversing these planes are tracked and key information is determined, see Section 2.3.
This approach gives access to at least as many properties as a real detection system for muon tomography may provide. Since the simulation allows determining all properties precisely on an eventbyevent base, it benefits from an infinite resolving power.^{1} In this regard, it provides a bestcase scenario and can be considered as a first important theoretical step towards substantiated feasibility studies.
2.2 Primary particles
Primaries are generated by a userdefined class based on the G4VUserPrimaryGeneratorAction class provided by Geant4.
The MonteCarlo approach is realized by using uniform random number generators with limits a, b () to determine for each event k the particle properties, i.e. particle type (muon (μ^{−}) or antimuon (μ^{+})), initial momentum information (, θ_{k}, φ_{k}) as well as initial position information (x_{k}, y_{k}). Details are described in the following.
2.2.1 Particle type
All primaries are muons or antimuons, assuming a charge ratio μ^{+} ∕μ^{−} equal to 1.28 [30]. The particle typefor the event k is determined assuming  if the generated random number is smaller than (or equal to) 1.28, the generated primary for this event will be an antimuon (μ^{+}), otherwise it will be a muon (μ^{−}).
2.2.2 Initial particle momentum
The description of the initial muon momentum is based on an empirical and parametric approach [31], according to which the muon intensity at any given momentum and angle to the vertical intensity I_{V} is given by (1)
Here, I_{V} is given by the Bugaev parametrisation [32]: (2)
The parameters c_{1} to c_{5} were determined in reference [31] by a fitting approach and are quantified as:
Due to the lack of a suitable random number generator to directly address the associated probability distribution in Geant4, the analytical description is discretized in a twodimensional pattern as follows.
Firstly, the polar component of the angular spectrum represented by θ_{in}, ranging from θ_{min}= 0° to an upper limit of θ_{max} = 25°, is split into n bins b_{θ,i} (i = 1, …, n). The quoted upper limit of 25° limits the contribution of trajectories that do not cross the full longitudinal length of the cask. All b_{θ,i} span over identical angular ranges. For each angular bin b_{θ,i}, the lower and upper limits are then given by and .
The function , for i = 1, …, n, can then be associated to each of these bins of the angular spectrum, where is the arithmetic mean angle of the specific bin. These assigned functions I_{i} are then integrated numerically within the kinetic energy range of T ∈ [1 GeV;1 TeV]. The calculated integrals quantify the relative weights w_{i} for each of the angular bins b_{θ,i}: (3)
Information on , and w_{i} are stored in a dedicated text file that is used as an input parameter to the simulation. uses these weights w_{i} to determine for each event k the proper angular bin b_{θ,k}. In a next step, the proper polar angle θ_{k} for the present event k is determined using . The azimuthal angle φ_{k} is determined for each event k using .
Secondly, the absolute momentum spectra are treated comparably. The momentum axis from p_{μ} (T = 1 GeV) to p_{μ} (T = 1 TeV) is discretized into m bins (j = 1, …, m) with increasing bin sizes. Each bin is characterized by the following integral: (4)
Here, and denote the lower and upper limits of the bin .
For each angular bin b_{θ,i}, values for , and v_{i,j} are stored in dedicated text files which are used as mandatory input information for the simulation code.
uses these weights v_{i,j} to determine for each event k the proper momentum bin . Finally, the proper p_{μ,k} is determined using . Within the present work, all numerical integrations were performed with GNU OCTAVE [33].
Figure 4 shows histograms of the simulated angular and momentum distributions for θ_{min} = 0°, θ_{max} = 25°, n = 10 and m = 172.
In addition to this distributionbased approach described above, the tool also allows performing simulations with monoenergetic and monodirectional muons.
Fig. 4
Sample distributions of the initial angle θ_{in} (left) and of the absolute muon momentum p (right) as determined with the incoming detector. The red curve shown in the left spectrum is proportional to , while the red curve shown in the right spectrum is proportional to equation (1) with θ = 12.5°. 
2.2.3 Initial particle position
The initial positions of the muons are selected as follows. The zcoordinate is a fixed value and selected so that it is ensured that the muon is created just above the incoming detector. The x and ycoordinates are selected from uniform distributions covering ranges from x_{min} to x_{max} () and y_{min} to y_{max} (), respectively.The associated limits can be specified in the input file of the simulation.
2.3 Physics and tracking
This code uses the modular physics list FTFP_BERT implemented in Geant4, which is recommended by the Geant4 developers for highenergy physics.
The incoming and outgoing detector volumes are linked to a dedicated tracking class. This allows tracking and storing information on the particle properties needed for the post processing to generate radiographic (or tomographic) images.
2.4 Cutting conditions to reduce computational time
Muons may interact with the materials of the generic model according to the chosen physics list in various ways. Some of these interactions are irrelevant to the present work. For example, some interactions may generate secondary particles – e.g. photonsor neutrons – whose tracking consumes computational time without any benefits or consequences to the generated images. To avoid an unnecessary increase of the computational time, the trajectories of such secondary particles are annihilated at the end of their first steps. It must be stressed that this does not concern the possible interactions of muons with matter  all interactions which may lead to such secondary particles still take place.
In addition and as already described in Section 2.2, the muon properties – with respect to the angle θ and kinetic energy T – were cut to limit the simulation of muon trajectories that are not useful for an analysis based on a twodetector setup. Finally, all muon trajectories are annihilated at the exit of the outgoing detector.
2.5 Validation
This section discusses the validation of two key aspects with respect to muon imaging. The first aspect deals with the energy loss of muons in matter and is discussed in Section 2.5.1. The second aspect deals with the angular scattering of muons after traversing matter with a known thickness, see Section 2.5.2. In both cases, the referenced thickness was specified to 1 mm.
2.5.1 Energy loss of Muons in various target materials
In this section, we compare simulated energy losses of muons in relevant target materials – uranium dioxide, polyethylene, stainless steel, ductile iron and zirconium alloy – to tabulated (or calculated) values of the mean differential energy loss [34,35].
The quoted references provide only for two of the listed compound materials – uranium dioxide and polyethylene – tabulated values. Hence, reference values were calculated for the other compounds – stainless steel, ductile iron and zirconium alloy – according to Bragg’s Rule of Stopping Power Additivity [36]:
Here, w_{j} denotes the mass fraction of the material j, dE/dx_{j} is the differential energy loss in the material j and dE/dx describes the mean differential energy loss by the muons in the compound material.
Strictly speaking, the simulated energy losses per distance ΔE are not identical to differential energy losses. To increase the comparability, we consider for both numbers a distance of 1 mm and cut the lower limit of the considered kinetic energy range to ensure that the (simulated) mean energy losses are smaller (or comparable) to five percent of the kinetic energy T, i.e. ΔE ≲ 0.05 T. The upper limit of the kinetic energy range is set to T = 1 TeV which equals the maximal kinetic energy of the considered muon primaries in the present study.
The number of simulated events was chosen in a way that the relative uncertainties of the mean values ΔE extracted from the simulated energyspectra are below 2%. The only exception is given for polyethylene, for which the limit is increased to 5%. In total, between 10^{5} and 5 ⋅ 10^{6} events were simulated for each kinetic energy and each material. The results are shown in Figure 5.
For the materials stainless steel, ductile iron and zirconium alloy, the agreement between simulated and tabulated values is very good and the majority of values deviate by less than 2%. Few exceptions were found for kinetic energies close to 1 GeV. It shall be mentioned that due to the high ΔE dispersion at high projectile energies, the mean values are quite sensitive to statistical outliers.
The largest deviations were found for polyethylene, for which relative deviations close to 7% were determined over a broad kinetic energy range.
For uranium dioxide, the agreement for kinetic energies between T = 80 MeV and T = 10 GeV is very good. However, with increasing kinetic energy a clear trend towards larger deviations of up to ~ 6% was observed.
In summary, the overall agreement between simulated and tabulated (or calculated) stopping power values is convincing. In general and with respect to the given kinetic energy limits, the models used in the simulation have a tendency to yield slighter larger values of the stopping power in the investigated materials compared to the reference tabulations of references [34,35].
Fig. 5
Comparison of simulated mean energy losses ΔE per distance and mean differential energy losses dE/dx (both in [MeV/mm]) of muons in various target materials. The top figure shows simulated mean energy losses (symbols) and compares them to tabulated reference values (solid lines). The red curve (boxes) shows the results for uranium dioxide, theorange curve (circles) corresponds to polyethylene, the blue curve (pentagons) corresponds to stainless steel, the green curve (uppointing triangles) corresponds to ductile iron and the violet curve (downpointing triangles) corresponds to zirconium alloy. The lower five figures show for each material the relative deviations between simulated energy loss per mm and tabulated mean stopping power values. 
2.5.2 Angular scattering of Muons in various target materials
In this section, we investigate simulated scattering angles θ^{S} of muons after traversing different target materials and compare the width of the associated distributions to an established semiempirical description. Each material had a thickness of 1 mm.
According to reference [37], the rootmean square of the scattering angle in the x and yplanes can be described using the following formula: (7)
Here, β is given by v∕c, p is the absolute muon momentum, Δz is the thickness of the irradiated material and X_{0} is the radiation length of the material. The rootmean square of the scattering angle θ^{S} in space is given by: (8)
In case of a compound material, the radiation length X_{0} can be calculated using the formula (9)
where w_{j} is the mass fraction of the component j and X_{0,j} is the radiation length of the component j. The radiation lengths were either taken directly from references [34,35] or calculated according to equation (9). An overview of the various radiation lengths X_{0} for the relevant materials is given in Table 3.
Simulations covering a broad range of kinetic energies were performed for all relevant materials with 10^{5} events. With respect to the energy range, the same restrictions as for the discussion of the energy loss were applied. The simulated results are shown and compared to calculated values in Figure 6.
The deviations between simulated and semiempirical values for are usually less than 3%. Only for UO_{2} larger deviations of about 6% occur. In all cases, the deviations as a function of the kinetic energy remain rather constant. Overall, the agreement between simulated and calculated values is very good.
Radiation lengths X_{0} for materials relevant to the present work.
Fig. 6
Comparison of rootmean square (RMS) values of scattering angles (in [mrad]) between semiempirical calculations and simulated values for muons in various target materials with a thickness of 1 mm. The top figure shows values derived from the simulation (symbols) and compares them to calculated values based on a semiempirical approach (solid lines). The red curve (boxes) shows the results for uranium dioxide, the orange curve (circles) corresponds to polyethylene, the blue curve (pentagons) corresponds to stainless steel, the green curve (uppointing triangles) corresponds to ductile iron and the violet curve (downpointing triangles) corresponds to zirconium alloy. The lower five figures show for each material therelative deviations between simulated and calculated values. 
3 Analysis and discussion
The analysis is split into two main sections. The first section investigates longitudinal scans of the generic model and analyses simulated radiographic images. Here, the focus lies on the effects of various angular acceptances on the image quality. In addition, we compare transmission radiographic images with scattering radiographic images, e.g. in terms of image contrast and resolving power. The presentation focuses on qualitative aspects and is described in more detail in Section 3.1.
The second section provides a detailed study that focuses on the central fuel assembly. In particular, it investigates the suitability of muonscattering radiography to make statements about the occupancy of individual fuelrod slots within a certain fuel assembly. This section focuses on quantitative aspects and is described in Section 3.2 in more detail.
3.1 Longitudinal scans covering the full cross section of the generic model
All simulations discussed in this section were performed with 5 × 10^{7} events. All primary muons were generated at a fixed height (z ≈ 6.1 m) just above the generic model, while the initial x and ycoordinates followed uniform distributions within the limits of − 1.25 m ≤ x, y ≤ 1.25 m. This area of the initial flux is sufficient to cover the full crosssection surface of the generic model.
The absolute muon momenta p were treated as described in Section 2.2, while different scenarios were simulated with respect to the initial muon direction. As a first scenario (reference case), we assumed monodirectional initial muons, i.e. p⃗ = −p ⋅ê_{z} (p > 0) and θ_{in} = 0°. In addition, we restricted the incoming muon spectrum to various angular ranges with respect to θ_{in}. In particular, we considered angular distributions of 0° to 1°, 0° to 5° and 0° to 25°.
Simulations were performed for two geometries of the generic model, which differed with respect to the occupancy of the fuel compartments with fuel assemblies. In the first case, all fuel compartments were empty (no fuel assemblies). In the second case, all fuel compartments with the exception of the central one were occupied with fuel assemblies.
Projection images were generated using the transmissionradiographic as well as the scatteringradiographic approach.
In case of the transmissionradiographic analysis, the major imaging information is the ratio of muons reaching the outgoing detector over the number of muons crossing a specific pixel of the incoming detector. The pixel size was specified as (1 × 1) cm^{2}.
In case of the scatteringradiographic analysis, the leading imaging information is the effective scattering angle θ_{eff} with respect to the initial direction, which is calculated eventbyevent according to: (10)
The basis is given by the position information provided by the incoming (x_{in}, y_{in}, z_{in}) and outgoing (x_{out}, y_{out}, z_{out}) detectors as well as the normalized muon direction d⃗ = (d_{x}, d_{y}, d_{z}) provided by the incoming detector. The projected distances Δ_{i} (i =x, y, z) are given by i_{out} − i_{in}.
Information on positions and directions was processed without any attempts to mimic resolution effects. Again, the pixel size was defined as (1 × 1) cm^{2}.
3.1.1 Transmissionradiographic images
Figure 7 shows transmission radiographic images for the two geometries and the different angular distributions with respect to θ_{in} as well as associated difference images.
It can be easily seen that the image quality in terms of resolution decreases with increasing angular acceptance with respect to θ_{in}. A good indicator is given in terms of the absorber rods that can be easily identified for the angular acceptance of 0° ≤ θ_{in} ≤ 1°. This, however, is not possible for an angular acceptance of 0° ≤ θ_{in} ≤ 5°, at least not with the number of simulated events. A similar picture can be drawn for larger structures such as fuel assemblies: for the angular acceptance of 0° ≤ θ_{in} ≤ 25°, the transmissionradiographic image based on the simulated statistics does not allow for a reliable conclusion on the occupancy of the central fuel compartment with a fuel assembly.
It is obvious that difference images between the two geometries benefit from an improved contrast that allows identifying structural differences (or similarities) more easily. This, for instance, holds with respect to the occupancy of the central fuel compartment. Even for the angular acceptance of 0° ≤ θ_{in} ≤ 25° the difference image provides weak evidence allowing the conclusion that for both geometries the occupancies of the central fuel compartment were identical.
Fig. 7
The figures in the two upper rows show transmissionradiographic images of longitudinal scans of a.) an empty generic model (top row) as well as b.) a generic model filled with 18 out of 19 possible fuel assemblies (middle row) where the central fuel compartment remained empty. The bottom row shows corresponding difference images, a − b. The color code represents the ratio of muons that were detected in both detectors over the total number of muons. Different assumptions on the initial muon directions are considered from left to right and range from a fixed initial muon direction with θ_{in} = 0° (left) to an angular distribution with respect to θ_{in} of 0° to 25° (right). The x and ycoordinates refer to the muon positions at the exit of the incoming detector. 
3.1.2 Scatteringradiographic images
Figure 8 shows scatteringradiographic images for the two geometries and the various angular acceptances as well as corresponding difference images. For each pixel, the median of the associated scatteringangle distribution is used to generate the scatteringradiographic images.
Similar to the transmission radiographic images, the image quality deteriorates significantly with the increasing angular acceptance with respect to θ_{in}. However, the image quality is much better and smaller structures can be resolved. For example, the individual walls of the fuel compartments can be easily identified for the angular acceptance of 0° ≤ θ_{in} ≤ 1° and  with limitations  also for the angular acceptance of 0° ≤ θ_{in} ≤ 5°. Unfortunately, also the scatteringradiographic images would not allow for a reliable statement on the occupancy of the central fuel compartment for an angular acceptance of 0° ≤ θ_{in} ≤ 25° based on the number of simulated events.
The improved resolution compared to the transmission radiographic images is also reflected in the difference images. With respect to the angular acceptance of 0° ≤ θ_{in} ≤ 25°, at least the difference plot provides clear evidence that the occupancies of the central fuel compartment are identical for both geometries.
Fig. 8
The figures in the two upper rows show scatteringradiographic images of longitudinal scans of a.) an empty generic model (top row) as well as b.) a generic model filled with 18 out of 19 possible fuel assemblies (middle row) where the central fuel compartment remained empty. The color code represents for each pixel the median of the associated scatteringangle distribution. The bottom row shows corresponding difference images, b − a. Different assumptions on the initial muon directions are considered from left to right and range from a fixed initial muon direction with θ_{in} = 0° (left) to an angular distribution with respect to θ_{in} of 0° to 25° (right). The x and ycoordinates refer to the muon positions at the exit of the incoming detector. 
3.1.3 General aspects
The blurring of the structures with increasing angular acceptance can be easily understood in terms of the longitudinal extension of the generic model.
In general, the simulated results show the advantages of the scatteringradiographic over the transmissionradiographic approach. Here, the improved resolving power is the most prominent indicator. Unlike most radiographic detectors, tomographic detection systems would be able to apply cut conditions based on the incoming muon flight directions which would allow us to realize different angular acceptances with respect to θ_{in}.
In summary, the simulations provide evidence that even without any reconstruction efforts to generate tomographic images it would be possible to make conclusions on the occupancy of specific fuel compartments with fuel assemblies within a reasonable amount of time. For example, in case of an angular acceptance of 0° ≤ θ_{in} ≤ 25° and based on reasonable assumptions (see Sect. 3.2), the number of simulated events would correspond to a measuring time of about 40 hours. This increases to ~ 190 hours and ~ 930 hours in case of 0° ≤ θ_{in} ≤ 5° and 0° ≤ θ_{in} ≤ 1°, respectively.
So far, these results cannot be extrapolated to smaller structures such as individual fuel rods. Qualitatively, it can be expected that for such a level of detail a narrow angular acceptance as well as a much larger muon flux per area would be required. The following section provides a more detailed discussion of this aspect.
3.2 Detailed longitudinal study of the central fuel assembly to detect missing fuel rods
In this section, we investigate if muonscattering radiography can be used to make reliable statements about the completeness of individual fuel assemblies. In particular, we investigate whether individual missing fuel rods in an otherwise complete fuel assembly can be detected. This investigation takes into account a few variable boundary conditions such as the angular acceptance with respect to θ_{in} as well as the number of events, which most often can be used to provide reasonable estimates on the required irradiation time.
All simulations were performed with up to 10^{8} events. The treatment of the absolute muon momenta was identical to the former Section 3.1. The considered angular acceptances with respect to θ_{in} of the primary muons ranged from 0° ≤ θ_{in} ≤ 0.25° to 0°≤ θ_{in} ≤ 2°. As in the former section, we also considered monodirectional incoming muons with θ_{in} = 0° to provide a reference scenario. The initial positions within the x, yplane were restricted to − 0.35 m ≤ x, y ≤ 0.35 m and the initial zcoordinate was fixed at z ≈ 6.1 m.
We performed simulations for two different geometries of the generic model. The first geometry ensured that all fuel compartments of the model were occupied with fuel assemblies and for each fuel assembly all individual slot positions were occupied according to the nominative layout, see Figure 3. The second geometry differs from the first one only by three vacated fuelrod positions along the diagonal (first, fifth and ninth position – see Fig. 3) within the central fuel assembly. The first fuel rod of interest (x_{id} = y_{id} = 1 (I)) is nominatively placed in the upper left corner of the fuel assembly and is surrounded by three fuel rods and the walls of the fuel compartment. The second fuel rod of interest (x_{id} = y_{id} = 5 (II)), is surrounded by two control rods and six fuel rods, while the third fuel rod of interest (x_{id} = y_{id} = 9 (III)), is surrounded by fuel rods on all sides.
Scatteringradiographic images were generated using a binning as indicated in Figure 9. Each bin covers (6.36 × 6.36) mm^{2} and is shifted in a way that for the central fuel assembly the center positions of the bins coincide with the nominative center positions within the (x, y)plane of the various fuel and control rods.
As in the previous section, the effective scattering angle θ_{eff} is calculated event by event according to equation (10). For each pixel we obtain a histogram that shows the absolute frequencies of the effective scattering angles θ_{eff}. Representative examples of such absolute frequency distributions are shown in Figure 10.
These frequency distributions are then normalized for each pixel (i, j) into a normalized probability distribution function (PDF) ρ_{i,j}(θ_{eff}):
In a final step, we use ρ_{i,j}(θ_{eff}) to calculate for each pixel i, j the cumulative distribution function (CDF) , i.e.:
This procedure is repeated for each simulation.
For illustration, Figure 11 shows for both geometries the PDF ρ_{II} as well as the CDF corresponding to the pixel of the slot position with x_{id} = y_{id} = 5 within the central fuel assembly. The analysis focuses on a pixelbased comparison between the two geometries described above. By computing the local distances between the empirical CDFs, the derived statistical measure can be used for the pairwise comparison of images [38,39]. In other words, we compare for each pixel i, j the CDFs and . Here, corresponds to the geometry for which all slot positions within the fuel assemblies are occupied, while corresponds to the geometry for which the fuelrod slot positions within the central fuel assembly are empty.
We consider for each pixel (i, j) in the (x, y)plane the twosample KolmogorovSmirnov test, (11)
which allows us to make quantitative statements on the agreement between and .
The null hypothesis  i.e. and describe identical distributions – is rejected at a statistical significance level α if the test statistic satisfies
Here n_{i,j} (m_{i,j}) describes the number of entries in the (i, j) pixel of the relevant vacated (occupied) spectrum. The condition can be reformulated as: (12)
We specified the significance level α of 0.1 in the following.
Figure 12 shows a heatmap of based on a simulation with 10^{8} events assuming monodirectional muons with θ_{in} = 0°. One can clearly identify areas for which the KolmogorovSmirnov test statistic indicates significant deviations. These areas coincide perfectly with areas for which the two geometries differ, i.e. with respect to the occupation of the slot positions x_{id} = y_{id} = 1 (I), x_{id} = y_{id} = 5 (II) and x_{id} = y_{id} = 9 (III) within the central fuel assembly.
Figure 13 summarizes for each of these three particular pixels and for different angular acceptances and quantifies the evolution as a function of simulated events. corresponds to the pixel of x_{id} = y_{id} = 1, while and correspond to the pixels of x_{id} = y_{id} = 5 and x_{id} = y_{id} = 9, respectively.
It is obvious that for all three pixels the significance ratios decrease in general with an increasing angular acceptance. The only exception is given in terms of the angular acceptances 0° ≤ θ_{in} ≤ 1.5° and 0° ≤ θ_{in} ≤ 2°, for which comparable significance ratios are observed. For all three pixels one can observe a clear trend towards smaller slopes of progression with an increasing angular acceptance and increasing number of simulated events. It is interesting to note that there are significant differences between the three pixels with respect to the trends with increasing event numbers and the achieved significance ratios. These effects are stronger for pixel II (III) compared to I. With respect to pixel I and based on the maximum number of events (10^{8}), the significance ratio exceeds one only up to an angular acceptance of 0° ≤ θ_{in} ≤ 0.5° and requiresat least 6 ⋅ 10^{7} simulated events in case of the latter. A different picture can be drawn for pixel II. Not only do we observe that less events (4 ⋅ 10^{7} events) are required for to exceed one for 0° ≤ θ_{in} ≤ 0.5°, we also observe significant deviations with respect to 0° ≤ θ_{in} ≤ 0.75°. For the latter, at least 5 ⋅ 10^{7} simulations events are required for to exceed one. In case of pixel III, we observe that also exceeds one for 0° ≤ θ_{in} ≤ 1.0°, requiring at least 8.5 ⋅ 10^{7} simulated events. The significant deviations between the individual pixels may be related to the larger numbers of neighbouring fuel rods in case of pixels II and III compared to pixel I.
For a specific simulation with known momentum distribution and angular acceptance, the measurement time can be estimated according to the following formula: (13)
Here, A is the area of the initial (x, y)plane, I_{μ} is the integrated muon flux at sea level (~ 1 muon∕cm^{2}∕min), is the share within the angular acceptance with respect to θ_{in}, I_{p} is the share of the considered momentum distribution with respect to the full momentum distribution and ϵ is the efficiency of the detector system, which is assumed to be 1 in the following. For a certain angular acceptance, can be estimated according to: (14)
Values of for various angular acceptances are listed in Table 4. In a similar way, a simplified estimate of I_{p} (≈ 0.8) can be determined by means of equation (1): (15)
Figure 14 summarizes for each of the three pixels of interest and quantifies the evolution as a function of the estimated measurement time.
In case of an angular acceptance of 0° ≤ θ_{in} ≤ 0.5°, at least ~ 1.8 years are required for a reliable statement for pixel II. Longer measurement times are required for pixels I (~ 2.6 years) and III (~ 2.2 years). Smaller measurement times are expected for the smaller angular acceptance of 0° ≤ θ_{in} ≤ 0.25°: Here, the estimated measurement times amount to ~1.8 years (I and II) and ~0.9 years (III), respectively.Overall and based on the assumptions of the present work, the statistical benefit of a larger angular acceptance does not translate into smaller measurement times; on the contrary, larger angular acceptances lead to extended measurement times.
Fig. 9
Illustration of the binning used for the analysis in Section 3.2. The green disks represent the cross sections of the various rods in the central fuel assembly. The grid indicates the limits of the bins within the (x, y)plane. 
Fig. 10
Absolute frequency distributions of the effective scattering angle θ_{eff} for different pixels and different angular acceptances. The red and black curves correspond to the pixel of the slot position (5,5) of the central fuel assembly. The black curve represents the case where this particular slot position is empty while the red curve corresponds to the case where the slot position is filled with a fuel rod. The blue curve corresponds to a pixel between two neighbouring fuel assemblies. The top figure shows results for a fixed angular acceptance (θ_{in} = 0°) while the bottom figure shows results for an angular acceptance of 0° ≤ θ_{in} ≤ 2°. All distributions correspond to simulations with 10^{8} events. 
Fig. 11
Probability density functions ρ_{II}(θ) (top) as well as cumulative probability functions (bottom) for the pixel corresponding to the slot position with x_{id} = y_{id} = 5 within the central fuel assembly for two different geometries. The red curve represents the geometry with a complete fuel assembly (all positions are occupied according to the nominative layout), while the black curve corresponds to the geometry where the specific slot position represented by this pixel is empty. The distributions represent results of a simulation with 10^{8} events assuming monodirectional muons with θ_{in} = 0°. 
Fig. 12
A heat map showing the significance ratio in the (x, y)plane. One can easily identify the locations of the three vacated fuelrod slot positions along the diagonal. The spectrum represents simulations with 10^{8} events assuming monodirectional muons. 
Fig. 13
Significance ratios as a function of angular acceptance with respect to θ_{in} and simulated events for the three pixels of interest. The top figure shows the results of D^{I}, while D^{II} (D^{III}) is shown in the middle (bottom) figure. 
Estimations for for various angular acceptances.
Fig. 14
Significance ratios for different angular acceptances with respect to θ_{in} as a function of the estimated measurement time for the three pixels of interest. The top figure shows the results of D^{I}, while D^{II} (D^{III}) is shown in the middle (bottom) figure. 
4 Summary and conclusion
In this work we investigated the theoretical feasibility to detect individual missing fuel rods in an otherwise fully loaded cask within a reasonable amount of time using cosmic muons. We used simulations based on the Geant4 toolkit and a generic model based on the CASTOR^{®} V/19 (see Fig. 2) loaded with 18x1824 PWR fuel assemblies (see Fig. 3). The detectors above and below a vertical (standing) model were mimicked by two rectangular planes with vanishing thickness and an area of 9 m^{2}. The gap between detectors and the model was assumed to be 10 cm. The muon characteristics are described in Section 2.2. The tool can be run with a realistic muon spectrum but also allows for performing simulations assuming monoenergetic and monodirectional muons. For validation purposes, the stopping powers (see Fig. 5) and scattering angles (see Fig. 6) of the relevant target materials have been investigated and compared against established reference values from the literature. We found good agreement between simulated and literature/analytical values.
We simulated longitudinal scans that covered the full cross section of the model (fullmodel simulation). In addition, we put a special emphasis on a single fuel assembly with missing fuel rods.
The focus of the fullmodel simulation was on the effect of angular acceptance criteria on the resolution. We compared monodirectional muons (θ_{in} = 0°) with different angular acceptance ranges with respect to θ_{in}: 0° ≤ θ_{in} ≤ 1°, 0° ≤ θ_{in} ≤ 5° and 0° ≤ θ_{in} ≤ 25°. The integrated flux was kept constant with 5⋅10^{7} muons, which corresponds to a radiation time of approximately 40 hours in case of 0° ≤ θ_{in} ≤ 25°. We performed two analyses and compared the empty model with a partially loaded one and the resulting difference between the two radiographic images. The partially loaded model was lacking the central fuel assembly. In the first analysis (transmission radiographic scans) we compared the ratio of muons in the lower and upper detectors. This radiographic analysis clearly indicated the decreasing resolution with increasing angular acceptance with respect to θ_{in} of the incoming muons (see Fig. 7). The central fuel assembly could still be identified as missing based on the difference plot of the acceptance criteria 0° ≤ θ_{in} ≤ 25°. In addition, we repeated the analysis with a focus on the effective scattering angles of the muons to receive scatteringradiographic images. The image quality was much better compared to the transmissionradiographic analysis and more details could be identified (see Fig. 8). Again, for the acceptance criteria of 0° ≤ θ_{in} ≤ 25° the difference plot can be used to identify the central fuel assembly as missing.
We further investigated the capability of muonscattering radiography to provide reliable statements on the completeness of a loaded fuel assembly, assuming up to 10^{8} events and angular acceptance criteria of 0° ≤ θ_{in} ≤ 0.25°, 0° ≤ θ_{in} ≤ 0.5°, 0° ≤ θ_{in} ≤ 0.75°, 0° ≤ θ_{in} ≤ 1.0°, 0° ≤ θ_{in} ≤ 1.5° and 0° ≤ θ_{in} ≤ 2.0° in addition to monodirectional muons with θ_{in} = 0°. Focussing on the central fuel assembly we investigated a fully loaded model and one where the central assembly was missing three fuel rods on the diagonal (see Fig. 3). All three missing rods were clearly identifiable after the simulation of 10^{8} events assuming monodirectional muons with θ_{in} = 0° (see Fig. 12). The significance ratio based on the KolmogorovSmirnov test statistic of each fuel rod (and thus the significance of the contrast of a missingrod) depends on the relative position of the missing fuel rod in the assembly, the angular acceptance criteria of incoming muons and the number of simulated events (see Fig. 13). The latter can be transformed accordingly into radiation or measuring times. Assuming significance ratios leading to a successful identification of missing fuel rods we found measuring times of approximately 2 years for the two inner missing rods and ~2.6 years for the corner rod assuming an angular acceptance of 0° ≤ θ_{in} ≤ 0.5°.
In summary, we have shown that the muonscattering radiography is capable of reliably visualizing the inside of a loaded and sealed model at a resolution scale of a single fuel rod. The assumptions we made led to timescales which seem feasible compared to the duration of the dry storage of spent nuclear fuel. In the next section we discuss how future work might even shorten the theoretically estimated measuring time and how image processing methods may be used to gain more detailed insights into the cask’s interior for the detection of individual missing fuel rods.
5 Outlook
The present work will serve as a starting point for further research activities that may evolve into several independent aspects.
The first aspect relates to the developed simulation tool itself. We intend to include the discussion of transversal scans, for which the detectors are located on the sides of the generic model. The simulation can also be extended to additional geometries which do not have to be limited to storage casks.
The second aspect treats the validation of the simulated data and comparisons to other simulations. We already addressed validation aspects in the present work that concerned the slowingdown process of the muons as well as the angular scattering. A logical next step would be given by a comparison of simulated to experimental results. For this it might be reasonable to start with less complex geometries or larger objectsofinterest which will require less time, both with respect to the computation and the experimental measurement. To take into account the uncertainty due to a particular choice of the simulation tool, it would be valuable to compare the results from different simulation tools for identical geometries for benchmark purposes.
The third aspect covers to sensitivity analyses. The present results indicate that the significance depends on the relative fuel rod position and, hence, the immediate surroundings of the considered rod. So far, we have arbitrarily selected three fuel rods and it may be worth the effort to repeat the analyses systematically, including the effects of rods in the immediate surroundings. The significance information of the second analysis part is derived from two predominantly identical geometries, which deviated only by the occupation of three fuel rod slots. It would be insightful to investigate the effects of additional differences such as slight misalignments. In addition, we intend to investigate the sensitivity of the simulation results to the physics models in Geant4. The present work focuses on the modular physics list FTFP_BERT.
The fourth aspect addresses the extension to tomographic images which allows using the complete information of scattered (or absorbed) muons within one visualization processed by various potential statistical muontomography reconstruction algorithms as discussed in [40]. For adaptively comparing muonscattering (or absorption) images of a cask with the purpose of automatically detecting changes within the interior of the cask, the imageprocessing strategy needs to be extended towards tomographic image reconstruction inherently designed for change detection or in combination with further image analysis methods.
We intend to make the developed simulation code available in the near future. If interested, please contact the main author for more details.
Author contribution statement
All the authors were involved in the preparation of the manuscript. All the authors have read and approved the final manuscript.
Acknowledgements
This research was partially funded by the Federal Ministry for the Environment, Nature Conservation and Nuclear Safety in Germany (BMU) under Contract 4720E03366.
References
 K.P. Ziock, G. Caffrey, A. Lebrun, L. Forman, P. Vanier, J. Wharton, IEEE Nuclear Science Symposium Conference Record, 2005, Fajardo, 2005, pp. 1163–1167 [CrossRef] [Google Scholar]
 V. Brdar, P. Huber, J. Kopp, Phys. Rev. Appl. 8, 2331 (2017) [CrossRef] [Google Scholar]
 G. Bonomi, P. Checchia, M. D’Errico, D. Pagano, G. Saracino, Prog. Part. Nucl. Phys. 112, 103768 (2020) [CrossRef] [Google Scholar]
 P.A. Zyla et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2020, 083C01 (2020) [CrossRef] [Google Scholar]
 M.P. De Pascale et al., J. Geophys. Res. A 98, 3501 (1993) [CrossRef] [Google Scholar]
 P.K.F. Grieder, Cosmic Rays at Earth (Elsevier Science, Amsterdam, 2001) [Google Scholar]
 L.W. Alvarez et al., Science 167, 832–839 (1970) [CrossRef] [PubMed] [Google Scholar]
 H.K.M. Tanaka, T. Kusagaya, H. Shinohara, Nat. Commun. 5, 3381 (2014) [CrossRef] [Google Scholar]
 N. Lesparre, D. Gilbert, J. Marteau, Y. Déclais, D. Carbone, E. Galichet, Geophys. J. Int. 183, 1348–1361 (2010) [CrossRef] [Google Scholar]
 K. Morishima et al., Nature 552, 386–390 (2017) [CrossRef] [PubMed] [Google Scholar]
 H.K.M. Tanaka, K. Nagamine, S.N. Nakamura, K. Ishida, Nucl. Instrum. Methods Phys. Res. A 555, 164–172 (2005) [CrossRef] [Google Scholar]
 K. Borozdin et al., Nature 422, 277 (2003) [CrossRef] [Google Scholar]
 A. Clarkson et al., J. Instrum. 10, P03020 (2015) [CrossRef] [Google Scholar]
 Decision Sciences https://www.decisionsciences.com/ourproduct/. [Google Scholar]
 A. Zenoni. 2015 4th International Conference on Advancements in Nuclear Instrumentation Measurement Methods and their Applications (ANIMMA) (IEEE, Lisbon, 2015) [Google Scholar]
 J.M. Durham et al., Phys. Rev. Appl. 9, 044013 (2018) [CrossRef] [Google Scholar]
 G. Jonkmans, V.N.P. Anghel, C. Jewett, M. Thompson, Ann. Nucl. Energy 53, 267 (2013) [CrossRef] [Google Scholar]
 A. Clarkson et al., Nucl. Instrum. Methods Phys. Res. A 746, 64–73 (2014) [CrossRef] [Google Scholar]
 S. Chatzidakis, M. Alamaniotis, L.H. Tsoukalas, Trans. Am. Nucl. Soc. 111, 369 (2014) [Google Scholar]
 F. Ambrosino et al., J. Instrum. 10, T06005 (2015) [CrossRef] [Google Scholar]
 D. Poulson, J.M. Durham, E. Guardincerri, C.L. Morris, J.D. Bacon, K. PlaudRamos, D. Morley, A.A. Hecht, Nucl. Instrum. Methods Phys. Res. A 842, 48 (2017) [CrossRef] [Google Scholar]
 D. Poulson, J. Bacon, M. Durham, E. Guardincerri, C.L. Morris, H.R. Trellue, Phil. Trans. Roy. Soc. A 377, 2137 (2018) [Google Scholar]
 S. Agostinelli et al., Nucl. Instrum. Methods Phys. Res. A 506, 250–303 (2003) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 J. Allison et al., IEEE Trans. Nucl. Sci. 53, 270–278 (2006) [Google Scholar]
 J. Allison et al., Nucl. Instrum. Methods Phys. Res. A 835, 186 (2016) [NASA ADS] [CrossRef] [Google Scholar]
 R. Brun, F. Rademakers, Nucl. Instrum. Methods Phys. Res. A 389, 81–86 (1997) [NASA ADS] [CrossRef] [Google Scholar]
 See also “ROOT” [software], Release v6.20/04, 01/04/2020 [Google Scholar]
 GNS Geselllschaft für NuklearSerivce mbH, Product Info Castor^{®} V/19 [Google Scholar]
 Bundesamt für Strahlenschutz (BfS), Radioaktive Frachten unterwegs  Atomtransporte und Sicherheit (brochure, 2000) [Google Scholar]
 CMS Collaboration, Phys. Lett. B 692, 83 (2010) [CrossRef] [Google Scholar]
 D. Reyna, arXiv:hepph/0604145v2 [Google Scholar]
 E.V. Bugaev, A. Misaki, V.A. Naumov, T. Sinegovskaya, S.I. Sinegovsky, N. Takahashi, Phys. Rev. D 58, 054001 (1998) [CrossRef] [Google Scholar]
 J.W. Eaton, D. Bateman, S. Hauberg, R. Wehbring, 2016 GNU Octave version 5.2.0 manual: a highlevel interactive language for numerical computations. https://www.gnu.org/software/octave/doc/v.5.2.0/ [Google Scholar]
 D.E. Groom, N.V. Mokhov, S.I. Striganov, At. Data Nucl. Data Tables 78, 183 (2001) [CrossRef] [Google Scholar]
 https://pdg.lbl.gov/2020/AtomicNuclearProperties/ [Google Scholar]
 W.H. Bragg, R. Kleeman, Philos. Mag. 10, 318–340 (1905) [CrossRef] [Google Scholar]
 G.R. Lynch, O.I. Dahl, Nucl. Instrum. Methods Phys. Res. B 58, 6 (1991) [CrossRef] [Google Scholar]
 E. Demidenko, in: Computational Science and Its Applications  ICCSA 2004. ICCSA 2004. Lecture Notes in Computer Science, edited by A. Laganá, M.L. Gavrilova, V. Kumar, Y. Mun, C.J.K. Tan, O. Gervasi (SpringerVerlag Berlin, Heidelberg, 2004), vol. 3046 [Google Scholar]
 Y. Tang, L. Zhang, X. Huang, Int. J. Remote Sens. 32, 5719 (2011) [CrossRef] [Google Scholar]
 S. Riggi et. al., Nucl. Instrum. Methods Phys. Res. A 728, 59 (2013) [CrossRef] [Google Scholar]
Cite this article as: Thomas Braunroth, Nadine Berner, Florian Rowold, Marc Péridis, Maik Stuke, Muon radiography to visualise individual fuel rods in sealed casks, EPJ Nuclear Sci. Technol. 7, 12 (2021)
All Tables
Comparison of key dimensions between the generic model and the CASTOR^{®} V/19 cask (in storage configuration).
All Figures
Fig. 1
Visualisation of the generic model and its upward orientation within the present work. The coordinate system as it is used in this study is provided in the lower right corner. The origin of the coordinate system coincides with the center position of the bottom part of the model. The red and green areas above and below the model indicate the incoming and outgoing detectors. 

In the text 
Fig. 2
Exploded view of the generic model showing its individual components: monolithic body (1), basket with 19 fuel compartments (2), trunnions (3), primary lid (4), polyethylene plate (5), secondary lid (6), protection plate (7), inner moderator rods and plugs (8), outer moderator rods and plugs (9), polyethylene plate (10), base plate (11) and a representative fuel assembly (12). See text for details. 

In the text 
Fig. 3
Left: topview on the arrangement of fuel rods (yellow), cladding tubes (green) as well as control rods (blue) in case of a complete fuel assembly. Each element can be specified unambiguously based on its slot position given by x_{id}, y_{id}. Right: same as left figure, but several elements along the diagonal are missing (x_{id} ∕y_{id}= 1, 5, 9). This configuration is investigated in more detail in Section 3.2. 

In the text 
Fig. 4
Sample distributions of the initial angle θ_{in} (left) and of the absolute muon momentum p (right) as determined with the incoming detector. The red curve shown in the left spectrum is proportional to , while the red curve shown in the right spectrum is proportional to equation (1) with θ = 12.5°. 

In the text 
Fig. 5
Comparison of simulated mean energy losses ΔE per distance and mean differential energy losses dE/dx (both in [MeV/mm]) of muons in various target materials. The top figure shows simulated mean energy losses (symbols) and compares them to tabulated reference values (solid lines). The red curve (boxes) shows the results for uranium dioxide, theorange curve (circles) corresponds to polyethylene, the blue curve (pentagons) corresponds to stainless steel, the green curve (uppointing triangles) corresponds to ductile iron and the violet curve (downpointing triangles) corresponds to zirconium alloy. The lower five figures show for each material the relative deviations between simulated energy loss per mm and tabulated mean stopping power values. 

In the text 
Fig. 6
Comparison of rootmean square (RMS) values of scattering angles (in [mrad]) between semiempirical calculations and simulated values for muons in various target materials with a thickness of 1 mm. The top figure shows values derived from the simulation (symbols) and compares them to calculated values based on a semiempirical approach (solid lines). The red curve (boxes) shows the results for uranium dioxide, the orange curve (circles) corresponds to polyethylene, the blue curve (pentagons) corresponds to stainless steel, the green curve (uppointing triangles) corresponds to ductile iron and the violet curve (downpointing triangles) corresponds to zirconium alloy. The lower five figures show for each material therelative deviations between simulated and calculated values. 

In the text 
Fig. 7
The figures in the two upper rows show transmissionradiographic images of longitudinal scans of a.) an empty generic model (top row) as well as b.) a generic model filled with 18 out of 19 possible fuel assemblies (middle row) where the central fuel compartment remained empty. The bottom row shows corresponding difference images, a − b. The color code represents the ratio of muons that were detected in both detectors over the total number of muons. Different assumptions on the initial muon directions are considered from left to right and range from a fixed initial muon direction with θ_{in} = 0° (left) to an angular distribution with respect to θ_{in} of 0° to 25° (right). The x and ycoordinates refer to the muon positions at the exit of the incoming detector. 

In the text 
Fig. 8
The figures in the two upper rows show scatteringradiographic images of longitudinal scans of a.) an empty generic model (top row) as well as b.) a generic model filled with 18 out of 19 possible fuel assemblies (middle row) where the central fuel compartment remained empty. The color code represents for each pixel the median of the associated scatteringangle distribution. The bottom row shows corresponding difference images, b − a. Different assumptions on the initial muon directions are considered from left to right and range from a fixed initial muon direction with θ_{in} = 0° (left) to an angular distribution with respect to θ_{in} of 0° to 25° (right). The x and ycoordinates refer to the muon positions at the exit of the incoming detector. 

In the text 
Fig. 9
Illustration of the binning used for the analysis in Section 3.2. The green disks represent the cross sections of the various rods in the central fuel assembly. The grid indicates the limits of the bins within the (x, y)plane. 

In the text 
Fig. 10
Absolute frequency distributions of the effective scattering angle θ_{eff} for different pixels and different angular acceptances. The red and black curves correspond to the pixel of the slot position (5,5) of the central fuel assembly. The black curve represents the case where this particular slot position is empty while the red curve corresponds to the case where the slot position is filled with a fuel rod. The blue curve corresponds to a pixel between two neighbouring fuel assemblies. The top figure shows results for a fixed angular acceptance (θ_{in} = 0°) while the bottom figure shows results for an angular acceptance of 0° ≤ θ_{in} ≤ 2°. All distributions correspond to simulations with 10^{8} events. 

In the text 
Fig. 11
Probability density functions ρ_{II}(θ) (top) as well as cumulative probability functions (bottom) for the pixel corresponding to the slot position with x_{id} = y_{id} = 5 within the central fuel assembly for two different geometries. The red curve represents the geometry with a complete fuel assembly (all positions are occupied according to the nominative layout), while the black curve corresponds to the geometry where the specific slot position represented by this pixel is empty. The distributions represent results of a simulation with 10^{8} events assuming monodirectional muons with θ_{in} = 0°. 

In the text 
Fig. 12
A heat map showing the significance ratio in the (x, y)plane. One can easily identify the locations of the three vacated fuelrod slot positions along the diagonal. The spectrum represents simulations with 10^{8} events assuming monodirectional muons. 

In the text 
Fig. 13
Significance ratios as a function of angular acceptance with respect to θ_{in} and simulated events for the three pixels of interest. The top figure shows the results of D^{I}, while D^{II} (D^{III}) is shown in the middle (bottom) figure. 

In the text 
Fig. 14
Significance ratios for different angular acceptances with respect to θ_{in} as a function of the estimated measurement time for the three pixels of interest. The top figure shows the results of D^{I}, while D^{II} (D^{III}) is shown in the middle (bottom) figure. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.