Momentum informed muon scattering tomography for monitoring spent nuclear fuels in dry storage cask

Development of an effective monitoring method for spent nuclear fuel (SNF) in a dry storage cask (DSC) is important to meet the increasing demand for dry storage investigations. The DSC investigation should provide information about the quantity of stored SNF, and quality assurance of materials should be possible without opening the cask. However, traditional nondestructive examination (NDE) methods such as x-rays are difficult to deploy for DSC investigation because a typical DSC is intentionally designed to shield against radiation. To address this challenge, cosmic ray muons (CRMs) are used as an alternative NDE radiation probe because they can easily penetrate an entire DSC system; however, a wide application of muons is often hindered due to the naturally low CRM flux (~104 muons/m2/min). This paper introduces a newly proposed imaging algorithm, momentum-informed muon scattering tomography (MMST), and presents how a limitation of the current muon scattering tomography technique has been addressed by measuring muon momentum. To demonstrate its functionality, a commercial DSC with 24 pressurized light water reactor fuel assemblies (FAs) and the MMST system were designed in GEANT4. Three noticeable improvements were observed for MMST system as a DSC investigation tool: (1) a signal stabilization, (2) an enhanced capability to differentiate various materials, and (3) statistically increased precision to identify and locate missing FAs. The results show that MMST improves the investigation accuracy from 79 to 98% when one FA is missing and 51% to 88% when one-half FA is missing. The advancement of the NDE technique using CRM for DSC verification is expected to resolve long-standing problems in increasing demand for DSC inspections and nuclear security.

direction with measurable angles because of consecutive Coulomb scatterings with nuclei and electrons.These unique features of CRMs enabled scientists to devise a new imaging technique, muon scattering tomography (MST) [13][14][15] .For decades, CRMs have been widely used for many engineering applications including monitoring SNF assemblies in DSCs [16][17][18][19][20][21][22][23][24] .Although the results are promising, the applications of muon tomography are often limited because of CRMs' wide energy spectrum and angular distribution.Therefore, without measuring the energy and incoming direction of CRMs, possibilities would be limited with respect to improving muon tomography image resolution.To measure muon energy in the field for muon tomography applications, a few technical approaches have been developed by using time-of-flight 25,26 , scattering angle analysis 27,28 , and Cherenkov radiation 29,30 .In this work, we use a Cherenkov muon spectrometry's energy measurement capability because it is field-deployable, technically simple, and easily coupled with muon detectors [31][32][33] .
Based on the measurement capability of scattering angle and energy analysis offered by muon detectors and spectrometries, we developed a new imaging algorithm that includes muon energy information in the existing MST imaging algorithm, called momentum-informed muon scattering tomography (MMST).Both MST and MMST use a point-of-closest-approach algorithm to locate scattering positions 34,35 .However, MST uses scattering angles for density mapping, whereas MMST uses M-values, which encode momentum and scattering angle parameters in a single variable.The MMST algorithm does not increase computational costs because an M-value is an independent variable in the objective function that simply replaces an old variable.To demonstrate the benefits of the new imaging algorithm, a DSC with 24 pressurized light water reactor (PWR) SNF assemblies was modeled, and CRMs were generated using a high-energy particle transport code, GEANT4 36,37 .Then, the computational simulation results of muon interactions with the DSC and SNFs were analyzed using MATLAB.Visual and statistical inspections with MST and MMST algorithms are performed for surveying the DSC when, one-half, one, and two fuel assemblies (FAs) are missing out of 24 FAs.
The objectives of this work are to (1) demonstrate the advantages of measuring muon energy in MST applications, (2) verify the functionality of the newly developed imaging algorithm, MMST, and (3) discuss the improvement of statistical reliability and image resolution when MMST is used for DSC investigations.The overall image quality is improved by adapting MMST algorithm and it enables to visually locate one and one-half missing fuel assemblies in the DSC.From the signal analysis results of scattering angles and M-values, three improvements were observed: (1) scattering angle and M-value signals are stabilized, (2) density mapping of structural materials and SNFs becomes distinctive, and (3) the detection capability of finding missing FAs in DSC is improved.The statistical signal analysis results show that the MMST can find and locate a missing FA with a 98% confidence level (CL) compared with 79% for MST.Additionally, the MMST algorithm successfully locates a position of a missing one-half FA with an 88% CL, which was challenging with MST because of its low CL (51%).

Cosmic ray muon
Except for neutrinos and protons, muons are the most abundant cosmic particles at sea level.Cosmic ray muons are known as secondary cosmic particles because they are produced by pion and kaon decays during extensive air showers caused by interactions between high-energy primary cosmic particle (i.e., protons, electrons, and heavy ions such as O + and N + ) and air molecules in the Earth's atmosphere 38 .Muons are similar to electrons, but they are approximately 207 times heavier and decay to electrons and neutrinos with a mean lifetime of 2.2 μs 39 .Although the best estimate of CRM flux at sea level is approximately 70 m −2 s −1 sr −1 , it varies depending on various geometrical factors such as zenithal angle, azimuthal angle, altitude, longitude, latitude, and solar activities 38 .

Cosmic ray muon energy spectrum
The muon energy spectrum was derived by exploiting three phenomena-production spectra from kaons and pions, energy loss in the atmosphere, and decay.The differential CRM energy spectrum at sea level is given by 38 : where dN μ /dE μ dΩ is the expected muon counts with respect to muon energy per solid angle, and φ is the muon zenith angle.The first and second terms on the right-hand side represent contributions of pions and kaons, respectively.The experiment results of CRM measurements in the muon momentum range of 0.1-10 GeV/c at φ = 0° are shown in Fig. 1.

Cosmic ray muon angular distribution
The angular distribution is another important CRM profile with the energy spectrum because it varies depending on the incoming flight direction of muons or the zenith angle.A cosine-squared model is widely used because of its simpleness and accuracy.
where I 0 is the vertical muon flux (φ = 0°).A new approximation model was developed to address two limitations of the cosine-squared model, (1) a non-zero muon flux at φ = 90° problem and (2) a point detector assumption.Therefore, a muon zenith angle distribution can be written as where (1)

Muon scattering tomography
When a muon interacts with matter, it loses energy and is deflected mainly due to the inelastic collisions with electrons and consecutive collisions with nuclei, or multiple Coulomb scattering (MCS) 38 .The muon deflection angle due to a single Coulomb scattering is described using the Rutherford formula: where dσ/dΩ is the scattering cross section into a unit solid angle, or differential scattering cross section, θ is the muon scattering angle, z is the atomic number of a target, e is the electron charge, ϵ 0 is the permittivity of the vacuum, and E k is the initial kinetic energy for the projectile (muon) in MeV.Unless a target object is extremely thin, a muon will undergo MCS and the sum of multiple random deflection angles can be approximated with a Gaussian distribution 41 .The number of scattered muons between θ to θ + dθ, can be approximated using a 2D radial Gaussian distribution, which is given by where θ is the muon scattering angle, and σ θ is the standard deviation.σ θ can be estimated by material properties and muon momentum, and it is given by 42 where X is the thickness of material, X 0 is the radiation length, p is the muon momentum, and βc is the muon velocity in terms of the speed of light, c.Radiation length, X 0 , is a nuclear property, and it can be found in the Particle Data Group library 43,44 .For composite and unconventional material, the X 0 value can be approximated within < 1% error by the following analytical formula 45 : where A is the atomic number. (

Position mapping method
To reconstruct images of objects using MST, muon scatterings must be located for position mapping.However, it is almost impossible to reconstruct the actual scattering trajectory of muons in a medium because the measured scattering angle is a result of consecutive random small deflections, although some studies successfully reduce uncertainty in the muon trajectory estimation using Bayesian analysis 46 and maximum likelihood expectation maximization 47 .To address this limitation, a simple and fast algorithm, a point-of-closest Approach (PoCA), was developed to efficiently locate a single scattering position 34 .However, there is a key assumption in the PoCA algorithm: a muon's trajectory changes once in the medium.To find the PoCA point, we need to measure incoming and outgoing muon trajectories.The incoming muon trajectory, or initial muon trajectory before interacting with target objects, is reconstructed using two muon trackers, as shown in Fig. 2. The outgoing muon trajectory, the final muon trajectory after interacting with target objects, is also reconstructed in the same manner.Both incoming and outgoing lines are extrapolated until the shortest distance between two lines in 3D is achieved, and the middle point of a distance between two lines is a PoCA point in 3D.Then, a PoCA voxel is assigned, which includes the PoCA point in a voxelated volume of interest, as shown in Fig. 3.
From the measurement in upper and lower muon trackers, the 3D coordinates of a muon are where (x, y, z) i is the 3D muon coordinate at the ith muon tracker.The extrapolated lines, l 1 and l 2 , using two muon coordinates at the upper and lower muon trackers are where e 1 and e 2 are the unit directional vectors for incoming and outgoing muon trajectories.K 1 and k 2 are the coefficient vectors to extend lines toward the region of interest.Then, the 3D PoCA point is given by This process is repeated for every recorded muon event, and PoCA voxels can contain multiple PoCA points.

Density mapping method
To complete images of objects using MST, the assigned PoCA voxels must be visually distinguished using colors or a density mapping.The assigned PoCA voxels are colored based on the muon scattering angle, which is an angular difference between l 1 and l 2 from Eqs. ( 12) and ( 13).The scattering angle, θ s , is given by where Δθ x and Δθ y are the angle difference between l 1 and l 2 in terms of x-or y-axis.If there exist multiple PoCA points with different scattering angles in a PoCA voxel, the averaged scattering angle is used.
where N v is the total number of recorded muon events in the PoCA voxel.

Modeling Cosmic ray muon simulations
To provide statistically reliable CRM samples for muon tomography and monitoring applications, a new CRM sampler was developed based on phenomenological models of muon energy spectra and angular distribution.
The CRM sampler generates muon samples from the described energy spectrum (Eq. 1) and angle distribution (Eq.3) with user-defined parameters-number of sample generations, maximum and minimum energy limits, and zenith angles, as well as their uncertainties.Because Eqs. ( 1) and ( 3) are not standard probability distributions, statistical algorithms are employed for generating random samples.The inverse transform method is used to obtain random samples based on probability profiles from both probability distribution functions, but the insignificant probabilities (< 5 × 10 −6 ) are dropped.The generated CRM samples are saved in a table format so that it can be exploited in GEANT4.The examples of muon energy spectrum and zenith angle distribution from the CRM samplers are shown in Fig. 4.

Spent nuclear fuel dry cask storage
SNFs are stored in various designs of storage depending on fuel type, capacity, usage, and manufacturer 51 .In general, the SNF dry cask consists of two major components: a stainless-steel canister and concrete overpack.
The outermost surface of a dry cask is made of thick concrete (800-1000 mm) to shield radiation from SNF decay and is designed to maintain the total radiation dose equivalent rate at the site less than 0.25 mSv per year at the controlled boundary of the system 52 .In this work, we used a commercial SNF dry cask canister model that stores up to 24 PWR FAs.Each FA includes 15 × 15 UO 2 fuel rods (ρ UO2 = 10.97 g/cm 3 ), which have a radius, pitch, and length of 5.35, 14.3, and 3658 mm, respectively.The overall dimension of each FA is 215 mm × 215 mm × 3,658 mm.The array of FAs is surrounded by an annular concrete shielding (ρ concrete = 2.3 g/ cm 3 ) with an inner and outer radius of 863.5 and 1685 mm.In simulations, we modeled 23½, 23, and 22 PWR FAs to represent one-half, one, or two middle FAs that are missing, as shown in Fig. 5. (15)

Momentum-informed muon scattering tomography
After CRM samples are stochastically generated from the muon sampler, muons interact with various structures, including muon trackers, spectrometry, and SNF DSC.The output data from GEANT4 simulations provides posterior muon energy and positions from the installed Cherenkov muon spectrometry and muon trackers, respectively.Incoming and outgoing muon trajectories were reconstructed from output data, and the PoCA points were determined from them.Then, PoCA voxel and scattering angle computed using Eqs.( 14) and ( 15) were saved in a vector.The 3D coordinates of PoCA voxel, P x,y, and z , and scattering angle, θ, for each muon event, or muon event vector, is written by  where N μ is a total number of recorded muon events and the subscript MST represents that Eq. ( 14) is used for the MST algorithm.In a typical MST algorithm, tomographic images of target objects are reconstructed based on Eq. ( 17).In this work, a Cherenkov muon spectrometry was incorporated into the typical MST imaging system to include muon momentum information.Therefore, Eq. ( 17) can be augmented by adding an extra element, p, which is the recorded muon momentum in GeV/c.
To simultaneously consider two parameters, θ and p in Eq. ( 18), we developed a new algorithm encoding momentum and scattering angle behaviors into a single variable (M-value).An expression to compute M-value from θ and p is mathematically derived based on the exponentially inverse-proportional relationship between most probable muon scattering angle and momentum (details can be found in 53 ).
where k and M are slope and y-intercept.Although k slightly varies from − 2.56 for heavy nuclei to − 2.23 for light nuclei, its variance is negligible and was considered as a constant, k = − 2.4, for the SNF DSC simulations.On the other hand, M explicitly depends on types and sizes of materials, and it is a useful parameter to represent material properties.To calculate an M or M-value for each muon event, Eq. ( 19) was modified using θ, instead of using mod(θ).Thus, a new function, M (p, θ), is given by when k is a constant, − 2.4, and Eq. ( 20) becomes In the MMST technique, M-value is facilitated to replace the scattering angle.Therefore, Eq. ( 18) can be updated and simplified by substituting θ and p with M, Therefore, PoCA voxels are mapped based on M density to reconstruct tomographic images.

Monitoring spent nuclear fuel dry storage cask
Two NDE technical approaches were used to monitor SNF in the DSC and its structure, (1) a signal analysis and (2) an image analysis.In the signal analysis approach, scattering angle and M-value distributions for the FAs and concrete overpacks were compared when 0, 1/2, 1, and 2 FAs are missing.The signal distributions were computed at the horizontal line, which crosses the surrounding air, concrete, airgap, and FAs as shown in Fig. 7.The signal amplitude distributions of scattering angles and M-values were averaged over 0-215 mm (size of one FA) in terms of y-axis, and the results are plotted along the x-axis from − 2000 to 2000 mm.Although both scattering angle and M-value signal amplitudes are theoretically proportional to the material density, scattering angle signal amplitudes of the airgap (38-50 mrad) were similar or greater than concrete overpacks in the MST signal analysis because of the strong interference signals around the airgap such as FAs (75-83 mrad).This ( 18)    detect missing FAs.The signal stability, or reduced signal fluctuation, can be achieved because the uncertainty in scattering angles is significantly reduced by including a muon momentum term as described in Eq. ( 21).Similarly, a difference of signal amplitudes for different materials becomes distinctive by using M-values instead of scattering angles.In particular, the capability of capturing the airgap signals between FA bundles in the center and concrete overpacks (highlighted in blue in Fig. 7) is noteworthy because it implies that the MMST method can find the missing FAs as well as possible cavities in the DSC structure.Therefore, the overall capability to identify the number and exact location of missing FAs in the DSC is significantly improved by using M-value signal distribution profile.The averaged scattering angle and M-value signal amplitudes are changed from 75 to 71 mrad, and − 0.01 to − 0.29 where two missing FAs exist.The signal amplitude reduction ratios were 5.0% and 13%, respectively.Similarly, the averaged scattering angle and M-value signal amplitudes decreased from 80 to 79 mrad, and 0.09 to − 0.03 where half a missing FA exists.The signal amplitude reduction ratios are 1.0% and 5.0%, respectively.The reconstructed muon tomographic images of the DSC with 0, 1/2, 1, and 2 missing FAs are shown in Fig. 8.In each case, four cross-sectional images are presented, which were reconstructed using the MST and MMST imaging algorithms with 10 5 and 10 6 muon samples, respectively.The overall image resolution is improved by increasing the number of muon samples (see and compare left and right columns in Fig. 8) and replacing scattering angles with M-values (see and compare upper and lower rows in Fig. 8).The best resolution can be achieved by using the MMST imaging algorithm with 10 6 muon samples: one-half, one, and two missing FAs were successfully visually identified and located as shown in Fig. 8.
The benefits of using MMST algorithm-which offers distinctive signal amplitudes for different materials and an improved capability to identify and locate missing FAs-are summarized in Fig. 9.The ranges of scattering angle and M-value signal amplitudes for various structural materials in the DSC are presented in different colors in Fig. 9. M-value signals for all materials are clearly separated in MMST signal analysis, whereas scattering angle signals for the airgap and concrete are overlapped, which means two materials are indistinguishable in MST signal analysis.The averaged signal amplitudes of scattering angles and M-values for 1/2, 1, and 2 missing FAs' are plotted with 1σ error bars.The scattering angle signals are overlapped with the fully loaded FAs signal area when one and two FAs are missing with less than 60% and 30% CLs.On the other hand, the M-value signals do not overlap with the fully loaded FAs signal area when one FA is missing (> 99% CL), and one-half FA can be found with an 80% CL.

Conclusion
This work shows that the measurements of cosmic-ray muon energy using MST can improve the image resolution and enhance our ability to find and locate missing FAs in DSCs.To mathematically encode muon energy information to the existing MST imaging algorithm, a new algorithm, MMST, was developed.To demonstrate the benefits of MMST in the application for monitoring SNF, a commercial SNF dry cask canister that stores up to 24 PWR FAs was designed in a GEANT4 environment.The array of FAs is surrounded by an annular concrete shielding.Additionally, muon trackers and a Cherenkov spectrometry were also modeled using GEANT4 to simulate muon interactions with DSC materials and measure muon scattering angles and energy as shown in Fig. 6.To test the functionality of the MMST algorithm, four scenarios-with fully loaded, two, one, and onehalf FAs missing-were considered and analyzed using two approaches, signal and image analyses.enables to find the position and quantity of missing FAs in DSC investigations.The rates of correctly identifying missing FAs were improved from 79 to 98% and 51% to 88% by replacing scattering angles with M-values when one FA and one-half FA is missing, respectively.These improvements also can be visually observed in the reconstructed images (Fig. 9).The advancement in the MST technique for DSC verification is expected to resolve long-standing problems in international nuclear safeguards and nuclear material accountancy 54 .With a traditional MST technique, approximately 90 days of measurement time is required to survey the entire DSC 17 .By utilizing muon energy information in the MST imaging algorithm, we found that the required measurement time is significantly reduced, and comparable results were acquired by the smaller number of muon samples to survey missing FAs in the DSC using MMST instead of MST.Given the current absence of fully licensed NDE techniques for monitoring SNF in DSCs, the MMST imaging algorithm emerges as a highly promising solution.Nevertheless, comprehensive field testing is required for it to attain licensure from authorized agencies like the Nuclear Regulatory Commission (NRC).Future investigations will focus on a sensitivity study of MMST, aimed at detecting missing fuel pins and monitoring material failures such as cracks, imperfections, and swelling in structural materials and fuel pins.Enhancing the resolution of MMST images necessitates: (1) improved muon momentum and position measurement capabilities, (2) enhanced particle tracking algorithms, and (3) increased muon flux.Enhancements in radiation detection instrumentation are expected to reduce statistical measurement uncertainties.Moreover, recent research demonstrates the feasibility of estimating the most probable muon trajectory within the target object using Bayesian theory-based algorithms 55 .While increased muon flux can be achieved through adaptation of a muon beam accelerator facility 56,57 deployment in field conditions remains impractical.

,Figure 1 .
Figure 1.Differential intensity of cosmic ray muons in the range of 0.1-10 GeV/c and the variation of muon flux with a zenith angle range of 0°-89° at sea level 54 .

Figure 2 .
Figure 2. Reconstruction of incoming and outgoing muon trajectories using two-fold upper and lower muon trackers.Four 8 by 8 muon trackers are shown in the example.

Figure 3 .
Figure 3. Visualized PoCA algorithm to locate a single scattering point using incoming and outgoing muon trajectories (left) and the assigned PoCA voxel in a 3D voxelated volume of interest (right).

Figure 5 .
Figure 5. Layouts of SNF assemblies in a dry cask when (a) they are fully loaded and when (b) two, (c) one, and (d) one-half of middle FA(s) are missing.

Figure 6 .
Figure 6.Overview of instrumentations and DSC for the momentum-informed muon scattering tomography system using a Cherenkov muon spectrometry for SNF imaging and monitoring.

Figure 7 .
Figure 7. Layouts of SNF assemblies, concrete overpacks, and airgap in the DSC when FAs are fully loaded, two middle FAs are missing, one middle FA is missing, and one-half FA is missing (from left to right).Averaged scattering angle and M-value signal distributions (0 < y [mm] < 215) are plotted along the x-axis (− 2000 < x [mm] < 2000).

Figure 8 .
Figure 8. Reconstructed muon tomographic images for SNF assemblies, concrete overpacks, and airgap in the DSC when two, one, and one-half FAs are missing (from left to right).For each case, four cross-sectional images are presented that were reconstructed using the MST and MMST imaging algorithms with 10 5 and 10 6 muon samples, respectively.The overall image resolution is improved by either increasing the number of muon samples from 10 5 to 10 6 (left and right columns) or replacing the MST algorithm with the MMST algorithm (upper and lower rows).

Figure 9 .
Figure 9. Visualized signal analyses in MST (left) and MMST (right) which present scattering angle and M-value signal amplitude distributions for fully loaded FAs (red), concrete shielding (gray), and airgap (blue) when one-half, one, and two FAs are missing.Signal stabilization, enhanced capability to distinguish various materials, and find missing FAs are observed.Error bars represent 1σ.

Table 1 .
Averaged scattering angle signal amplitudes in MST for various structural materials in the DSC and various numbers of missing FAs.*Averaged signal amplitudes from 6 FAs in a row except missing FAs.**Airgap between FAs and concrete overpacks which is highlighted in blue in Fig.7.

Table 2 .
Averaged M-value signal amplitudes in MMST for various structural materials in the DSC and various numbers of missing FAs.