Muon Radiography to Visualise Individual Fuel Rods in Sealed Casks

Cosmic-ray muons can be used for the non-destructive 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 Kolmogorov-Smirnov 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.


Introduction
The operation of nuclear power plants generates high-level radioactive wastes which need to be stored and disposed of. A well-established 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 high-level 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 three-dimensional 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 anti-muons (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 non-invasive 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 two-dimensional 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 three-dimensional 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, Monte-Carlo 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 stor-age 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 muon scattering radiography to a Westinghouse MC-10 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 high-Z. Clarkson et al. [18] performed Geant4 simulations of a scintillating-fibre 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 one-month period of measurement. Poulsen et al. [21] were the first to apply filtered back projection 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 Ref. [16] for a numerical study using Geant4 to distinguish different loads of a cask. The study indicates that a one-week 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 Monte-Carlo 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 scatter-angle 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 Sec. 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. Sec. 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 Sec. 4. We end with a short outlook in Sec. 5.

Simulation and Data Processing
Simulations were performed with a dedicated tool based on the Monte-Carlo toolkit Geant4 [23,24,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 (Sec. 2.1), the treatment of primary particle properties (Sec. 2.2), aspects related to physics and tracking (Sec. 2.3) as well as optimization strategies (Sec. 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 eventby-event into ROOT [26,27] container files, which allows performing post-processing in a flexible manner. Finally, Sec. 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.

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 Refs. [28,29]. A visualisation of the generic model was generated with the Geant4 OpenGL interface and is depicted in Fig. 1. The individual components of the generic model can be identified in the exploded view shown in Fig. 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 bottom-nozzle, fuel rods as well as control rods. The fuel rods consist of nuclear fuel (UO 2 ) and cladding tubes (Zirconium alloy). A complete 18x18-24 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 Fig. 3. The modelled top and base components of the assembly are simplified.
They are assumed to be box-like, 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 command-based 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.

Box-Like Object
Instead of the generic model, it is possible to generate a box-like 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, particularly useful for validation purposes.

Coordinate System
The z-axis coincides with the symmetry axis of the generic model and is oriented upwards, i.e. its orientation is selected so that the z-coordinate of the model's top is larger  than the z-coordinate 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 right-handed euclidean space. The orientation of the individual axes is indicated in Fig. 1. The angle θ reflects the angle of a given vector r and the inverted z-axis. 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 .

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  (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.
≈ 10 cm. Muons traversing these planes are tracked and key information is determined, see Sec. 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 event-by-event base, it benefits from an infinite resolving power 1 . In this regard, it provides a best-case scenario and can be considered as a first important theoretical step towards substantiated feasibility studies.

Primary Particles
Primaries are generated by a user-defined class based on the G4VUserPrimaryGeneratorAction class provided by Geant4.
The Monte-Carlo approach is realized by using uniform random number generators with limits a, b (U[a, b]) to determine for each event k the particle properties, i.e. particle type (muon (µ − ) or anti-muon (µ + )), initial momen-tum information (p µ k , θ k , ϕ k ) as well as initial position information (x k , y k ). Details are described in the following.

Particle Type
All primaries are muons or anti-muons, assuming a charge ratio µ + /µ − equal to 1.28 [30]. The particle type for the event k is determined assuming U[0, 2.28] -if the generated random number is smaller than (or equal to) 1.28, the generated primary for this event will be an anti-muon (µ + ), otherwise it will be a muon (µ − ).

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 Here, I V is given by the Bugaev parametrisation [32]: The parameters c 1 to c 5 were determined in Ref. [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 two-dimensional 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 contributions of trajectories that geometrically can only be detected by one of the detectors. All b θ,i span over identical angular ranges. For each angular bin b θ,i , the lower and upper limits are then given by θ ll ., n, can then be associated to each of these bins of the angular spectrum, where θ i 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 : Information on θ ll i , θ ul i and w i are stored in a dedicated text file that is used as an input parameter to the simulation. U[0, i w i ] 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 U[θ ll k , θ ul k ]. The azimuthal angle ϕ k is determined for each event k using U[−π, π]. , 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 Sec. 3.2. 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 b pµ,j (j = 1, ..., m) with increasing bin sizes. Each bin b pµ,j is characterized by the following integral: Here, p ll µ,j and p ul µ,j denote the lower and upper limits of the bin b pµ,j . For each angular bin b θ,i , values for p ll µ,j , p ul µ,j and v i,j are stored in dedicated text-files which are used as mandatory input information for the simulation code. U[0, j v i,j ] uses these weights v i,j to determine for each event k the proper momentum bin b pµ,k . Finally, the proper p µ,k is determined using U[p ll µ,k , p ul µ,k ]. Within the present work, all numerical integrations were performed with GNU Octave [33]. Fig. 4 shows histograms of the simulated angular and momentum distributions for θ min = 0 • , θ max = 25 • , n = 10 and m = 172. In addition to this distribution-based approach described above, the tool also allows performing simulations with mono-energetic and mono-directional muons.

Initial Particle Position
The initial positions of the muons are selected as follows. The z-coordinate is a fixed value and selected so that it is ensured that the muon is created just above the incoming detector. The x-and y-coordinates are selected from uniform distributions covering ranges from x min to x max (U[x min , x max ]) and y min to y max (U[y min , y max ]), respectively. The associated limits can be specified in the input file of the simulation.

Physics and Tracking
This code uses the modular physics list FTFP BERT implemented in Geant4, which is recommended by the Geant4 developers for high-energy 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.

Cutting Conditions to Reduce Computational Time
Muons may interact with the material 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. photons or 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 Sec. 2.2, the muon properties -with respect to the angle θ and kinetic energy T -were cut to avoid the simulation of muon trajectories that are not useful for an analysis based on a two-detector setup. Finally, all muon trajectories are annihilated at the exit of the outgoing detector.

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 Sec. 2.5.1. The second aspect deals with the angular scattering of muons after traversing matter with a known thickness Sec. 2.5.2. In both cases, the referenced thickness was specified to 1 mm.

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 alloyto 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 polyethylenetabulated values. Hence, reference values were calculated for the other compounds -stainless steel, ductile iron and zirconium alloy -according to Bragg 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 energy-spectra 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 Fig. 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 between T = 80 MeV and T = 10 GeV is very good. However, with increasing kinetic energies 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 simulation appears to slightly overestimate the stopping power of muons in the investigated materials.

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 semi-empirical description. Each material had a thickness of 1 mm. According to Ref. [37], the root-mean square of the scattering angle θ S px,y in the x-and y-planes can be described using the following formula: σ θ S px,y = 13.6 MeV β·p·c ∆z X 0 · 1 + 0.038 · ln ∆z X 0 ·β 2 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 root-mean square of the scattering angle θ S in space is given by: In case of a compound material, the radiation length X 0 can be calculated using the formula 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 Refs. [34,35] or calculated according to Eq. 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 Fig. 6. The deviations between simulated and semi-empirical σ θ S 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.

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 Sec. 3.1.
The second section provides a detailed study that focuses on the central fuel assembly. In particular, it investigates the suitability of muon scattering radiography to make statements about the occupancy of individual fuel-rod slots within a certain fuel assembly. This section focuses on quantitative aspects and is described in Sec. 3.2 in more detail.

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 y-coordinates 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 cross-section surface of the generic model. The absolute muon momenta p were treated as described in Sec 2.2, while different scenarios were simulated with respect to the initial muon direction. As a first scenario (reference case), we assumed mono-directional 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 assem-blies). 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 transmission radiographic as well as the scattering radiographic approach.
In case of the transmission radiographic 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 scattering radiographic analysis, the leading imaging information is the effective scattering angle θ eff with respect to the initial direction, which is calculated event-by-event according to: 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 .
Transmission Radiographic Images Fig. 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 transmission radiographic 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.
Scattering Radiographic Images Fig. 8 shows scattering radiographic images for the two geometries and the various angular acceptances as well as corresponding difference 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 scattering radiographic 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.

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 scattering radiographic over the transmission radiographic 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 Sec. 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.

Detailed Longitudinal Study of the Central Fuel Assembly to Detect Missing Fuel Rods
In this section, we investigate if muon scattering 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 mono-directional incoming muons with θ in = 0 • to provide a reference scenario. The initial positions within the x, y-plane were restricted to −0.35 m ≤ x, y ≤ 0.35 m and the initial z-coordinate 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 Fig. 3. The second geometry differs from the first one only by three vacated fuel-rod positions along the diagonal (first, fifth and ninth position -see Fig. 3) within the central fuel assembly. The first relevant fuel rod I (x id = y id = 1) 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 relevant fuel rod II (x id = y id = 5) is surrounded by two control rods and six fuel rods, while the third relevant fuel rod III (x id = y id = 9) is surrounded by fuel rods on all sides. Scattering radiographic images were generated using a binning as indicated in Fig. 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 last section, the effective scattering angle θ eff is calculated event by event according to Eq. 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 Fig. 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) F Θi,j (θ eff ), i.e.: This procedure is repeated for each simulation. For illustration, Fig. 11 shows for both geometries the PDF ρ II as well as the CDF F Θ II corresponding to the pixel of the slot position with x id = y id = 5 within the central fuel assembly. The analysis focuses on a pixel-based 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 pair-wise comparison of images [38,39]. In other words, we compare for each pixel i, j the CDFs F occ Θi,j (θ eff ) and F vac Θi,j (θ eff ). Here, F occ Θi,j (θ eff ) corresponds to the geometry for which all slot positions within the fuel assemblies are occupied, while F vac Θi,j (θ eff ) corresponds to the geometry for which the fuel-rod slot positions within the central fuel assembly are empty. We consider for each pixel (i, j) in the (x, y)-plane the two-sample Kolmogorov-Smirnov test,  which allows us to make quantitative statements on the agreement between F occ Θi,j (θ eff ) and F vac Θi,j (θ eff ). The null hypothesis -i.e. F occ Θi,j (θ eff ) and F vac Θi,j (θ eff ) 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: We specified the significance level α of 0.1 in the following. Fig. 12 shows a heatmap ofD i,j based on a simulation with 10 8 events assuming mono-directional muons with θ in = 0 • . One can clearly identify areas for which the Kolmogorov-Smirnov test statistic indicates the 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. Fig. 13 summarizesD i,j for each of these three particular pixels and for different angular acceptances and quantifies the evolution as a function of simulated events.D I corresponds to the pixel of x id = y id = 1, whileD II andD III correspond to the pixels of x id = y id = 5 and x id = y id = 9. It is obvious that for all three pixels the significance ratios D 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 significant ratioD I exceeds one only up to an angular acceptance of 0 • ≤ θ in ≤ 0.5 • and requires at 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 D II 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 forD II to exceed one. In case of pixel III, we observe thatD III 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: 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), I θin 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. For a certain angular acceptance, I θin can be estimated according to: Values of I θin for various angular acceptances are listed in Tab. 4. In a similar way, a simplified estimate of I p (≈ 0.8) can be determined by means of Eq. 1: p(T =1 GeV) dp I(p µ , 0) ∞ 0 dp I(p µ , 0) With respect to 0 • ≤ θ in ≤ 0.5 • ,D I exceeds 1.0 for 6 · 10 7 events. Based on the above estimates, this event number corresponds to a measurement time of ∼ 2.5 years. Slightly fewer events (4 · 10 7 and 5 · 10 7 ) are required for the pixels II and III, which leads to estimated measurement times of ∼ 1.7 years and ∼ 2.1 years.

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 18x18-24 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 mono-directional 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 (full-model simulation). In addition, we put a special emphasis on a single fuel assembly with missing fuel rods. The focus of the full-model simulation was on the effect of angular acceptance criteria on the resolution. We compared mono-directional 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 scattering radiographic images. The image quality was much better compared to the transmission radiographic 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 muon scattering 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 mono-directional 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 mono-directional muons with θ in = 0 • (see Fig. 12). The significance ratioD based on the Kolmogorov-Smirnov test statistic of each fuel rod (and thus the significance of the contrast of a missing rod) 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 ratiosD i,j > 1 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.5 years for the corner rod assuming an angular acceptance of 0 • ≤ θ in ≤ 0.5 • . In summary, we have shown that the muon scattering 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.

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 relates to the validation of the simulated data and comparisons to other simulations. We already addressed validation aspects in the present work that concerned the slowing-down 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 objects-of-interest 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. A third aspect relates 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. A 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 muon tomography reconstruction algorithms as discussed in [40]. For adaptively comparing muon scattering (or absorption) images of a cask with the purpose of automatically detecting changes within the interior of the cask, the image processing strategy needs to be extended towards tomographic image reconstruction inherently designed for change detection or in combination with further image analysis methods.