© T. Braunroth et al., Published by EDP Sciences, 2021

1 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⁻² s⁻¹ sr⁻¹ [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⁻² min⁻¹ 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 non-invasive 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 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³ 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 reference [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 scattering-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 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 Monte-Carlo 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 event-by-event into ROOT [26,27] container files, which allows performing post-processing 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 bottom-nozzle, fuel rods as well as control rods. The fuel rods consist of nuclear fuel (UO₂) 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 Figure 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.

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.

Fig. 3 Left: top-view 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 xid, yid. Right: same as left figure, but several elements along the diagonal are missing (xid ∕yid= 1, 5, 9). This configuration is investigated in more detail in Section 3.2.

2.1.2 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, which is particularly useful for validation purposes.

2.1.3 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 Figure 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.

2.1.4 Detectors

Detector systems are mimicked by two rectangular detector planes (vanishing thickness, area of (3 × 3) m²), 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 event-by-event base, it benefits from an infinite resolving power.¹ In this regard, it provides a best-case scenario and can be considered as a first important theoretical step towards substantiated feasibility studies.

2.2 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 () to determine for each event k the particle properties, i.e. particle type (muon (μ⁻) or anti-muon (μ⁺)), 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 anti-muons, 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 anti-muon (μ⁺), 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₁ to c₅ 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 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 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 distribution-based approach described above, the tool also allows performing simulations with mono-energetic and mono-directional 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 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 () 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 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.

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 two-detector 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 energy-spectra are below 2%. The only exception is given for polyethylene, for which the limit is increased to 5%. In total, between 10⁵ and 5 ⋅ 10⁶ 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 (up-pointing triangles) corresponds to ductile iron and the violet curve (down-pointing 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 semi-empirical description. Each material had a thickness of 1 mm.

According to reference [37], the root-mean square of the scattering angle in the x- and y-planes 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₀ is the radiation length of the material. The root-mean square of the scattering angle θ^S in space is given by: (8)

In case of a compound material, the radiation length X₀ 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₀ 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⁵ 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 semi-empirical values for are usually less than 3%. Only for UO₂ 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.

Fig. 6 Comparison of root-mean square (RMS) values of scattering angles (in [mrad]) between semi-empirical 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 semi-empirical 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 (up-pointing triangles) corresponds to ductile iron and the violet curve (down-pointing 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 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 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⁷ 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 Section 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 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 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².

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: (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².

3.1.1 Transmission-radiographic 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 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.

Fig. 7 The figures in the two upper rows show transmission-radiographic 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 y-coordinates refer to the muon positions at the exit of the incoming detector.

3.1.2 Scattering-radiographic images

Figure 8 shows scattering-radiographic images for the two geometries and the various angular acceptances as well as corresponding difference images. For each pixel, the median of the associated scattering-angle distribution is used to generate the scattering-radiographic 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.

Fig. 8 The figures in the two upper rows show scattering-radiographic 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 scattering-angle 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 y-coordinates 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 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 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 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⁸ 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 Figure 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 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.

Scattering-radiographic images were generated using a binning as indicated in Figure 9. Each bin covers (6.36 × 6.36) mm² 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 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 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 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, (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⁸ events assuming mono-directional muons with θ_in = 0°. One can clearly identify areas for which the Kolmogorov-Smirnov 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⁸), the significance ratio exceeds one only up to an angular acceptance of 0° ≤ θ_in ≤ 0.5° and requiresat least 6 ⋅ 10⁷ 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⁷ 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⁷ 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⁷ 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²∕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 108 events.

Fig. 11 Probability density functions ρII(θ) (top) as well as cumulative probability functions (bottom) for the pixel corresponding to the slot position with xid = yid = 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 108 events assuming mono-directional 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 fuel-rod slot positions along the diagonal. The spectrum represents simulations with 108 events assuming mono-directional 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 DI, while DII (DIII) is shown in the middle (bottom) figure.

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 DI, while DII (DIII) 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 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². 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 mono-energetic 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⁷ 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⁸ 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⁸ events assuming mono-directional muons with θ_in = 0° (see Fig. 12). The significance ratio based on the Kolmogorov-Smirnov 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 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.

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 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.

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 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.

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

All Tables

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: top-view 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 xid, yid. Right: same as left figure, but several elements along the diagonal are missing (xid ∕yid= 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 (up-pointing triangles) corresponds to ductile iron and the violet curve (down-pointing 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 root-mean square (RMS) values of scattering angles (in [mrad]) between semi-empirical 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 semi-empirical 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 (up-pointing triangles) corresponds to ductile iron and the violet curve (down-pointing 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 transmission-radiographic 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 y-coordinates 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 scattering-radiographic 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 scattering-angle 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 y-coordinates 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 108 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 xid = yid = 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 108 events assuming mono-directional 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 fuel-rod slot positions along the diagonal. The spectrum represents simulations with 108 events assuming mono-directional 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 DI, while DII (DIII) 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 DI, while DII (DIII) is shown in the middle (bottom) figure.
In the text