Correlation between the Storage Pattern of Spent Fuel Casks and the Neutron Multiplication Factor

Mosebetsi J. Leotlela ^1,2,* 

1. School of Physics, University of the Witwatersrand, 1 Jan Smuts Avenue, Braamfontein, Johannesburg 2017, South Africa

2. Neutronix Nuclear Consulting Service (Pty.) Ltd., 6 Chemico Street, Crystal Park Benoni 1515, South Africa

* Correspondence: Mosebetsi J. Leotlela 

Academic Editor: Fabian Ifeanyichukwu Ezema

Received: April 21, 2025 | Accepted: November 17, 2025 | Published: December 19, 2025

Journal of Energy and Power Technology 2025, Volume 7, Issue 4, doi:10.21926/jept.2504018

Recommended citation: Leotlela MJ. Correlation between the Storage Pattern of Spent Fuel Casks and the Neutron Multiplication Factor. Journal of Energy and Power Technology 2025; 7(4): 018; doi:10.21926/jept.2504018.

© 2025 by the authors. This is an open access article distributed under the conditions of the Creative Commons by Attribution License, which permits unrestricted use, distribution, and reproduction in any medium or format, provided the original work is correctly cited.

1. Introduction

When containers of fissile material—such as spent fuel casks in the spent fuel storage facility—are placed close together, the neutron interaction among them significantly influences the neutron multiplication factor, k_eff, of the system quite considerably compared to the k_eff of a single fissile system in isolation. This is due to the interaction of the neutron flux from adjacent casks.

The k_eff of the entire spent fuel facility at any given point in time and space is thus a function of the spatial arrangement of its subsystems and is further exacerbated by the storage array/matrix adopted in storing the casks.

Consequently, time-dependence (see Ref. [1]) and spatial-dependence of the k_eff should be assessed as part of the safety analyses of spent fuel storage facilities—the spacetime dependence of both the neutron multiplication factor and the dose rate is a crucial component of safety assessments and should be considered a standard aspect of nuclear safety analysis. This assessment should precede the placement of spent fuel casks in storage facilities with a view to selecting a storage matrix that will optimise storage space by reducing the amount of space used, while keeping the k_eff below the regulatory limit of 0.95 [1,2].

Failure to analyse the impact of the spatial dependence of the k_eff on a nuclear subsystem may result in exceeding the regulatory limit of the k_eff emanating from the spent fuel storage facility and consequently lead to nuclear criticality accidents [3].

The effect of spatial arrangement on the k_eff of a nuclear system was initially published in 1955, using fissile fuel pellets [4,5,6,7]. It was discovered that, for specific storage matrices, the k_eff increases, whereas in others it may decrease. These studies were later followed by another investigation, which was conducted in 2012 using a simulation of spent fuel casks, and the results were published in a PhD thesis (see Refs. [3,8]). In the present study, a discrete random variable sampling technique was used in which the spent fuel casks were separated at regular intervals defined by a mathematical equation (to be discussed in the methodology section). The effect of changing the storage matrix and the separation gap (d) between spent fuel casks—measured by k_eff—was analysed to determine the relationship $(x,f(x))$ between the separation gap and k_eff [8].

The novelty of this work lies in establishing a direct link between the solid angle Ω and the k_eff of the storage matrix. Any change in the dimensions of the storage matrix results in a change in the magnitude of Ω, which, based on the correlation equation between Ω and k_eff, allows the calculation of the new k_eff at a new distance between casks. In addition, this work also derives a correlation equation between the distance between casks (d), or between rows of casks (D) in a storage matrix, and the k_eff of the system.

The number of casks and storage matrices simulated was arbitrarily based on the projected number of casks that would be required to sustain the existing operating nuclear power reactors commissioned in the 1980s, taking into consideration their size, the reload frequency, and the minimum number of fuel assemblies unloaded. This was based on the assumption that they would be relicensed for an extended operational licence.

2. Methodology

The simulation was performed using STARBUC, a burnup module of the Standardized Computer Analyses for Licensing Evaluation (SCALE) code. SCALE computer code is a suite of nuclear safety analysis tools developed by Oak Ridge National Laboratory (ORNL). It is a modular code consisting of several functional modules—some based on Monte Carlo techniques (e.g., KENO), while others are based on deterministic methods—all coupled together to allow the user to perform complex calculations.

Because of its robustness, it is highly trusted and widely used in the nuclear industry worldwide. It is used by both nuclear regulators and nuclear reactor operators alike to perform nuclear criticality safety analysis, reactor physics, and radiation shielding.

Because several factors affect the precision of Monte Carlo results, some of which include the number of histories simulated, the sample size of the particles simulated, the tally type, and the variance reduction techniques applied, a large number of neutron generations—5000 neutron generations and 1500 neutrons per generation—were used.

In addition, given that the number of neutron generations are directly related to the precision and the variance of the results, if the number of neutron generations is low it will result in large variance values, taking into consideration that the estimated relative error squared (R²) is proportional to $\frac{1}{n} $ where n is the total number of histories that are calculated in the simulation. Furthermore, the low neutron generation number (which may be referred to colloquially as a sample size) is related to the number of histories each neutron would be subjected to and the computation time, T, required to complete each simulation, also known as a run. Thus, the relationship between T and R² is that R²T should be constant. Therefore, because there is a relationship between R² and T, it is possible to reduce the variance associated with the short computation time in this simulation. This would increase the number of histories quite substantially compared to if a low number of neutron generations, say 1000, were applied. The sigma was set at 0.0005 so that the cut-off time could be reached either when all 5000 neutron generations had been run or when a sigma of 0.0005 was reached, whichever occurred first.

The number of generations skipped (nsk) before testing for convergence was put at 500, which is high enough to stabilise any fluctuation before convergence [3,8].

The analysis was based solely on spent fuel exposed to burnup in Table 1. Furthermore, since the primary objective was to determine what the effect of varying storage matrices/patterns would be on the k_eff, the only nuclide set that was taken into account for burnup credit analysis was the actinide nuclide set. The project was divided into two main categories: thirty casks and four casks. The casks were further subdivided into groups: four casks which were split into a 2 × 2 array and a 1 × 4 array and the thirty casks in a 2 × 15 and a 3 × 10 array. The choice of these arrays was based on the arbitrary dimensions of the spent fuel cask storage building.

Table 1 Burnup regimes (MWD/MTU).

In addition to the division of the project into these geometric categories, each category was run in three burnups designated A, B, and C (MWD/MTU) (as reflected in Table 1), to assess whether the combination of the burnup and the storage pattern would influence the k_eff of the system.

The enrichment (ϵ) in all cases was 3.25 wt.%, of which part of the input deck is shown. For a complete input file, which shows all simulation parameters, the reader is referred to Appendix 1.

Except for the Zigzag storage matrix, the dimensions of which differ from the others, all casks under study were modelled as stored in a cask storage building—refer to Figure 1—with arbitrary dimensions as indicated in the input deck below:

com = cask storage building.

Figure 1 Cask storage building.

To determine the way in which fission and absorption are affected by the change in geometry of the storage matrix, i.e., the spatial location of the cask, fission and absorption were taken from the output file and plotted alongside their respective distance from the adjacent casks. These, i.e., fission and absorption, were taken from node 9 of unit 1 of the output file, bearing in mind that unit 1 was divided into 18 axial profiles and node 9 was therefore in the middle of unit 1 and thus had the highest nuclear reaction rates.

The analysis was divided into two categories—four casks and thirty casks—which will be described in detail in their respective section.

2.1 Four Casks

To determine the influence of storage pattern/matrix on the k_eff, for four casks, the four casks were stored in two different matrices: Linear storage matrix (i.e., a 1 × 4 storage array) as depicted in Figure 2 [8].

Figure 2 1 × 4 cask storage matrix (25 cm apart).

2.1.1 Linear Storage Matrix

To meet the objective of determining the impact of separation distance on the k_eff in the linear storage matrix (1 × 4), the distance, d, between adjacent casks was only varied along the x-axis. The separation distance was increased by a constant distance of 25 cm, starting from d = 0 (i.e., casks in contact), taking into consideration the cask diameter (Ф) of 288.6 cm. The x-coordinate of the cask at the origin was fixed for all separation gaps at x = 0, y = 0, and z = -180, while the x-coordinates for the second cask (column 2) were calculated using the following equation: (Eqn. 1)

\[ x_d=|d|+\Phi \tag{Eqn.1} \]

where

X_d = the x-coordinate at distance d from the adjacent cask.

d = the distance between adjacent casks.

Ф = 288.6 cm, which is the diameter of the cask.

x_d for the 3^rd cask (in column 3) was obtained by -|x_d|, while x_d for the 4^th cask (in column 4) was obtained by 2(x_d) indicated in column 3.

The x-coordinate at

d = 25 (d₂₅) for column 2 is Ф + 25= 313.6 as in Eqn. 1.

d₅₀ = 313.6 + 50 = 363.6.

d₇₅ = 363.6 + 75 = 438.6.

⁞

d_n = x-coordinate before + |d|, i.e. d_n = x_(n-25) + |d|.

where x_(n-25) is the x-coordinate at distance (d_n-25).

Thus, the sequence for x-coordinate at d_n in column 2 may be summarised as

\[ _{\mathrm{column}2}x_{(d_n)}=x_{=|n-25|}+|d_n| \tag{Eqn.2} \]

The x-coordinates in columns 3 and 4 were represented by Eqn. 3 and Eqn. 4 respectively.

\[ _{\mathrm{column}3}x_{(d_n)}=-_{\mathrm{column}2}x_{(d_n)} \tag{Eqn.3} \]

\[ _{\mathrm{column}4}x_{(d_n)}=2\times\left ( {_{\mathrm{column}2}x_{(d_n)}}\right ) \tag{Eqn.4} \]

Thus, for the following the sequences from Eqn. 1 to Eqn. 4 the coordinates of the linear cask storage array (refer to Figure 2) are summarised in Table 2.

Table 2 Spatial positions of four casks in a 1 × 4 storage matrix in the cask storage building.

2.1.2 Square Storage Matrix

With regard to the square storage matrix (i.e. 2 × 2 matrix), the distances between adjacent casks were changed along both axes, i.e. along the x and the y-axis, with the initial separation gap between adjacent casks set at 0 (zero) where all casks were in contact and with |d| increasing by a regular interval of 25 cm for every space until the maximum distance was reached. Following the same technique, the distance between each cask and the (x, y, z) coordinates of a 2 × 2 cask (refer to Figure 3).

Figure 3 2 × 2 Cask storage matrix.

To obtain the (x, y, z) coordinates of the spatial position of each cask in the storage building, |d| between adjacent casks was added to Ф.

Thus, obtaining the x-coordinate for the 2 × 2 matrix for |d|₂₅ = 25 cm was performed as follows: Ф + |d|₂₅ = 313.6, and y = |x|, -|x|, z = -180.2. Similarly, the x-coordinates for d₅₀ (at 50 cm) were Ф + |d|₅₀ = 338.6 and y = |x|, -|x|, z = -180.2 etc. The respective |d| and (x, y, z) coordinates of the 2 × 2 cask storage matrix are summarised in Table 3.

Table 3 Spatial position of the four casks in the 2 × 2 storage matrix in the cask storage building.

2.2 Thirty Casks

To determine the effect on the k_eff of the storage matrix with 30 casks, the same technique was applied as to the four casks: with the casks being arranged in 2 × 15 and 3 × 10 arrays. Each one of them was subjected to the burnup regimes A, B, and C listed in Table 1.

2.2.1 3 × 10 Storage Matrix

The 3 × 10 array was subdivided into three subgroups according to the height of elevation of the middle row of casks 10 cm (Z_h = -170) and 20 cm (Z_h = -160). The 0 cm (Z_h = -180.2) elevation is regarded as the base position as the casks rest on a slab (i.e., elevation = 0). Contrary to the d of four casks, in the 30 casks, d referred to the distance between the rows of casks, instead of distance between individual casks as was the case with the four casks, the 3 × 10 vertically misaligned casks are studied in section 3.2.1.

2.2.2 2 × 15 Storage Matrix

About the 2 × 15 cask storage facility, the two rows of casks were separated by distance D^a where D was added to the cask diameter (i.e. D + Ф) which was shifted along the y-axis. The 2 × 15 storage matrix was divided into two categories: the linearly aligned storage matrix (i.e., the traditional storage matrix) and the horizontally misaligned (i.e., 2 × 15 Zigzag storage matrix).

The 2 × 15 Traditional Linear Storage Matrix. In the 2 × 15 traditional storage matrix, the casks were separated by a regular interval of 50 cm, starting with an initial gap of 100 cm and a maximum distance of 500 cm. Thus, using the same technique, for D = 100 cm, the x, y, z coordinates defining the spatial locations of the casks were obtained by D + Ф, where Ф is the diameter of the cask (refer to Figure 4). Subsequently, y = 100 cm + 288.6 cm, and z was kept constant at z = -180.

Figure 4 2 × 15 traditional storage matrix 100 cm apart.

Thus, the coordinates of the reference cask where D = 100 cm are (Cask #1) (x, y, z) = (0, 388.6; -180.2). From there on, the y-coordinate of cask # n is obtained by adding 388.6 to the y-coordinate of the prior cask to get the y-coordinate of the next cask. Therefore,

Cask #2: 388.6 + 388.6 = 777.2

Cask #3: 777.2 + 388.3 = 1165.8

Cask #4: 1165.8 + 388.6 = 1554.4

⁞

Thus, based on this sequence, the y-coordinate of cask n(y_n) is obtained by Eqn. 5, which is summarized in Table 4 and incorporated in input files for the 2 × 15 storage array (refer to Appendix 1).

\[ y_n=y_{n-1}+388.6 \tag{Eqn.5} \]

where

\[ y_n=y_{co-ordinates\text{ of cask n}} \]

\[ y_{n-1}=y_{co-ordinates\text{ of cask n-1}} \]

Table 4 Coordinates of casks in a 2 × 15 storage matrix of a small sample of distances along the y-axis (D_y).

2 × 15 Zigzag Storage Matrix. In addition to the 2 × 15 traditional storage array indicated in Figure 4, the 2 × 15 storage array was rearranged to form a Zigzag pattern to determine what the effect of the horizontal misalignment of casks would be in the k_eff of the system (refer to Figure 5).

Figure 5 2 × 15 Zigzag storage matrix.

3 × 10 Storage Matrix (Vertically Misaligned Row). In the case of the Vertically Misaligned 3 × 10 storage matrix, the separation was achieved by shifting the two outer rows outwards and keeping the middle row in the same place at y = 0 but raising its z-axis to z = -160.2, to achieve vertical misalignment. The objective of this was to determine whether any misalignment in the axial profile would influence the k_eff.

To find the x-coordinates of the 3 × 10 storage array for a given D, Eqn. 5 was modified to $x_{n-1}+388.6$.

where

\[ x_n=x_{co-ordinates\text{ of cask n}} \]

\[ x_{n-1}=x_{co-ordinates\text{ of cask n-1}} \]

The coordinates of the 3 × 10 storage array for D = 100, 150 and 200 and the middle row is elevated by 20 cm and are listed in Table 5 (and depicted in Figure 6). They are applied in the respective source files.

Table 5 Coordinates of the vertically misaligned (20 cm elevation, i.e. z = -160.2) 3 × 10 storage matrix where the y-coordinates of the outer rows are separated by D_y = 100, 150, and 200.

Figure 6 The 3 × 10 storage matrix, 100 cm apart and elevated by 20 cm (i.e. Z_h = -160.2 cm).

3. Results

The results are presented in two categories in line with the simulations comprising 4 and 30 casks and are further divided into subgroups according to the storage matrix.

In all cases, a polynomial and a linear correlation between k_eff and d (D in the case of 30 casks) were plotted and the corresponding linear and polynomial relationships were derived. Because the same data were used to plot linear and polynomial correlations, the linear correlation coefficient (R) was low (in many cases below 0.5), while R² for the polynomial correlation was >0.8, which is confirmation that one cannot have a perfect fit for both linear and polynomial correlations with the same data.

3.1 Four Casks

The two storage matrices were compared based on the k_eff as a function of d and the corresponding linear regression (R) derived from the correlation equation to determine whether the k_eff increases or decreases.

3.1.1 Linear (1 × 4) Versus Square (2 × 2) Storage Matrix

The results were grouped into three categories in accordance with the burnup regime simulation that was performed (burnup A, B, or C MWD/MTU).

Low Burnup (Burnup Regime A). The results in Figure 7a indicate that for burnup A, the k_eff of the linear storage array decreases quite rapidly (refer to Table 6) for the slope B and (correlation coefficient R) as the distance between adjacent casks increases. In contrast, the k_eff of the 2 × 2 matrix increases steadily with B = 2.07E-08 (refer to Table 6).

Figure 7 k_eff of 4 Casks as a function of the storage matrices in three different burnup regimes.

Table 6 Linear fit values for parameters of the 1 × 4 vs the 2 × 2 storage matrix at burnup A.

Midrange Burnup (Burnup Regime B). When burnup is increased from burnup regime A to burnup regime B (refer to Figure 7b and Table 7), it is observed that the k_eff of both the 1 × 4 and 2 × 2 storage matrices increase; However based on their respective B and R values, the 1 × 4 matrix increases more rapidly (refer to Figure 7b and Table 7).

Table 7 Linear fit values for parameters of the 1 × 4 vs the 2 × 2 storage matrix at burnup B.

High Burnup (Burnup Regime C). The results performed at higher burnup (burnup C) indicate that the k_eff for the 1 × 4 storage matrix increases at a rapid rate with slope, B = 1.72E-06, and R = 0.20942. In contrast, the k_eff for the 2 × 2 storage matrix decreases with B = -5.65E-07 and R = -0.14125 (refer to Figure 7c and Table 8).

Table 8 Linear fit values for parameters of the 1 × 4 vs the 2 × 2 storage matrix at burnup C.

The equation for the linear correlation 1 × 4 and the 2 × 2 arrays is as follows:

\[ y=A+Bx \tag{Eqn.6} \]

The rapid decrease in k_eff for the 1 × 4 storage matrix (refer to Figure 7a) is a testament to the neutron leakage associated with the storage array. It is also consistent with the analysis reached when using Eqn. 10 and Eqn. 11—to be discussed later. This shows that for 1 × 4 array, as D increases, there will be a situation where $B_g^2>B_m^2$ which enables neutron leakage and because of that k_eff will decrease.

In a 2 × 2 square matrix, as observed from Figure 7a, the k_eff increases gradually with an increase in D—refer to the slope of the linear regression line in Table 6. This increase is because the storage array design minimizes neutron leakage. Also, the void formed between the four casks creates a flux trap--which, because of scattering, increases the fission density and consequently increases k_eff. This analysis is consistent with the solid angle method, as the representative unit will be in the middle of the void; therefore, an increase in D will have an insignificant effect on H. As a result, the increase in k_eff is driven primarily by either $B_g^2<B_m^2$ or $B_g^2=B_m^2$.

When the burnup increases between B and C, there is an increase in fission products, which absorb neutrons and decrease the k_eff.

In deriving the mathematical relationship between k_eff and d, given that we are dealing with an equation in three dimensions, the results suggest that the relationship between k_eff and d for the 1 × 4 storage matrix at burnup A is given by Eqn. 7, and the variance for various parameters is summarised in Table 9.

\[ \begin{aligned} y_{1\times4(20\,MWD/MTU)}=k_{eff} &=0.88422+9.21246\times10^{-5}d-1.32537\times10^{-5}d^2 \\&+47913\times10^{-7}d^3-1.58339\times10^{-8}d^4+2.13535\times10^{-10}d^5 \\&-1.65512\times10^{-12}d^6+7.32959\times10^{-15}d^7 \\&-1.72149\times10^{-17}d^8+1.66199\times10^{-20}d^9 \end{aligned} \tag{Eqn.7} \]

Table 9 Variance for parameters of the 1 × 4 storage matrix (burnup A (MWD/MTU)).

For the 1 × 4 storage matrix at burnup B, the results suggest that the relationship between k_eff and d for the 1 × 4 storage matrix is given by Eqn. 8 and the variance of the parameters is listed in Table 10.

\[ \begin{aligned} y_{1\times4\,cask\,storage\,matrix(30\,MWD/MTU)} &=0.86337-5.0402\times10^{-4}d+5.14274\times10^{-5}d^2 \\&-2.01435\times10^{-6}d^3+4.0448\times10^{-8}d^4 \\&-4.6477\times10^{-10}d^5+3.17736\times10^{-12}d^6 \\&-1.27647\times10^{-14}d^7+2.77936\times10^{-17}d^8 \\&-2.5287\times10^{-20}d^9 \end{aligned} \tag{Eqn.8} \]

Table 10 Variance parameters for 1 × 4 storage matrix (burnup B MWD/MTU).

Finally, for the same storage matrix, i.e., 1 × 4, at burnup C, the relationship between k_eff and d is given by Eqn. 9 with the parameters displayed in Table 11.

\[ \begin{aligned} y_{(1\times4;(40\,MWD/MTU))} &=0.8539+3.57277\times10^{-4}d-3.27356\times10^{-5}d^{2} \\&+1.00406\times10^{-6}d^3-1.47871\times10^{-8}d^4 \\&+1.18375\times10^{-10}d^5-5.28977\times10^{-13}d^6 \\&+1.24677\times10^{-15}d^7-1.23556\times10^{-18}d^8 \\&+8.66878\times10^{-23}d^9 \end{aligned} \tag{Eqn.9} \]

Table 11 Variance for the parameters of the 1 × 4 storage matrix (burnup C MWD/MTU).

3.1.2 Factors Influencing the k_eff Trend of the Storage Four Casks

The investigation analysed four casks in linear (1 × 4) and in square (2 × 2) storage matrices, as shown in Figure 2 and Figure 3, respectively, under three different burnups as indicated in Table 1. Upon analysing the trends of the effect storage pattern on k_eff by various burnup regimes as represented by Figure 7, it is noted that the factors that determine the trend are [2]:

The Material and Fission Density. In reactor physics, nuclides are typically grouped into three nuclide sets based on their fractional contribution to neutron absorption. These sets are further classified into three categories, each characterised by a specific range of fractional contribution to neutron absorption. For actinides, the categories are: dominant absorbers, moderate absorbers, and least absorbers. Furthermore, within each category, individual nuclides contribute differently to the total absorption. Consequently, their presence in the material composition of the fissile system—such as spent fuel casks—will result in a corresponding fractional contribution to neutron absorption, which in turn will lead to a proportional change in k_eff. The three categories are [3]:

• Actinides (²³⁴U, ²³⁵U, ²³⁸U, ²³⁸Pu, ²³⁹Pu, ²⁴⁰Pu, ²⁴¹Pu, ²⁴²Pu and ²⁴¹Am),
• Actinides + Minor Fission Products (²³⁴U, ²³⁵U, ²³⁸U,²³⁸Pu, ²³⁹Pu, ²⁴⁰Pu, ²⁴¹Pu, ²⁴²Pu, ²⁴¹Am, ²⁴³Am, ²³⁷Np,¹³³Cs,¹⁴³Nd, ¹⁵¹Sm and ¹⁵⁵Gd) and
• Actinides + Principal Fission Products (²³⁴U, ²³⁵U, ²³⁸U, ²³⁸Pu, ²³⁹Pu, ²⁴⁰Pu, ²⁴¹Pu, ²⁴²Pu, ²⁴¹Am, ²⁴³Am, ²³⁷Np, ⁹⁹Tc,¹³³Cs, ¹⁴³Nd ¹⁴⁵Nd, ¹⁴⁷Sm, ¹⁵⁰Sm, ¹⁵¹Sm, ¹⁵²Sm, ¹⁵¹Eu, ¹⁵³Eu, ¹⁵⁵Gd) with their respective nuclides in brackets [9,10,11,12].

It is well known that when the burnup of fissile material is increased, as shown in the burnup regimes listed in Table 1, the material density of some nuclides will increase while others will decrease [3,8]. Considering the actinide nuclide group, an increase in burnup results in an increase in the material density of all nuclides in this group except for ²³⁵U and ²³⁴U, which decrease. It is also important to note that increases or decreases in various nuclides occur at different rates, dictated by their respective sensitivity coefficients.

Regarding the responses of fission products (minor and principal) to increased burnup, it has been found that all increase with burnup. It is therefore expected that when spent fuel is subjected to increasing burnup regimes A to C (refer to Table 1), the response will follow the same trend. Thus, the ranking of the above nuclide sets based on the total absorption cross section ($\sigma_{total}^a$) of the three nuclide sets listed above is as follows (in order of decreasing magnitudes (Principal Fission Products > Minor Fission Products > Actinides)). Thus, under a similar geometric arrangement (i.e., similar storage matrices), simply changing the burnup regime is expected that the k_eff would respond accordingly.

The significance of the influence of $\sigma_{total}^a$ is noticeable in the slope of the linear fit graph of 1 × 4 configuration depicted in Figure 7, where it is observed that for burnup A, there is a rapid decreases in k_eff (with a slope B = -2.48E-06, refer to Table 6) as the distance between casks increases. However, when the burnup increases to burnup regime B, the slope of the linear fit changes to 3.77E-07 (refer to Table 7), indicating a gradual increase in k_eff with an increase in distances between adjacent casks, with the y-intercept changing from 0.88372 in burnup A to 0.86334 in burnup B. Finally, when the fuel burnup is increased to Burnup regime C, there is a noticeable increase in k_eff with distance between casks, as shown by the slope of the linear regression line of the 1 × 4.

Furthermore, in addition to changes in slopes of linear regression curves, the k_eff generally decreases as the burnup regime changes from A to C, as has been reported previously [3].

When studying the effect of a change in Burnup on the k_eff with respect to the square matrix (i.e., 2 × 2 storage array), it is noted that a shift in burnup from A to B has absolutely no effect on k_eff. This is confirmed by the fact that the slope (i.e., the value of B in the linear regression curve) is the same, i.e., B = 2.07E-08 (refer to Table 6 and Table 7) for both burnup regimes (i.e., burnup regimes A and B). However, when one studies the effect of a change in area occupied by the casks with respect to Burnup regime C, it is noted that there is a rapid decrease in k_eff as the distance between adjacent casks increases, which is confirmed by the slope, which is visibly negative (B = -5.65E-07 as shown in Table 8).

The indifference between k_eff between burnup A and B in the same spatial arrangement is because there not much difference in ratio of fission rate-to-absorption rate ($\frac{\sigma_{f}}{\sigma_{\gamma}+\sigma_{a}}$) of the two burnups. Furthermore, because of the configuration of the casks, they form a “void” among them, which results in a flux trap where high-energy neutrons are trapped and undergo thermalisation and eventually contribute to the thermal non-leakage probability, which helps in keeping the k_eff reasonably similar. With regard to the effect of spatial arrangement on the k_eff in respect of burnup regime C, there is a decrease in k_eff as the distance between casks increases. This may be ascribed to the large number of nuclides created as a result of the rise in burnup, then $\frac{\sigma_{f}}{\sigma_{\gamma}+\sigma_{a}}$ decreases rapidly, giving rise to the results displayed in Figure 7.

Influence of the Type of Storage Matrix on k_eff. The geometry of the reactor and of the type of storage matrix influences the k_eff, because of the dimensions of the spent fuel storage facility or reactor design—such as its height, shape, and radius of the cylinder, in the case of cylindrical design—which determine whether neutron leakage will occur or not. The relationship between the geometrical configuration of the spent fuel storage facility (i.e., the storage matrix of spent fuel casks in a short or long-term storage facility), the material composition of fissile material, and the k_eff of the system is described by the geometric buckling ($B_g^2$)—named “geometric” since it is a function of the spatial arrangement of the fissile material—which is compared to the material buckling ($B_m^2$) [4,5,6,13]. The $B_m^2$ is a function of the material composition and is an indication of the neutron production and absorption rate. The relationship between these parameters is as follows: If $B_g^2>B_m^2$, the neutron leakage is too high; if $B_g^2=B_m^2$ the system is just stable, i.e., neutron production is self-sustaining; however, if $B_g^2<B_m^2$, the system is critical with a net increase in neutron production.

Because casks in spent-fuel storage facilities are arranged in arrays, the solid-angle $\Omega $ method is generally more applicable than alternative approaches [14]. The method is based on the idea that the effective multiplication factor (k_eff) of fissile material in an array depends on the k_eff of a single, most-representative unit together with the probability that a neutron emitted from that unit will escape and interact with another unit—an interaction that can raise the array’s overall k_eff. That interaction probability is a function of the solid angle subtended by the surrounding units as viewed from the central most representative unit. Based on data from numerous experiments, it was established that the correlation between $\Omega $ is given by

\[ \Omega_{allowable}=9-10k_{eff} \tag{Eqn.10} \]

where $ \mit{\Omega}_{allowable}$ represents the allowable solid angle that may be subtended at the most representative fissile unit of the array, and the k_eff is the effective multiplication factor of the unreflected fissile unit in the array. Given that the casks are cylindrical, then $ \mit{\Omega}_{allowable}$ is given by Eqn. 11, read in conjunction with Figure 8:

\[ \Omega_{allowable}=\frac{LD}{H\sqrt{(L/2)^2+H^2}} \tag{Eqn.11} \]

where

L = Length of the cylinder,

D = Diameter of the cylinder,

H = distance from point P to the surface of the cylinder.

Figure 8 Solid Angle Approximate Formula for Point to Cylinder.

Point P (refer to Figure 8) denotes the centermost unit in the array. According to this method, the total solid angle subtended at point P is obtained by summing the individual solid angle contributions from all other units in the array.

In the context of 4 casks in a linear storage configuration (i.e., 1 × 4 storage matrix depicted in Figure 2) and square matrix (Figure 3), it is observed that for burnup regime A, the k_eff of the 1 × 4 matrix decreases with an increase in distance (refer to Figure 7a)—which is consistent with the increase in H in solid angle method which result in the decrease in k_eff. On the contrary, in the case of the 2 × 2 matrix, for the same burnup regime, the k_eff is increasing, albeit gradually. This also indicates that for a 1 × 4 storage matrix, there is a significant neutron leakage, which implies that $B_g^2>B_m^2$. In contrast, in the case of a 2 × 2 storage matrix, the k_eff increases gradually, which may be ascribed to the fact that;

• In a 2 × 2 array, a void is formed in the middle of the four casks due to the storage pattern, which acts as a flux trap, scattering neutrons until they reach thermal energy. Because of this, there is an increase in fission density, which leads to a rise in k_eff.
• In a 2 × 2 array, the centermost unit, positioned at point P, is in the centre of the square formed by the four casks at each corner of the square—i.e., P is inside the void (refer to Figure 3). So, for a 2 × 2 array, all four casks are equidistant from P—implying that H is of equal length for all casks. Because of that, H has no effect on the $\Omega_{allowable}$, which explains the gradual increase in k_eff. This response is different from the 1 × 4 matrix, where the four casks have significantly different magnitudes of H, which makes its contribution in $\Omega_{allowable}$ observable.

Furthermore, because of the abovementioned factors, it may be inferred that $B_g^2<B_m^2$.

When studying the response of k_eff on the same storage patterns (i.e. linear and square storage pattern) with respect of burnup regime B and C, the results (refer to Figure 7b and Figure 7c respective), it is observed that in the case of burnup regime B, the k_eff of both 1 × 4 and 2 × 2 storage matrices increases. However, the k_eff of 1 × 4 increases much more rapidly than that of the 2 × 2 configuration, as indicated by their respective slopes of 3.77E-07 and 2.07E-08 (refer to Table 7). The rapid increase in k_eff of the 1 × 4 matrix can be attributed to the fact that there are two casks on the extreme ends (i.e., Figure 2) that have the largest H, which, according to the solid angle method, is inversely proportional to $\Omega_{allowable}$. As such increasing H results in a decrease in $\Omega_{allowable}$ which when applied to Eqn. 10, will result in an increase in k_eff.

In the case of burnup regime C, it is observed from Table 8 that, based on the slopes of 1.72E-06 and 5.65E-07 for 1 × 4 and 2 × 2, respectively, and Figure 7c, the k_eff of the 1 × 4 matrix increases with distance, whereas that of the 2 × 2 matrix decreases. This may be attributed to two factors: in a 1 × 4 matrix, the value of H is significantly different between casks, and increases as d increases. With a 2 × 2 matrix, on the other hand, all four casks are equidistant from P, making the value of H the same for all of them; even if d is increased, the magnitude of H for all casks is the same. Therefore, the value of H plays a major role in 1 × 4, because as H increases, the solid angle also increases and consequently the k_eff. This is contrary to what a 2 × 2 array experiences; all 4 casks have the same value of H, making an insignificant contribution to the solid angle. Furthermore, since there is an increase in burnup from A to B and C, the accumulation of fission products occurs due to an increase in burnup (refer to (section 3.1.2)). Fission products act as neutron poisons and reduce reactivity, thus decreasing the k_eff. Furthermore, the central void between the four casks also plays a significant role in the decrease in k_eff since neutrons that previously were responsible for high fission density—and subsequently in an increase in k_eff on the lower burnup—have now been absorbed by the accumulation of fission products, which explains the decrease in k_eff.

3.2 Thirty Casks

In addition to the spatial analysis of four casks performed above, it was prudent that an analysis of more casks be performed because nuclear power plants that were commissioned in the 1980s are running out of spent fuel storage space; as such, more casks will be required. To that effect, a spent fuel storage facility must be built to accommodate the casks from the plant. Also, because many of these plants (i.e., nuclear power plants commissioned in the 1980s) are already applying to their respective nuclear regulators for licenses to operate beyond their existing operating licences.

Accordingly, this section focuses on the analysis of the relationship between k_eff and D in 30 casks for a number of storage matrices subjected to the same burnup as was used in the four-cask analysis, i.e., burnup regimes A, B, and C (refer to Table 1, actual power of various regimes). The 30 casks were divided into five categories as follows:

i. 2 × 15 traditional storage matrix: In this case, the 30 casks are divided into two rows, each of 15 casks, and rows are separated by distance D, which is increased by a constant D_y = 50 cm for every run (refer to Figure 4 and Table 4).

ii. 3 × 10, Z_h = -160 cm storage matrix: In this storage pattern, the 30 casks are arranged in three rows, each of ten casks. The middle row is lifted by 20 cm from Z_h = -180 to Z_h = -160 cm (refer to Figure 6). It is worth noting that Z_h = -180 is the datum line where the cask is resting on the concrete slab (i.e. Z_h = 0 cm).

iii.3 × 10, Z_h = -170 cm storage matrix: The same arrangement as in the 3 × 10, Z_h = -160 cm storage matrix, except that Z_h = -170 cm.

iv.3 × 10, Z_h = -180 cm storage matrix: The same arrangement as in the 3 × 10, Z_h = -160 cm storage matrix, except that Z_h = -180 cm (or Z_h = 0 cm).

v. Zigzag storage matrix: This is a 2 × 15 storage matrix^b where the casks in the two rows are horizontally misaligned (refer to Figure 5).

3.2.1 The Effect of Burnup on the Correlation between D and k_eff

The role of burnup on the material density was described in section 3.1.2, and is still equally applicable to 30 casks and will therefore not be repeated, as it was already described in detail in that section. The focus of this section will therefore be much more on how, for a given burnup level, the selected storage matrix influences the k_eff.

The objective of changing D by a regular amount for various storage matrices is to determine whether there is a link/correlation between D and k_eff and not so much to show the effect of burnup on k_eff, which had already been done (Ref. [3]). The results of the analysis indicate that for a given burnup regime, the Zigzag storage matrix will have the lowest k_eff irrespective of the value of D (refer to Figure 8).

Low Burnup (Burnup A).

The results in Figure 9a and Table 12 indicate that there is a correlation between D and the k_eff for various storage matrices, and this correlation between the two variables is defined by polynomial equations Eqn. 12-Eqn. 16. The overall trend is quantified by a linear correlation (R) derived from the same data points used for the polynomial fit. The resulting linear relationship describing the general/net trend of k_eff with increasing D, along with its associated statistical parameters, is presented in Table 12.

Figure 9 Comparison of the k_eff of 30 casks of different storage matrices within the same burnup.

Table 12 Statistical data for the linear correlation coefficient (R) for burnup A.

(i) 3 × 10, Z_h = -180: k_eff Decreases as D Increases (for the net trend, refer to Figure 9a and Table 12).

\[ \begin{aligned} y_{3\times10,Z=-180,Burnup\,A}=k_{eff} &=3.35957-0.08797\mathrm{D}+0.0013\mathrm{D}^{2} \\&-1.04997\mathrm{E}-5\mathrm{D}^3+5.07811\mathrm{E}-8\mathrm{D}^4 \\&-1.51211\mathrm{E}-10\mathrm{D}^5+2.71556\mathrm{E}-13\mathrm{D}^6\\&-2.6967\mathrm{E}-16\mathrm{D}^7+1.1367\mathrm{E}-19\mathrm{D}^8 \end{aligned} \tag{Eqn.12} \]

(ii) Zigzag: k_eff Decreases as D Increases.

\[ \begin{aligned} y_{Zigzag,Burnup\,A}=k_{eff} &=2.20447-0.04489\mathrm{D}+6.33028\mathrm{E}-4\mathrm{D}^2 \\&-4.87872\mathrm{E}-6\mathrm{D}^3+2.26049\mathrm{E}-8\mathrm{D}^4 \\&-6.47772\mathrm{E}-11\mathrm{D}^5+1.12532\mathrm{E}-13\mathrm{D}^6 \\&-1.08672\mathrm{E}-16\mathrm{D}^7+4.47683\mathrm{E}-20\mathrm{D}^8 \end{aligned} \tag{Eqn.13} \]

(iii) 3 × 10, Z_h = -160: k_eff Increases as D Increases (refer to Figure 9a and Table 12).

\[ \begin{aligned} y_{3\times10,Z=-160(Burnup\,A)}=k_{eff} &=1.82774-0.033\mathrm{D}+4.80244\mathrm{E}-4\mathrm{D}^2 \\&-3.81389\mathrm{E}-6\mathrm{D}^3+1.81391\mathrm{E}-8\mathrm{D}^4 \\&-5.30915\mathrm{E}-11\mathrm{D}^5+9.37027\mathrm{E}-14\mathrm{D}^6 \\&-9.1459\mathrm{E}-17\mathrm{D}^7+3.79048\mathrm{E}-20\mathrm{D}^8 \end{aligned} \tag{Eqn.14} \]

(iv) 3 × 10, Z_h = -170: k_eff Increases as D Increases (refer to Figure 9a and Table 12).

\[ \begin{aligned} y_{3\times10,Z_{h}=-170,Burnup\,A}=k_{eff} &=-1.05072+0.06695\mathrm{D}-9.64759\mathrm{E}-4\mathrm{D}^2 \\&+7.59689\mathrm{E}-6\mathrm{D}^3-3.59154\mathrm{E}-8\mathrm{D}^4 \\&+1.04782\mathrm{E}-10\mathrm{D}^5-1.84815\mathrm{E}-13\mathrm{D}^6 \\&+1.80678\mathrm{E}-16\mathrm{D}^7-7.51365\mathrm{E}-20\mathrm{D}^8 \end{aligned} \tag{Eqn.15} \]

(v) 2 × 15 Traditional Storage Matrix: The k_eff Increases as D Increases (for the net trend refer to Figure 9a and Table 12).

\[ \begin{aligned} y_{2\times15,Burnup\,A}=k_{eff} &=-0.94199+0.06372\mathrm{D}-9.23956\mathrm{E}-4\mathrm{D}^2 \\&+7.30421\mathrm{E}-6\mathrm{D}^3-3.45768\mathrm{E}-8\mathrm{D}^4 \\&+1.00753\mathrm{E}-10\mathrm{D}^5-1.77095\mathrm{E}-13\mathrm{D}^6 \\&+1.7222\mathrm{E}-16\mathrm{D}^7-7.11429\mathrm{E}-20\mathrm{D}^8 \end{aligned} \tag{Eqn.16} \]

Medium Burnup (Burnup B).

For the mid burnup (burnup B), the results in Figure 9b and Table 13 reveal that there is a correlation between D and the k_eff across all storage matrix configurations. This correlation is described by the polynomial fit—which is an indication of an instantaneous response rather than general or net overall k_eff trend—and is given by Eqn. 17-Eqn. 21. The net trend of the k_eff on the other hand is described by a linear correlation (R), with the corresponding statistical parameters listed in Table 13. These results indicate (based on the negative value of R) that, as D increases the k_eff for 2 × 15, 3 × 10, Z_h = -160 and 3 × 10, Z_h = -180 matrices decreases whereas for 3 × 10, Z_h = -170 and Zigzag matrices the k_eff increases. The statistical data for the linear correlations (R)—which indicates the general/net trend of the k_eff with change in D—for the storage matrices is listed in Table 13.

Table 13 Statistical data for the linear correlation coefficient (R) for burnup B.

(i) 2 × 15 Traditional Storage Matrix. k_eff decreases as D increases (B = -2.45E-06, R = -0.54853) and the relationship between k_eff and D is defined by Eqn. 17.

\[ \begin{aligned} y_{2\times15(Burnup\,B)}=k_{eff} &=-0.05041+0.03101\mathrm{D}-4.35001\mathrm{E}-4\mathrm{D}^2 \\&+3.32008\mathrm{E}-6\mathrm{D}^3-1.51741\mathrm{E}-8\mathrm{D}^4 \\&+4.2753\mathrm{E}-11\mathrm{D}^5-7.28336\mathrm{E}-14\mathrm{D}^6 \\&+6.88368\mathrm{E}-17\mathrm{D}^7-2.77143\mathrm{E}-20\mathrm{D}^8 \end{aligned} \tag{Eqn.17} \]

(ii) 3 × 10, Z_h = -160. For 3 × 10, Z_h = -160, the results indicate that k_eff decreases as D increases, with B and R equal to -3.50E-07 and -0.08319 respectively (refer to Figure 9b and Table 13). The relationship between k_eff and D is captured by Eqn. 18, and the statistical regression data for polynomial correlation are summarised in Table 14.

\[ \begin{aligned} y_{3\times10,Z=-160,Burnup\,B}=k_{eff} &=1.77076-0.0302\mathrm{D}+4.1518\mathrm{E}-4\mathrm{D}^2 \\&-3.10261\mathrm{E}-6\mathrm{D}^3+1.38778\mathrm{E}-8\mathrm{D}^4 \\&-3.82631\mathrm{E}-11\mathrm{D}^5+6.37991\mathrm{E}-14\mathrm{D}^6 \\&-5.904\mathrm{E}-17\mathrm{D}^7+2.32889\mathrm{E}-20\mathrm{D}^8 \end{aligned} \tag{Eqn.18} \]

Table 14 Regression data for polynomial equations that define the correlation between k_eff and D for various storage matrices for Burnup B.

(iii) 3 × 10, Z_h = -180. Similarly, for 3 × 10, Z_h = -180, the results indicate that k_eff decreases as D increases, with B and R equal to -2.96E-07 and -0.42633 respectively. The relationship between k_eff and D is given by Eqn. 19.

\[ \begin{aligned} y_{3\times10,Z_{h}=-180,Burnup\,B}=k_{eff} &=0.11877+0.0258\mathrm{D}-3.72933\mathrm{E}-4\mathrm{D}^2 \\&+2.94944\mathrm{E}-6\mathrm{D}^3-1.4005\mathrm{E}-8\mathrm{D}^4 \\&+4.09863\mathrm{E}-11\mathrm{D}^5-7.2368\mathrm{E}-14\mathrm{D}^6 \\&+7.0659\mathrm{E}-17\mathrm{D}^7-2.92825\mathrm{E}-20\mathrm{D}^8 \end{aligned} \tag{Eqn.19} \]

(iv) In the 3 × 10, Z_h = -170. The correlation between k_eff and D has been found to have a positive R-value, i.e., k_eff increases as D increases, with B = 5.26E-07 and R = 0.13348. The statistical data are listed in Table 13. The relationship between k_eff and D has been found to obey Eqn. 20.

\[ \begin{aligned} y_{3\times10,Z_{h}=-170,Burnup\,B}=k_{eff} &=-0.34518+0.04314\mathrm{D}-6.39899\mathrm{E}-4\mathrm{D}^2 \\&+5.17247\mathrm{E}-6\mathrm{D}^3-2.50227\mathrm{E}-8\mathrm{D}^4 \\&+7.44719\mathrm{E}-11\mathrm{D}^5-1.33631\mathrm{E}-13\mathrm{D}^6 \\&+1.32601\mathrm{E}-16\mathrm{D}^7-5.58667\mathrm{E}-20\mathrm{D}^8 \end{aligned} \tag{Eqn.20} \]

(v) Zigzag. The k_eff of this storage matrix displays a positive correlation coefficient between k_eff and D, with B = 2.53E-07 and R = 0.06048 (refer to Table 13).

The relationship between k_eff and D is given by Eqn. 21; for Variance values refer to Table 14.

\[ \begin{aligned} y_{Zigzag,Burnup\,B}=k_{eff} &=1.08698-0.0079\mathrm{D}+1.16782\mathrm{E}-4\mathrm{D}^2 \\&-9.56768\mathrm{E}-7\mathrm{D}^3+4.75354\mathrm{E}-9\mathrm{D}^4 \\&-1.46735\mathrm{E}-11\mathrm{D}^5+2.74933\mathrm{E}-14\mathrm{D}^6 \\&-2.85968\mathrm{E}-17\mathrm{D}^7+1.26476\mathrm{E}-20\mathrm{D}^8 \end{aligned} \tag{Eqn.21} \]

High Burnup (Burnup C).

The k_eff in burnup regime C is the lowest of all burnups studied, i.e., lower than burnup regimes A and B, which is consistent with decreases in k_eff as burnup is increased. An analysis of the trend between k_eff and D indicates that.

(i) 3 × 10, Z_h = -180. The net trend of the k_eff in this storage matrix decreases as D increases, B = -4.93E-07 and R = -0.097 (refer to Table 15 and Figure 9c).

The corelation between k_eff and D is given by polynomial equation described by Eqn. 22-Eqn. 26:

\[ \begin{aligned} y_{3\times10,Z_{h}=-180,Burnup\,C}=k_{eff} &=0.15982+0.02183\mathrm{D}-2.78452\mathrm{E}-4\mathrm{D}^2 \\&+1.88647\mathrm{E}-6\mathrm{D}^3-7.42906\mathrm{E}-9\mathrm{D}^4 \\&+1.73572\mathrm{E}-11\mathrm{D}^5-2.32547\mathrm{E}-14\mathrm{D}^6 \\&+1.59187\mathrm{E}-17\mathrm{D}^7-3.9746\mathrm{E}-21\mathrm{D}^8 \end{aligned} \tag{Eqn.22} \]

Table 15 Statistical data for the linear correlation coefficient (R) for burnup C.

(ii) Zigzag Storage Matrix. The Zigzag storage pattern follows the same trend as the 3 × 10, Z_h = -180, i.e., k_eff decreases as D increases, with B = -8.7E-08 and R = -0.03165.

The equation for the relationship between k_eff and D is given by Eqn. 23:

\[ \begin{aligned} y_{Zigzag,Burnup\,C}=k_{eff} &=0.82826+0.00219\mathrm{D}-5.36226\mathrm{E}-5\mathrm{D}^2 \\&+6.00639\mathrm{E}-7\mathrm{D}^3-3.66596\mathrm{E}-9\mathrm{D}^4 \\&+1.29766\mathrm{E}-11\mathrm{D}^5-2.66033\mathrm{E}-14\mathrm{D}^6 \\&+2.93041\mathrm{E}-17\mathrm{D}^7-1.34159\mathrm{E}-20\mathrm{D}^8 \end{aligned} \tag{Eqn.23} \]

(iii) 3 × 10, Z_h = -160. The results in Table 15 and Figure 9c indicate that the correlation coefficient between k_eff and D is positive, with B = 5.64E-07 and R = 0.11414. The relationship between k_eff and D is given by Eqn. 24.

\[ \begin{aligned} y_{3\times10,Z=-160,Burnup\,C}=k_{eff} &=-0.0561+0.03356\mathrm{D}-5.16396\mathrm{E}-4\mathrm{D}^2 \\&+4.33334\mathrm{E}-6\mathrm{D}^3-2.17255\mathrm{E}-8\mathrm{D}^4 \\&+6.67996\mathrm{E}-11\mathrm{D}^5-1.23376\mathrm{E}-13\mathrm{D}^6 \\&+1.25547\mathrm{E}-16\mathrm{D}^7-5.40571\mathrm{E}-20\mathrm{D}^8 \end{aligned} \tag{Eqn.24} \]

For the errors and other statistical data for this equation, refer to Table 16.

Table 16 Regression data for polynomial equations that define the correlation between k_eff and D for various storage matrices for Burnup C.

(iv) 3 × 10, Z_h = -170. This storage matrix also suggests that there is a positive correlation between k_eff and D, indicating that the k_eff increases as D increases. The slope (B) and correlation coefficient (R) are 5.25E-07 and 0.10782, respectively. The relationship between k_eff and D is given by Eqn. 25:

\[ \begin{aligned} y_{3\times10,Z=-170,Burnup\,C}=k_{eff} &=2.01701-0.04107\mathrm{D}+6.02497\mathrm{E}-4\mathrm{D}^2 \\&-4.8144\mathrm{E}-6\mathrm{D}^3+2.30019\mathrm{E}-8\mathrm{D}^4 \\&-6.75475\mathrm{E}-11\mathrm{D}^5+1.19527\mathrm{E}-13\mathrm{D}^6 \\&-1.16947\mathrm{E}-16\mathrm{D}^7+4.85968\mathrm{E}-20\mathrm{D}^8 \end{aligned} \tag{Eqn.25} \]

The reader is referred to Table 16 for the statistical data pertaining to this equation.

(v) 2 × 15 Traditional Storage Matrix. This storage pattern also reveals that as D increases, the k_eff will also increase, with B and R being 1.46E-06 and 0.43526 respectively. The mathematical relationship between k_eff and D is defined by Eqn. 26

\[ \begin{aligned} y_{2\times15,Burnup\,C}=k_{eff} &=-0.21353+0.0372\mathrm{D}-5.38824\mathrm{E}-4\mathrm{D}^2 \\&+4.26143\mathrm{E}-6\mathrm{D}^3-2.02224\mathrm{E}-8\mathrm{D}^4 \\&+5.9211\mathrm{E}-11\mathrm{D}^5-1.04831\mathrm{E}-13\mathrm{D}^6 \\&+1.02909\mathrm{E}-16\mathrm{D}^7-4.29905\mathrm{E}-20\mathrm{D}^8 \end{aligned} \tag{Eqn.26} \]

The relationship between k_eff and D of a 2 × 15 storage matrix at burnup A is described by a polynomial equation of the 8^th order presented as in Eqn. 12-Eqn. 16 with coefficients summarised in Table 17.

Table 17 Regression data for polynomial equations that define the correlation between k_eff and D for various storage matrices for Burnup A.

The equations relating the k_eff to D for burnup B are summarised by Eqn. 17 to Eqn. 21 with the coefficients listed in Table 14.

The polynomial equations summarising the relationship between D and k_eff for burnup C are Eqn. 22-Eqn. 26, with their statistical data summarised in Table 16.

The Effect of Combination of Storage Matrix with Burnup Credit in Preserving Spent Fuel Storage. The application of burnup credit to recover storage capacity that would otherwise be lost under fresh fuel scenario was described in Refs. [3,11]. In these Refs, i.e. Refs. [3,11] it was demonstrated that a nuclear criticality safety analysis that takes into consideration the combined effects of burnup credits of major actinides and principal fission products results in the recovery of a large of spent fuel storage capacity. This is because the accumulative neutron absorption cross section of these two nuclide sets is large in comparison to the combined fission cross section. As a result of this, the neutron absorption event exceeds the nuclear fission event. Therefore, when these are taken into account in nuclear criticality calculation, the nuclear operator saves between 11-15% of storage space. This percentage can be increased quite significantly if burnup credit is used in combination with an appropriate storage matrix as indicated in Figure 10—especially if used conjunction with the zigzag storage matrix in an appropriate vertical misaligned configuration.

Figure 10 Effect of combination of spatial arrangement and burnup credit on k_eff.

To determine the effect of raising the middle row with an offset (Z_h), the behaviour of k_eff with increasing Z_h was analysed by lifting the middle row of 3 × 10, Z_h = -180. The analysis was conducted at a constant burnup and constant D (burnup regime B, and D = 100, respectively). The results of the analysis are tabulated in Table 18 and Figure 11.

Table 18 Z-axis elevation.

Figure 11 3 × 10 (30 MWD/MTU) k_eff as a function of the offset of the z-axis.

\[ \begin{aligned} y=k_{eff} &=3209519.78682-189095.02433z_h-4944.41622(z_h)^2 \\&-75.30854(z_h)^3-0.73632(z_h)^4-0.00479(z_h)^5 \\&-2.07667\mathrm{E}-5(z_h)^6-5.77636\mathrm{E}-8(z_h)^7 \\&-9.35916\mathrm{E}-11(\mathrm{z}_h)^8-6.73002\mathrm{E}-14(\mathrm{z}_h)^9 \end{aligned} \tag{Eqn.27} \]

Where y is k_eff and x is the z-coordinate, and the values of the coefficients of this equation (i.e. Eqn. 27) are summarised in Table 19.

Table 19 Regression data of a polynomial equation that shows the correlation between k_eff and Z_h.

3.3 Factors Influencing the k_eff Trend of the Thirty Casks Storage Matrix

3.3.1 Material and Fission Density

The material data discussed in section 3.1.2, which influence the k_eff of storage pattern of four casks, are still applicable to storage of 30 casks in various storage matrices and are discussed in detail in the following paragraph (i.e. paragraphs 3.1.2 and 4.1). This buildup leads to an increase and overlap of resonance absorption cross section, which ultimately results in a decrease in the resonance escape probability of neutrons and, consequently, in a decrease in k_eff.

3.3.2 Choice of a Storage Matrix

The key factors that influence the k_eff when viewed from different storage matrix perspectives, include changes in the effectiveness of spatial shielding—which in some cases can provide partial shielding or complete shielding—due to their spatial configuration within the storage matrix and the self-shielding effects that result from the buildup of fission products caused by burnup.

Furthermore, increasing (D) between adjacent rows of casks increases the area and volume of storage space they occupy, which in turn decreases the minimum critical mass of fissile material required to sustain the chain reaction [15,16,17].

To understand the influence of changing D in a given storage on the k_eff, the solid angle method described by Eqn. 10 and Eqn. 11 in conjunction with Figure 8.

Comparison of Storage Arrays on the Basis of k_eff: Burnup Regime A. In applying the solid angle method (refer to Figure 8), to analyse the reasons for the behaviour of the k_eff in response to increasing the in D, it was observed that for the traditional 2 × 15 storage matrix the most representative central unit is located midway between the two rows, bordered by casks number 7 from each of the two rows.

In that way, H between point P and the most central representative unit represented by Figure 8, increases as D, i.e., the distance between the two rows of casks increases, which in turn decreases $\mathrm{\Omega}_{allowable}$. By applying Eqn. 10, this will give the corresponding k_eff associated with the minimum allowable solid angle and the corresponding distance, D, between the two rows [6,14].

Upon studying the response k_eff to increase in D under the same burnup conditions (i.e., Burnup regime A) and the same storage matrix (i.e. 2 × 15 traditional storage matrix (refer to Figure 4)), it is observed that the magnitudes of k_eff of 2 × 15, traditional storage matrix, 3 × 10 with Z_h = -160 (refer to Figure 6 and Figure 9), 3 × 10 with Z_h = -170, all decrease, with their respect linear regression slope of 3.34E-07, 3.78E-07, 2.11E-06, respectively. The k_eff of 3 × 10 with Z_h = -180 and 2 × 15 Zigzag, on the other hand, were found to be decreasing as D increases, with the respective linear regression slopes being -2.62E-06 and 2.91E-06, respectively (refer to Table 12). The increase in k_eff is because of the increase in H, as shown in Figure 8 and described by Eqn. 11, whenever H refocuses to the middle (i.e., L/2) of the representative cask. The same reasoning—an increase in H—also explains the increase in the k_eff of the traditional 2 × 15 storage matrix.

When the 2 × 15 Zigzag storage matrix is studied, it is observed that k_eff decreases rapidly as D increases. This may be attributed to the fact that the Zigzag configuration is longer than the 2 × 15 traditional matrix (refer to Figure 5). Because of its length–from north to south—H increases with an increase in D, but at a much faster rate because in the Zigzag pattern the distance between adjacent casks is also larger than in the traditional 2 × 15 storage matrix. Furthermore, because of the staggered nature of the array, the 2 × 15 Zigzag casks storage matrix has a higher neutron resonance escape probability than the traditional 2 × 15 storage matrix.

Regarding 3 × 10, with Z_h = -180, the rapid decrease in k_eff is because of the alignment of high burnup regions of adjacent casks. Furthermore, because the casks are positioned on their base level—i.e., flat on the concrete slab—this results in the shorter H and thus the largest $\mathrm{\Omega}_{allowable}$.

Correlations between k_eff and Storage Matrix in the Medium Burnup (Burnup Regime B). When the burnup is increased from lower burnup (burnup regime A) to medium and higher burnup (burnup B and C, respectively), it is observed that when comparing similar storage matrices—for example, a 2 × 15 storage matrix—it is observed that the k_eff decreases as burnup increases. This may be attributed to an accumulation of fission products, which are formed as burnup increases and therefore absorb neutrons that would otherwise contribute to fission and increase the fission density.

The effect of burnup on k_eff was described in Ref. [3], and the combination of burnup and storage pattern is presented in Figure 9 and Figure 10.

When the response of k_eff to changes of various storage matrix is compared under similar medium burnup regime, it is observed that the 2 × 15 traditional storage matrix, 3 × 10 at Z_h = -160 and 3 × 10 at Z_h = -180, show a decrease with an increase in D—based on the linear fit data. However, the 2 × 15 traditional storage matrix has the fastest rate of decrease, with the slope of -2.45E-06 (refer to Table 13). On the other hand, the 3 × 10 at Z_h = -170 and 2 × 15 Zigzag storage matrix show an increase in k_eff with an increase in D. In addition to that, the 2 × 15 Zigzag storage matrix has the lowest k_eff.

The 2 × 15 Zigzag storage matrix has the lowest k_eff because it has the largest H, and also because of the high resonance escape probability of neutrons in this storage matrix. The decrease in k_eff observed on 2 × 15 traditional storage matrix and 3 × 10 at Z_h = -160 is because of an increase in H as D increases. The k_eff of the 2 × 15 traditional storage decreases at a faster rate than 3 × 10 at Z_h = -160, because the resonance escape probability is higher in the 2 × 15 traditional matrix than it is for the 3 × 10 at Z_h = -160 storage matrix. This, combined with the fact that the array for a 2 × 15 traditional matrix is longer (in the north—south direction) than the 3 × 10 at Z_h = -160. Because of this, H increases much more rapidly with increasing D compared to the 3 × 10 at Z_h = -160 matrix.

Correlations between k_eff and Storage Matrix in the High Burnup (Burnup Regime C). On studying the results of the relationship between k_eff and D from Figure 9c, it was observed that the 2 × 15 Zigzag storage matrix has the lowest k_eff among all arrays. In addition to that, the k_eff decreases as D increases. Based on the regression data listed in Table 15, the slope of its linear regression line was found to be -8.7E-08.

The explanation for the low k_eff in the 2 × 15 Zigzag storage matrix is that, in this configuration, the casks are misaligned/staggered, offering a greater opportunity for neutron resonance escape probability. This, combined with the fact that in the 2 × 15 Zigzag storage matrix, the array is horizontally long, requires H (refer to Figure 8, Eqn. 10 and Eqn. 11) to increase significantly as D increases.

Upon further analysis of the behavior of k_eff of other storage patterns, it was found that the k_eff of the following storage arrays increases as D increases: 2 × 15 Traditional storage matrix, 3 × 10 at Z_h = -160, and 3 × 10 at Z_h = -170. In contrast, the k_eff of 3 × 10 at Z_h = -180 was found to decrease as D increases. The rate of change as D increases is represented by the gradient (denoted as parameter B) of their respective linear regression lines, which are listed in Table 15.

On analysing the rationale for the increase in k_eff for the 2 × 15 Traditional storage matrix, it was found that it may be attributed to a horizontally long array, which requires H to increase as D increases. With regard to the increase in k_eff of 3 × 10 at Z_h = -160 and 3 × 10 at Z_h = -170, this can be ascribed to the fact that the storage is compact and therefore experiences a higher fission density. In addition to that, because the casks are elevated by 20 and 10 cm, respectively, it implies that in order for P to maintain the focal point at L/2 as shown in Figure 8, P must refocus by increasing both the vertical and horizontal distances of H as D increases.

On the 3 × 10 at Z_h = -180, on the other hand, because the casks are resting on their base level, the increase in D has no vertical effect on H; only the horizontal distance will increase as D increases.

4. Analysis

When the number of casks is increased from four, where the storage pattern is either 1 × 4 or 2 × 2, to 30, where the storage pattern is primarily 2 × 15 or 3 × 10, and within these patterns there are variations as described in the results, the neutron interaction among neighbouring surrounding casks and therefore k_eff at any given time is found to be a function of:

i. fissile material density occupied (g/cm³).

ii. distance of non-multiplying reflector, i.e., the wall.

iii. number of casks.

iv. the type of storage pattern selected.

Considering the mutual interaction/interference of γ-photon and neutron flux among adjacent casks, according to Scotta et al. [18], the probability that an incident neutron with energy E and solid angle Ω will interact with the target molecule and is scattered with energy $E'$ with a solid angle $\mathrm{\Omega}'$, may be expressed by a double differential scattering cross section given by

\[ \frac{d^2\sigma}{d\mit{\Omega} dE}\left(E\to E^{\prime},\mit{\Omega}\to\mit{\Omega}^{\prime}\right)=\frac{\sigma_b}{4\pi kT}\sqrt{\frac{E^{\prime}}{E}}e^{-\frac{\beta}{2}}S(\alpha,\beta) \tag{Eqn.28} \]

Where Ω is defined by

\[ \mit{\Omega}\equiv\iint_S\frac{\hat{n}\delta A}{d^2} \tag{Eqn.29} \]

and

σ_b = characteristic bound cross section of the material.

T = temperature.

k = Boltzmann constant.

δA = a small area of a shape.

$\hat{n}$ = unit vector.

d² = distance between the adjacent casks (emitter and the receptor).

$S(\alpha,\beta)$ is the thermal scattering function where α and β are dimensionless momentum and energy transfer defined by Eqn. 30 and Eqn. 31, respectively:

\[ \alpha=\frac{E^{\prime}+E-2\sqrt{E^{\prime}E}\mathrm{cos}\theta}{AkT} \tag{Eqn.30} \]

where A is the ratio of the scattering/target molecule and the neutron mass, and cos θ is the cosine of the scattering angle in the laboratory system [19],

\[ \beta=\frac{E^{\prime}-E}{kT} \tag{Eqn.31} \]

Thus, it stands to reason that a change in geometry, i.e., the spatial arrangement of spent fuel casks, will have a corresponding effect on the k_eff because of the change in Ω and A. The solid angle Ω, which is the basis of the inverse square law, is defined by Eqn. 29.

Given that the spent fuel cask is a cylinder, and the area of a cylinder (A_cyl) is defined by

\[ A_{cyl}=2\pi rh+2\pi r^2 \tag{Eqn.32} \]

where for this particular cask r = 144.3 cm is the radius and h = 423.6 cm is the height of the cask Eqn. 29, may be modified as follows:

\[ \mit{\Omega}_{cyl}\equiv\iint_S\frac{\hat{n}\delta A_{cyl}}{d^2} \tag{Eqn.33} \]

By substituting Eqn. 32 into Eqn. 33 the solid angle of the cylinder ($\Omega_{cyl}$) is obtained.

\[ \mit{\Omega}_{cyl}=\iint_S\frac{\hat{n}\delta(2\pi rh+2\pi r^2)}{d^2} \tag{Eqn.34} \]

Given that for the casks under study (i.e., Castor X/28), r and h are 144.3 cm and 422.4 cm, respectively. Subsequently, by substituting the z and r into Eqn. 34 then $\Omega_{cyl}$ is obtained.

\[ \mit{\Omega}_{cyl{\_}CastorX/{28}}=\iint_S\frac{\hat{n}\delta(2\pi(144.3)422.4+2\pi(144.3)^2)}{d^2} \tag{Eqn.35} \]

Thus, from Eqn. 35, it is evident that Ω is directly proportional to the unit vector and inversely proportional to the distance between casks that interfere with each other.

4.1 Increase in Resonance Cross Section as a Result of Irradiation

The decrease in k_eff observed after the fuel assemblies have been exposed to a certain degree of burnup (see Table 1) or after a period of storage is due to the respective radiation-induced or decay-induced build-up of certain nuclides, actinides, and fission products (refer to Ref. [19]). When the fuel is ‘fresh’ from the manufacturer, it is reasonable to expect that for the UO₂ fuel, the fuel component with a significant amount will be ²³⁵U and ²³⁸U, in which case the cross-section linking will be as indicated in Figure 12b. After the fuel has been irradiated, there will be a build-up of actinides, some of which include ²³⁹Pu, ²⁴⁰Pu, ²⁴¹Pu, and ²⁴²Pu, and fission products (refer to Figure 12a and Figure 12b respectively). The build-up of these nuclides will increase the effective cross section of the system as a result of the interaction/cross-linking of the resonance cross section of isotopes with cross sections that lie in close proximity to one another in the energy spectrum (refer to Figure 11) [20].

Figure 12 Cross-linking of resonance cross sections.

Because of this cross-linking, the net effect on self-shielding by neutron absorption is much higher than if the resonance peaks were farther apart. Had the peaks been farther apart, the neutron would have had a much higher probability of escaping resonance capture than if they were closer together.

Therefore, it is because of the cross-linking of the resonance cross section of actinides and fission products that the k_eff will decrease when spent fuel is stored for long periods or it is subjected to a certain degree of irradiation [19,21].

4.2 Spatial Shielding

The effectiveness of shielding casks from adjacent cylindrical containers is a function of the magnitude of a solid angle that subtends a given point of interest, which in turn is a function of the distance between the neighbouring casks. In an array of storage cylinders containing fissile material, the angular position of the casks determines the magnitude of the area of the cask ‘seen’ by other casks and therefore the area available for neutron interaction. Depending on the amount of shadowing/shielding each cylinder provides, the shielding effectiveness is described as either a grey Dancoff factor or a black Dancoff factor when the shielding is partial or complete, respectively [22,23].

4.2.1 Four Casks

Considering the spatial arrangement of the four casks (refer to Figure 2 and Figure 3) and comparing their respective mutual spatial shielding they provide for one another, in the case of Figure 3, not a single cask is completely shielded from the neutron flux of either, whereas in Figure 2, the middle casks completely shield the outer casks from one another. Because of this arrangement, the storage arrangement in Figure 3 tends to form a flux trap and trap the heat flux between the four casks. The k_eff of the two storage matrices plotted against the separation distances between the casks shows the following (refer to Figure 7):

Burnup Regime A. When the fuel is underburned (refer to Table 20), cumulative burnup = 7537.74 MWD) the k_eff for the linear storage array (i.e. 1 × 4) starts higher than that of the 2 × 2 array when the casks are in contact, and decreases rapidly as the distance between adjacent casks is increased (refer to Figure 7a and Table 6) to below that of the 2 × 2 array. By contrast, in the case of the 2 × 2 array, the k_eff increases steadily (B = -2.48E-06 and 2.07E-08 for 1 × 4 and 2 × 2, respectively). The decrease in k_eff with increased distance in the 1 × 4 array is consistent with the total absorption and total fission of the system, which increase as the distance increases, but absorption is higher than fission (refer to Figure 12a).

Table 20 Implicit burnup regime A read from node 9.

Burnup Regime B. When the fuel burnup at node 9^c increased to an accumulative burnup of 15649.2 MWD/MTU (refer to Table 21), the k_eff tends to increase as the separation distance between casks increases. Contrary to what was observed with burnup A, in this case, i.e., burnup regime B, the k_eff of the 1 × 4 array increases with an increase in the separation distance with B = 3.77E-07, whereas with the 2 × 2 storage array the k_eff increases at a much slower rate with B = 2.07E-08 (refer to Figure 12b and Table 7). The increase in k_eff in the 1 × 4 and 2 × 2 arrays is consistent with the expected decrease in system total fission and system total absorptions, which is the required downward trend for the k_eff to increase, as observed in Figure 12b.

Table 21 Implicit burnup regime B read from node 9.

Burnup Regime C. When the fuel burnup at node 9 increases to an accumulative burnup of 19365.5 MWD/MTU (refer to Table 22), it is observed that when the casks are in contact, the k_eff of the 1 × 4 array is below that of the 2 × 2 array (at 0.8538 and 0.8541 respectively) but the k_eff of the 1 × 4 array increases rapidly (B = 1.72E-06), whereas the k_eff of the 2 × 2 array decreases (B = -5.65E-07) (refer to Figure 7c). These changes in the k_eff in the 1 × 4 and 2 × 2 arrays are in line with an expected decrease in system total fission and system total absorption, which are the conditions needed for the k_eff to increase, as observed in Figure 12c.

Table 22 Implicit burnup regime C read from node 9.

The increase in burnup from burnup regime A to burnup regime C leads to the buildup of Fission Products, which will result in the overlap of the resonance absorption cross-section of various nuclides. Because of this, the neutrons will have a lower resonance escape probability and consequently a decrease in k_eff. This factor is equally applicable to both 4 and 30 casks.

4.2.2 Thirty Casks

Considering Figure 13, it is noted that there is no shielding between casks 1 and 4 and casks 1 and 3. The neutron flux interaction between casks 1 and 3 and casks 1 and 4 is completely unhindered/unshielded, and therefore, this pair, i.e., casks 1 and 3 and casks 1 and 4, can interact directly with each other with either one acting as the radiation source. However, when one looks at the interaction between cask 1 and cask 5, it is noted these two, i.e. casks 1 and cask 5, are completely shielded from each other by cask 3 which represent a black Dancoff Factor, and casks 1 and 2 are partly shielded from each other by casks 3, 5 and 6 (representing the grey Dancoff factor).

Figure 13 Partial and complete shielding.

Because the casks are cylinders with radius r and distance d apart, for the interaction situation where there is no shielding between interacting fissile cylinders, as is the case of casks 1 and 3, casks 1 and 4, casks 2 and 6, casks 4 and 6 and casks 3 and 5, the Dancoff factor^d of such cylinders may be derived from Pershagen [24]. The shielding/shadowing factor among them is given by the Dancoff factor C and is described by Eqn. 36, which is in accordance with Pershagen's formalism [24]. According to Rodrigues et al. [23], for perfectly absorbing material, as is the case between casks 1 and 3 and casks 1 and 4, C may be calculated using escape probability as follows:

\[ C=1-\frac{P_F}{P_0} \tag{Eqn.36} \]

where P_F is the probability that a neutron origination in the fuel region will make it to the moderator without collision/absorption.

P₀ is the escape probability from fuel.

For partially absorbing/shielding material such as casks 3, 5 and 6 which only provide partial shielding for casks 2 from 1, the Dancoff factor may be calculated by

\[ C_1=\frac{\left(\Sigma P_{1,3}\right)-P_{1,1}{^{\prime}}}{1-P_{1,1}{^{\prime}}} \tag{Eqn.37} \]

where

P_1,3 = the probability of a neutron escaping isotropically from 1 and making its first collision in fuel 3.

$P_{1,1}{^{\prime}}$ = is the self-collision probability in 1 (i.e., neutron originating from cask 1 and being absorbed in cask 1).

The probability of absorption is further increased by the burnup, which, as stated earlier, often leads to a build-up of nuclides with a large cross section for the absorption of thermal neutrons, e.g., ¹⁵⁵Gd. The adverse effect of increased burnup is a decrease in the k_eff. The relationship between D and the k_eff under different burnup regimes is described in the following sections.

Burnup Regime A. The results show that the k_eff of all storage matrices fluctuates with an increase in the distance between adjacent rows. In addition, it is also observed that the 2 × 15 matrix and the 3 × 10 matrix are both hovering around the same k_eff band (circa 0.882 and 0.885, with 3 × 10, Z_h = -180, which dictates the upper and lower limit band width). Their max and min k_eff are as follows (refer to Figure 9a):

i. 3 × 10, Z_h = -180 (max k_eff) = 0.885 at D = 300 and (min k_eff) = 0.882 at D = 50.

ii. 3 × 10, Z_h = -170 (max k_eff) = 0.884 at D = 350 and (min k_eff) = 0.882 at D = 200.

iii. 3 × 10, Z_h = -160 (max k_eff) = 0.884 at D = 100 & at 450; and (min k_eff) = 0.883 at D = 250.

iv. 2 × 15 traditional storage array: (max k_eff) = 0.884 at D = 248.773 and (min k_eff) = 0.882 at D = 300.120.

The k_eff of the Zigzag storage matrix is significantly lower than that of the 2 × 15 and all variations of the 3 × 10 storage matrices (refer to Figure 9a). One of the critical factors leading to the decrease in k_eff as D increases in the Zigzag storage matrix is partial shielding (refer to Figure 5), an effect that was first reported by Dowson and Abbey in Ref. [6]. Because of partial shielding, some neutrons do not interact with the target fissile material to cause fission; instead, they are absorbed or leak out of the system. The mathematical relationship between escape of neutrons and partial shielding is described in section 4.3.3, and the equation that provides the relationship between the escape probability and the Dancoff factor due to partial shielding is given by Eqn. 37. With regard to the general trend in the relationship between D and k_eff, it is observed that:

a. The k_eff of 3 × 10, Z_h = -180; decreases as D increases with B = -2.62E-06 (refer to Table 12). This (i.e., the decrease in k_eff as D increases) is a result of a decrease in fissile material density (g/cm³) as the volume increases, owing to an increase in D, which is consistent with the increase in the number of fissions and the number of absorptions depicted in Figure 13.

b. For 3 × 10, Z_h = -160 and 3 × 10, Z_h = -170; the results show that when the height (also referred to as the offset) of the middle row of the 3 × 10 array is increased by 10 and 20 cm, i.e. Z_h = -170 and Z_h = -160 respectively, and D is increased, it was discovered that the k_eff increases which is contrary to when all the casks had a zero offset, i.e. Z_h = -180. However, 3 × 10 Z_h = -170 (B = 2.11E-06) increases at a much faster rate^e than 3 × 10, Z_h = -160 (B = 3.78E-07). This increase in k_eff with an increase in Z_h is because of the following:

• Absence/lack of spatial shielding on the fraction of the casks that have been elevated (refer to Figure 6).
• Decrease in energy self-shielding (refer to Figure 11) on the vacant bottom part left by casks that have been elevated.
• Misalignment of highly reactive regions (bottom and top of the casks) in the middle row and the two outer rows. By lifting the middle row, the bottom reactive row becomes aligned with the less reactive middle part of the casks in the outer rows. The amount of under-reacted material on the top and bottom parts of the fuel assembly, exposed to the over-reacted middle part of the fuel assembly in the two outer rows, determines the k_eff. It is because of that Z_h = -170 changes at a higher rate than Z_h = -160.
• When the offset, i.e., Z_h, is increased from 0 to 20 as in 3 × 10, Z_h = -180, to 3 × 10, Z_h = -160, respectively, the k_eff will increase; thus, increasing the offset beyond that will result in a decrease in k_eff (refer to Figure 10, Table 18, and Table 19). Having studied the fluctuation of k_eff with increases in offset (Z_h), it has been discovered that the lowest k_eff occurs when Z_h = 16.75, followed by Z_h = 41.59, where the k_eff are 0.83354 and 0.83386, respectively.
• The k_eff of the Zigzag matrix is significantly lower than that of all the other storage matrices and decreases even more as D increases. The contributory factors which make the k_eff of the Zigzag storage matrix to be lower than the others include the fact that each individual cask on either row is at an angle with respect to its counterpart on the other row of casks which affects the solid angle as the number of neutrons emanating one particular casks in a given row with a requite energy reaching the other cask on the opposite row is affected by the angular position of the casks with respect to one another, which has an effect on the neutron flux. The reactivity effects of the interaction neutrons between two containers have also been quantified by Kotaro Tonoike et al. [25] where it was discovered that if a single unit containing fissile material is located near a similar unit with the same shape and both contain fissile material, then the relationship between the change in reactivity and the solid angle (Ω) (refer to Eqn. 28 to Eqn. 35) may be given by ${\frac{\Delta\mathbf{k}}{\mathbf{k}}}\Omega={\frac{\Omega}{1/\mathbf{M}^{2}\mathbf{B}^{2}+1-\Omega}}$, where M² is the neutron migration area, B² is the geometric buckling of a single unit and Ω is the solid angle between the units [25]. This equation i.e, ${\frac{\Delta\mathbf{k}}{\mathbf{k}}}\Omega={\frac{\Omega}{1/\mathbf{M}^{2}\mathbf{B}^{2}+1-\Omega}}$ is consistent with Eqn. 10 when applied in combination with Eqn. 11.

Burnup Regime B. When the burnup is increased from A to B, there will be an increase in the material density of some nuclides. For example, the actinides that will increase in material density, which are included in the actinide nuclide set, are ²³⁸Pu, ²³⁹Pu, ²⁴¹Pu, and ²⁴¹Am. The fission products that will increase as a result of burnup are ²⁴³Am, ²³⁷Np, ¹⁴³Nd, ¹³³Cs, ¹⁵⁵Gd, ¹⁵¹Sm, ⁹⁹Tc, ¹⁵³Eu, ¹⁴⁷Sm, ¹⁴⁵Nd, ¹⁵⁰Sm, ¹⁵²Sm [19] because of the build-up of these nuclides and the fact that the effective total cross section is due to cross-linking of respective individual neutron absorption and fission cross-section, the resonance escape probability of a neutron decreases because of the compounding effect of coupling/cross-linking.

Because of differences in burnup between burnup regimes A, B, and C, the neutron flux will be affected and consequently the k_eff.

As stated earlier, the results show that the k_eff of the respective storage matrices fluctuates; however, the peaks and troughs of burnup B differ from those of burnup A (median, circa 0.86 and 0,88, respectively) and are not synchronised with one another. The min and max k_eff of the respective storage matrices at their corresponding locations are as follows;

i. 3 × 10, Z_h = -180 (max k_eff) = 0.864 at D =150 and (min k_eff) = 0.862 at D = 300.

ii. 3 × 10, Z_h = -170 (max k_eff) = 0.864 at D = 450 and (min k_eff) = 0.862 at D = 150.

iii. 3 × 10, Z_h = -160 (max k_eff) = 0.864 at D = 200 and (min k_eff) = 0.8625 at D = 400.

iv. 2 × 15 traditional storage array: (max k_eff) = 0.8641 at D = 100 and (min k_eff) = 0.86247 at D = 450.

The general trend was obtained by performing a linear fit across all storage matrices, and the results are summarised in Figure 9b and Table 13. The analysis proves the following:

(i) 3 × 10, Z_h = -180. The k_eff decreases with an increase in D, with B = -2.96E-07. This is driven by the rise in total fission, together with the increase in the system's total absorption.

(ii) 3 × 10, Z_h = -170. k_eff increases with an increase in D, B = 5.26E-07. The factors behind the behaviour of the k_eff are the increase in both systems' total fissions.

(iii) 3 × 10, Z_h = -160. The k_eff decreases as D increases with slope B = -3.50E-07. The cause of the behaviour of the k_eff as D increases is the decrease in both systems' total fissions together with the reduction of the system's total absorption.

(iv) 2 × 15. There is a rapid decrease in k_eff as D increases (B = -2.45E-06). This is driven by a reduction in the amount of fission, combined with a decrease in the system total absorptions.

(v) Zigzag. The k_eff decreases as D increases, with B = 2.53E-07. The driving force behind this behaviour is the increase in both system total fission and system total absorption, combined with the fact that H increases quite significantly as D increases.

Furthermore, the Zigzag matrix was found to be significantly lower k_eff than the other storage matrices, as was the case in burnup regime A. This is caused by the long storage matrix, which results in an increase in H as D increases, and also because of the higher neutron resonance escape probability due to the staggered configuration of the Zigzag storage pattern.

Burnup Regime C.

When the explicit burnup^f was increased from burnup B to the implicit burnup regime C, which was taken at node 9 (refer to Table 22), the k_eff will change accordingly as follows:

(i) 3 × 10, Z_h = -180. k_eff (max) is 0.8554 at D = 250 and k_eff (min) is 0.8533 at D = 200. When looking at the correlation between k_eff and D, it has been discovered that k_eff decreases as D increases. The slope (B) of the equation describing the relationship between k_eff and D was found to be -4.93E-07 (refer to Figure 9c, Table 15, and Eqn. 20). The behaviour of k_eff as D increases in this burnup regime is driven by the following factors

• A decrease in the system's total fissions.
• A decrease in the system's total absorptions.

(ii) 3 × 10, Z_h = -170. The maximum k_eff (k_eff max) was found to be 0.8551, occurring at D = 300, and (k_eff min) 0.85319, at D = 250 (refer to Figure 9c). The linear correlation performed to determine the relationship between k_eff and D in this storage matrix indicates that k_eff increases as D increases with B = 5.25E-07 (refer to Figure 9, Table 15, and Eqn. 23). It has been found that the underlying causes of the behaviour of k_eff and D in this storage matrix are:

• a decrease in the system's total fissions.
• a decrease in the system's total absorptions.

Of these two, the dominant factor was found to be the system's total fissions.

(iii) 3 × 10, Z_h = -160. Like other storage matrices, the k_eff of this storage matrix also fluctuates with the location of the casks in space. The maximum k_eff (k_eff (max)) was found to be 0.8550 at D = 450, and k_eff (min) was found to be 0.8530 at D = 300. The linear correlation analysis indicates that k_eff increases with an increase in D, with B = 5.64E-07 (refer to Figure 9, Table 15, and Eqn. 6).

(iv) After further investigation to determine the underlying factors, they were found to be the system's total fissions and the system's total absorptions, which varied with D as follows;

• Increase in systems' total fissions.
• Increase in systems' total absorptions.
• Total fission was higher than total absorptions and also increasing at a much faster rate than total absorptions.

(v) 2 × 15 Traditional Storage Matrix. k_eff (max) was found to be circa 0.8642 at D = 100, and the k_eff (min) was 0.8625 at D = 450 (refer to Figure 9c). A linear correlation analysis indicates that k_eff increases with an increase in D with B = 1.46E-06 (refer to Figure 9, Table 15, and Eqn. 6). The factors which were found to be the cause of the behaviour of k_eff relative to D were the system's total fissions and system's total absorptions, which themselves varied with D as follows:

• Increase in systems total fissions (refer to Figure 14c).
• Increase in systems total absorptions (refer to Figure 14c).
• Total fission was higher (larger y-intercept, relative to that of total absorptions) than total absorptions and also increased faster than total absorptions (refer to Figure 14c).

Figure 14 Systems fission and absorption for thirty (30) casks as a function of the separation gap between adjacent rows.

(vi) Zigzag Storage Matrix. The k_eff of this storage matrix is significantly lower than that of the 2 × 5 matrix and all variations of 3 × 10. One reason for this is the deliberate departure from the size of the building that was used for the 2 × 15 array and all variations of 3 × 10, increasing the dimensions so as to accommodate the Zigzag pattern which could include all 30 casks in this storage pattern within the same building, as used in the analysis of 2 × 15 and all variations of 3 × 10. The dimensions (arbitrary) of the building for the Zigzag storage matrix were as follows:

com = "cask storage building"

Because of the increase in volume, the lattice density^g in this matrix was lower than in the case of the 2 × 15 array and all variations of the 3 × 10. The result of an increase in volume is a decrease in backscattering from the walls because the space in which this simulation took place was larger.

Furthermore, because the casks, when viewed in columns, were not in direct line of sight of each other, i.e., they were situated at angles to each other, thus providing only a partial shielding for each, a grey Dancoff factor [7]. Therefore, the combined effect of this, i.e., the decrease in lattice density, partial shielding, and the reduction in backscattering, leads to a decrease in k_eff, as observed in Figure 9c.

The k_eff (max) was found to be 0.8522 at D = 200 and k_eff (min) 0.8510 at D = 300. A linear correlation analysis indicates that the k_eff decreases with an increase in D, with B = -8.7E-08. The relationship between D and the system's total fission and system total absorption was found to decrease with an increase in D. The y-intercepts for the system's total fission and system total absorption were found to be 6.449E-5 and 4.750E-5, respectively (refer to Figure 14).

4.3 Fissile Material Density

4.3.1 Four Casks

The area occupied by the four casks was approximated by calculating the area of the rectangles enclosing these casks (see Figure 15). The casks in the 2 × 2 array, which are in contact, i.e., d = 0 (refer to Figure 15a), are enclosed by a square whose sides are formed by 2× the diameter of the cask (i.e., 2 × Ф = 2 × 288.6 = 577.2 cm²). The sides of the square were increased by the value of d as the distance between the two casks was increased by d. Thus, the area of the square (A_square) in Figure 15a is $A=(577.2)^2$ given that $A_{square}=L^{2}$. For the 1 × 4 storage array (refer to Figure 15b), the approximate area was calculated by calculating the area of the rectangle circumscribing the cylinders (refer to Figure 15b), where the length (L) of the rectangle is 4× Ф and the breadth (B) is Ф. The volume of the region circumscribing the cask was calculated considering the height of the rectangle = Z-axis = 423.6 cm.

Figure 15 Approximation of area occupied by four casks.

Given that the separation distance was in one direction only, L was increased by d while the breadth was kept constant. When storage arrays 1 × 4 and 2 × 2 are compared based on the area over which the fissile material is distributed, it is noted that the area of the rectangles that circumscribe the casks increases accordingly. Thus, for the same separation distance, e.g., 250 cm, the 2 × 2 storage array is spread over a much larger area than the 1 × 4 array (3.81E06 cm² and 7.30E5 cm², respectively), which makes resonance self-shielding less effective for the 2 × 2 than the 1 × 4 array. As a result of this, for burnup regime A, there is a decrease in k_eff for 1 × 4 (B = -2.48E-06) with increasing d, whereas for the 2 × 2 array, the k_eff increases (B = 2.07E-08) as d increases. For burnup A, the optimum burnup where k_eff (1 × 4) = k_eff (2 × 2) is circa 0.8834 and occurs at d = 132.828 cm.

When burnup is increased from A to B, it is observed that the y-intercept of 1 × 4 is higher than that of 2 × 2 (circa 0.8633 and 0.8631 respectively) and contrary to the behaviour that was observed in burnup A, 1 × 4 increases with an increase in d (B = 3.77E-07), whereas 2 × 2 is relatively insensitive to the rise in burnup, as the gradient (B = 2.07E-08) is the same as it was in burnup regime A. For burnup B, the optimum burnup where k_eff (1 × 4) = k_eff (2 × 2) is indeterminate.

For burnup regime C with the same geometry as in burnup regimes A and B, it is observed that the y-intercept of 1 × 4 is lower than that of 2 × 2 (0.8536 and 0.8542, respectively). It has also been noted that the k_eff for 1 × 4 increases as d increases (B = 1.72E-06), whereas for 2 × 2 the k_eff decreases gradually with an increase in d (B = -5.65E-07). The optimum distance where k_eff (1 × 4) = k_eff (2 × 2) is 0.8541, found when d = 201.515 cm.

The geometry and material density for all three burnup regimes are the same; the factor in the behaviour of k_eff relative to change in d is a change in resonance self-shielding brought about by the build-up of actinides and fission products as a result of irradiation, which will decrease the resonance escape probability. The key factors leading to a decrease in k_eff when d changes are the fast fission factor (which is the fission process caused by neutrons in the fast energy range), the increase in resonance escape probability, as well as the thermal utilisation factor (which is an indication of how effectively neutron are absorbed by the fuel) in a 2 × 2 array there is a much better chance of increase in k_eff caused fast fission factor and thermal utilisation factor than in the 1 × 4 configuration. Furthermore, spent fuel generates a substantial amount of decay heat as various isotopes decay into lower-energy isotopes. This temperature change causes a shift in the neutron absorption cross-section due to Doppler broadening, which in turn leads to a reduction in the fission rate, as the neutrons required for fission are increasingly absorbed. This effect is particularly significant in storage matrices with limited ventilation—such as the 2 × 2 matrix, which has poorer ventilation compared to the 1 × 4 matrix. As a result, the rise in temperature causes a decrease in the effective neutron multiplication factor (k_eff), as observed in the 2 × 2 configuration.

The increase in k_eff in the 1 × 4 storage matrix as burnup increases from burnup regime B to C is caused by the buildup of ²³⁹Pu and ²⁴¹Pu, which are fissile. Because of the increased fuel material density caused by this buildup, this will increase the fission density, thereby increasing the fast fission factor and thermal utilisation factor. In burnup regime A (see Figure 7a), the burnup is too low for the buildup of these nuclides (i.e., ²³⁹Pu and ²⁴¹Pu) to be significant, which is why the k_eff decreases as observed.

4.3.2 Thirty Casks

Similarly, the storage of 3 × 10 is compact compared to that of the traditional 2 × 15 storage matrix. Thus, the material density for 2 × 15 is spread over a wider area than in the 3 × 10 storage array. Therefore, the integral resonance self-shielding is less of a factor in 2 × 15 than in the 3 × 10 storage array.

The available surface area over which the fissile material is spread becomes even more significant in the Zigzag array—because it decreases the critical mass—than it is in the traditional storage array, primarily because of the absence of the cask in the direct line of sight of the cask in the opposite row.

When the response of k_eff to the perturbation of the D is studied under different burnups, the linear correlation coefficients are found to be as follows:

Burnup Regime A.

(i) 2 × 15 Traditional Storage Array. The k_eff increases with an increase in D (B = 3.34E-07) and the y-intercept was found to be 0.8834 and k_eff (2 × 15) = k_eff (3 × 10, Z_h = -170) = 0.8835, which was at D = 317.797 cm, k_eff (2 × 15) = k_eff (3 × 10, Z_h = -180) = 0.8835 at D = 424.939 cm.

(ii) 3 × 10, Z_h = -160. The k_eff for 3 × 10, Z_h = -160 has been found to increase as D increases with B = 3.78E-07 and the y-intercept = 0.88379. The k_eff (3 × 10, Z_h = -160) was found to be equal to k_eff (3 × 10, Z_h = -180) = 0.8839 at D = 297.354 cm and also k_eff (3 × 10, Z_h = -160) = k_eff (3 × 10, Z_h = -170) = 0.8839 occurring at D = 525.348 cm.

(iii) 3 × 10, Z_h = -170. When D is increased, it is noted that the k_eff also increases. The linear correlation of the k_eff for 3 × 10, Z_h = -170 was found to be the same as that for 3 × 10, Z_h = -180, i.e. k_eff (3 × 10, Z_h = -170) = k_eff (3 × 10, Z_h = -180) = 0.88364 at D = 387.064 cm and k_eff (3 × 10, Z_h = -170) = k_eff (2 × 15) = 0.88349 occurring at D = 321.765 cm.

(iv) 3 × 10, Z_h = -180. The k_eff for this storage matrix was found to decrease with an increase in D (R = -0.4263 and B = -2.62E-06). It has also been found that k_eff (3 × 10, Z_h = -180) = (3 × 10, Z_h = -160) = 0.8838 at D = 298.316; k_eff (3 × 10, Z_h = -180) = k_eff (3 × 10, Z_h = -170) = 0.8836 at D = 389.947 and (3 × 10, Z_h = -180) = k_eff (2 × 15) = 0.88352 at D = 426.984 cm.

(v) Zigzag. The k_eff is generally significantly lower than that of 2 × 15 traditional storage matrices and all variations of 3 × 10. The linear correlation coefficient indicates that k_eff decreases with an increase (R = -0.52714) in D, with B = -2.91E-06 and y-intercept = 0.8813.

Burnup Regime B.

(i) 2 × 15 Traditional Storage Array. The k_eff for the 2 × 15 array in this burnup regime decreases with increasing D (B = -2.45E-06, R = -0.548). It has also been found that k_eff (2 × 15) = k_eff (3 × 10, Z_h = -170) = 0.8635 at D = 185.401; k_eff (2 × 15) = k_eff (3 × 10, Z_h = -180) = 0.8632 at D = 311.063; k_eff (2 × 15) = 3 × 10, Z_h = -160) = 0.8630 at D = 390.909.

(ii) 3 × 10, Z_h = -160. The k_eff of the 3 × 10, Z_h = -160, decreases gradually (B = -3.50E-07) with an increase in D (R = 0.08319, y-intercept = 08631).

(iii) 3 × 10, Z_h = -170. The k_eff increases with D (B = 1.26E-06). Furthermore, it has also been discovered that the k_eff (3 × 10, Z_h = -170) = k_eff (2 × 15) = 0.8635 at D = 187.325 cm.

(iv) 3 × 10, Z_h = -180. The k_eff decreases as D increases with B = -2.62E-06 and R = 0.4263. It has also been observed that k_eff (3 × 10, Z_h = -180) = k_eff (2 × 15) = 0.8632 at D = 311.063.

(v) Zigzag. The k_eff of this storage matrix is visibly lower than that of the 2 × 15 traditional storage matrix and all variations of the 3 × 10 storage matrices. The k_eff increases with D, albeit at a slower rate (B = 2.53E-07, R = 0.06048). The linear fit did not intersect any other linear fit line.

Burnup Regime C

(i) 2 × 15 Traditional Storage Array. The k_eff of this storage matrix increases as D increases, with B = 1.46E-06 and R = 0.435. It has also been found that k_eff (2 × 15) = k_eff (3 × 10, Z_h = -170) = 0.8539 at D = 137.662; k_eff (2 × 15) = k_eff (3 × 10, Z_h = -160) = 0.8541 at D = 324.65; k_eff (2 × 15) = k_eff (3 × 10, Z_h = -180) = circa 0.8542 and D = circa 398.725 cm.

(ii) 3 × 10, Z_h = -160. The k_eff increases with an increase in D at a rate of B = 1.26E-06. The analysis also proves that the k_eff (3 × 10, Z_h = -160) = k_eff (2 × 15) = 0.8541 at D = 319.6; k_eff (3 × 10, Z_h = -160) = k_eff (3 × 10, Z_h = -180) = 0.542 at D = 452, k_eff (3 × 10, Z_h = -170) = 0.8538 at D = 122 cm.

(iii) 3 × 10, Z_h = -170. The k_eff increases at a rate of B = 1.26E-06 and k_eff (3 × 10, Z_h = -170) = k_eff (2 × 15) = 0.8539 at D = 126.359.

(iv) 3 × 10, Z_h = -180. The k_eff decreases with an increase in D, with B = -4.93E-07. It has also been found that k_eff (3 × 10, Z_h = -180) = k_eff (2 × 15) = 0.8541 at D = 318.037; k_eff (3 × 10, Z_h = -180) = k_eff (3 × 10, Z_h = -160) = 0.85426 at D = 442.616 cm.

(v) Zigzag. The k_eff of this storage matrix decreases as D increases (B = -0.87E-08 and R = 0.0316).

4.3.3 The Role of the Storage Matrix in Influencing the k_eff of the Fissile System

The two key factors that influence the k_eff of the fissile system—and ultimately the selection of the storage matrix—are the fast non-leakage probability and the thermal non-leakage probability. However, these factors do not act in isolation; they must be considered together with other parameters, such as shielding, material composition, etc, to understand how they collectively influence the k_eff of the storage matrix [1].

When examining Figure 10, it is observed that for all burnup regimes (i.e., Burnup A, B, and C), the 2 × 15 Zigzag storage matrix has the lowest k_eff (refer to Figure 5). For burnup regime A, when the traditional 2 × 15 storage matrix is compared with the 2 × 15 Zigzag, it is observed that the k_eff of the 2 × 15 Zigzag decreases as the D increases, whereas that of the traditional 2 × 15 increases (refer to Table 12 for the slopes of the linear fit). The difference in k_eff between these storage matrices is because the Zigzag storage matrix allows greater neutron escape—i.e., high neutron leakage—whereas in the traditional storage matrix, there is less leakage because of the strong coupling between adjacent reactors.

4.4 Reflectors

The design of castor X/28 includes polyethylene rods on the outer metal shell for neutron thermalisation. Furthermore, in this analysis, the cask is modelled as if it were filled with water to determine the effects of abnormal operating conditions when flooded.

4.4.1 Proximity of a Non-Multiplying Reflector

Because the casks are modelled as stored within the cask storage building, the walls of the building act as a non-multiplying reflector. This has the propensity to produce backscattering neutrons, which leak out of the casks and are scattered from the walls back into other casks and thus add positive reactivity to the fissile material. However, if the reflector, i.e., the walls, is moved further away from the casks, a decrease in k_eff results, as observed in the Zigzag storage pattern relative to the traditional 2 × 15 and 3 × 10 storage matrices.

Furthermore, in the case of the Zigzag storage pattern, the absence of casks in direct line-of-sight in the opposite row in k_eff because of the reduced interaction rates between opposite casks.

4.4.2 High-Multiplying Reflector

Because the cask is made up of polyethylene rods and is also flooded with water, the neutron flux of thermal neutrons backscattered to the fuel assembly will be maximised thereby increasing k_eff. However, given that the concentration of hydrogen is time-dependent, it follows that the k_eff will also fluctuate between any given points.

4.4.3 Alignment of High Neutron Flux Region with Underburned Region

Considering that the neutron flux is unevenly distributed along the length of the cask (refer to Figure 16), when the there is no off set included as in 3 × 10, Z_h = -180, the high neutron flux density of the two casks will be perfectly aligned as is the case with 1 and 1’ and 3 and 3’ and the low flux region will also be aligned as in 2 and 2’. However, when the adjacent casks are misaligned by including an offset as in the case of 3 × 10, Z_h = -160 and 3 × 10, Z_h = -170, the neutron-flux rich region (depleted of fissile material) on one cask in the opposite casks inadvertently becomes aligned with the underburned, highly reactive region, rich in fissile material of its neighbour as is the case with (1 and 2’) and (2 and 3’) of casks A and B respectively (refer to Figure 16).

Figure 16 Two misaligned adjacent casks with axial flux distribution showing the cause of fluctuation of k_eff with changes in Z_h.

Because of this (i.e., the flux-rich region of one cask being aligned with the underburned/highly reactive region of the neighbouring cask), the fission rate will increase, as confirmed by the increase in k_eff in Figure 10. This, i.e., the rise in k_eff, will continue as long as Z_h is increased up to the point where a dip in neutron flux (between nodes 3 and 6) is reached (refer to Figure 16). Increasing Z_h of cask B further will result is a misalignment between the neutron-rich region of one cask and the underburned (neutron depleted) region of a neighbouring cask (2 and 3’). This, i.e., the misalignment, will lead to a decrease in k_eff (refer to Figure 10). The increase in k_eff will resume again later when the second peak in neutron flux is reached (between nodes 6 and 13).

5. Conclusion and Forward Looking

There is a direct correlation between the k_eff and the storage pattern of spent fuel casks in a spent fuel storage facility, which is defined by a polynomial equation. This relationship is directly influenced by the solid angle between a neutron-transmitting cask and a receiving cask. Furthermore, as shown in section 3.1.2, the calculation of the k_eff in an array of fissile components, such as spent fuel casks, takes into consideration the diameter and length of the cask, as well as the distance between the transmitting and receiving casks. The integration of these parameters is shown in Figure 8, from which Eqn. 10 and Eqn. 11 are derived, and the correlation between the solid angle is based.

Furthermore, the behaviour of k_eff relative to changes in D (i.e., the distance between the rows of casks) is not cast in stone; it depends on the burnup the material is exposed to. For example, for a 2 × 15 traditional storage matrix, the k_eff increases with an increase in D for burnup regimes A and C, and yet for burnup regime B, it decreases as D increases. This may be ascribed to the ratio of fission to absorption rate, and also the self-shielding resulting from the buildup of nuclides with large neutron-absorption cross sections due to burnup, combined with changes in the solid angle of the reactive unit [1].

Similarly, in the case of four casks, the optimum distance, which is independent of the storage array, is a function of the burnup (refer to 1 × 4 versus 2 × 2 storage matrices). For example, with the 1 × 4 storage array the k_eff decreases with an increase in D for burnup regime A but increases with an increase in D for burnup regimes B and D. For the 2 × 2 storage matrix, on the other hand, the k_eff rises with an increase in D for burnup regimes A and B but decreases as D increases when the burnup regime is increased to C.

In the 3 × 10 storage array, positioning the casks vertically in a staggered manner, such as by adding an offset of a few centimetres in one of the rows, as is the case with 3 × 10, Z_h = -160 and 3 × 10, Z_h = -170 where the middle row was raised by 10 and 20 cm respectively, resulted in an increase in k_eff as the casks are separated by D cm. The rise in k_eff is a result of H in Figure 8, having to refocus itself—both horizontally and vertically—to maintain its focus at L/2 all the time. Because of that, H increases in both directions (vertically and horizontally), the minimum solid angle decreases as it is inversely proportional to H (refer to Eqn. 11). The k_eff will follow accordingly based on its correlation with the solid angle (refer to Eqn. 10).

However, when an analysis of the behaviour of k_eff with respect to an increase in Z_h is performed, with shorter intervals but over longer offset, it was discovered that the k_eff fluctuated with an increase in offset, and the most economically optimum offset (i.e. offset that results in the lowest k_eff) is Z_h = 16.75 cm (refer to Figure 10 and Table 19). Based on this, it is recommended that an offset be taken into account in cask storage, but only after a cask storage matrix analysis has been performed and an optimum offset has been found, which will result in the highest decrease in k_eff.

For a given burnup regime, the sensitivity of the k_eff to the change in separation distance (D) between spent fuel casks is a function of the magnitude of the area over which the fissile material is dispersed. As burnup increases, the concentration of fission products also increases, which results in a decrease in fission density because of neutrons required for fission being absorbed. The change in D affects the geometric buckling $B_g^2$, which increases as a result of increase in H as D increases. An analysis has to be performed to determine whether the system is critical, subcritical, or supercritical.

• If $B_g^2=B_m^2$, the system is stable and self-sustaining,
• However, if $B_g^2<B_m^2$, the system is critical with a net increase in neutron production, and the k_eff increases.
• If $B_g^2>B_m^2$, the neutron leakage is too high, which is the case with 2 × 15 Zigzag storage matrix.

Given that there is a linear relationship between burnup and flux, selection of a spent fuel storage pattern must take into consideration the burnup, because when the burnup increases, there will be more cross-linking/interference between the flux from adjacent casks, which will reduce resonance escape probability and decrease the k_eff.

Thus, the calculation of the relationship between k_eff and the distance between adjacent casks and, hence, the selection of an appropriate storage pattern, should be performed prior to the placement of spent fuel casks within the spent fuel storage facility, with a view to determining which storage matrix will optimise storage space. This (i.e., analysis of storage matrix) must be submitted to and approved by the nuclear regulatory authority to prevent k_eff exceeding the regulatory limit. This is very important, as some storage matrices often lead to an increase in k_eff as the distance between casks increases.

Furthermore, because of the interaction and coupling of neutron fluxes from different casks, it is advisable to store casks with different burnups separately to prevent thermal neutron flux from underburned fuel assemblies from coupling with that from other casks, particularly if the spent fuel casks have different designs and are designed for different enrichment levels. This interaction may alter the orientation and the relationship between $B_g^2$ and $B_m^2$, potentially resulting in a self-sustaining system, at worst, a critical system.

Author Contributions

The author did all the research work for this study.

Competing Interests

The authors have declared that no competing interests exist.

Data Availability Statement

The dataset generated and/or analysed during the current study are included in published article and its supplementary information files.

AI-Assisted Technologies Statement

ChatGPT was employed to improve the readability and linguistic clarity of the English text. All scientific content, data interpretation, and conclusions were developed independently by the author. The authors have thoroughly reviewed and edited the AI-assisted text to ensure its accuracy and accept full responsibility for the content of the manuscript.

Additional Materials

The following additional materials are uploaded at the page of this paper.

1. Appendix 1: Trad_Matrix_100.inp.