Analytical Sensitivity Analysis of a Spent Nuclear Fuel Cask

Like its counterparts in physics, engineering, biology, earth science, population dynamics, economics, and many other fields, nuclear engineering is a discipline dominated by the exercise of computational horsepower, and the perpetual search for agreement between the results of those computations and experiments. For example, flagship computational physics simulation software packages (hereafter referred to as “codes”) have by now long served as the tools underpinning the design and evaluation of nuclear reactors, various types of detectors, radiation shields, radiation sources, the characterization and assessment of various radiation fields, and multitudinous other applications.

The extensive use and thoroughly proven effectiveness of codes for this incredibly diverse set of applications have in a sense relegated to a secondary role the analytical studies pioneered in the early days of nuclear engineering by luminaries such as Fermi [1], Goodman [2], Weinberg and Wigner [3], and Case and Zweifel [4]. These venerable efforts and associated analytical solutions of entities such as neutron transport equations are now most often relegated to use as educational tools in the classroom setting but are otherwise sometimes dismissed as not having substantial relevance or practical application to real-world scenarios.

Highlighting the well-known limitations of mathematical models (e.g., differential equations) featuring analytical solutions is almost always honestly conceived, in that any such model is not likely to be applicable outside of a highly specialized set of idealized circumstances. Indeed, as noted by Zwillinger [5],

Such is one of the principal drivers in the genesis and widespread use of computational science codes. Analytical solutions to differential equations are rarely if ever attainable in generalized circumstances, and so codes are often designed with the explicit intention of rectifying this deficiency within some predefined scope of applications but at the cost of producing approximate solutions.

Despite this now firmly established paradigm, in the past 25 years existing and new analytical solutions within a variety of nuclear science and engineering contexts have enjoyed something of a renaissance as “test problems” for “benchmarking” the very codes that once replaced them. This newer use of analytical modeling has proceeded in tandem with the rigorous codification of integrated verification, validation, and uncertainty quantification practices throughout the broader computational science community. These practices are increasingly viewed as essential for establishing both the pedigree and predictive capability of any computational science code.

The notion of code benchmarking is essentially encoded within the broader notion of quantitative code verification, as discussed thoroughly by Oberkampf et al. [6], Roache [7], and Oberkampf and Roy [8]. In short, one outcome of code verification is enhanced assurance that a code solves its underlying equations in a manner consistent with the expectations of the code developer. As noted by Oberkampf et al. [6], code verification proceeds “by comparing computational solutions with highly accurate solutions,” including but not limited to “analytic solutions for simplified physics.” With this philosophy in mind, several compendia of code verification test problems for nuclear engineering applications have been developed. For example, in the context of neutron transport and its surrogates, collections of analytical and semi-analytical solutions intended at least in part for application as test problems owe to Ganapol [9], Myers [10], Sood et al. [11], and most recently Tsyfra and Czyzycki [12], Ramsey et al. [13], Giron et al. [14], and Pittman et al. [15].

Looking beyond even their immense utility in quantitative code verification studies, analytical results within simplified mathematical models are valuable for a variety of other purposes. From a practical standpoint, known functional forms (as opposed to numerical tables) are optimal for both exhaustive insight into essential phenomenology and broader application, especially if those forms are differentiable or integrable. To this point, as noted by Barenblatt [16],

and moreover, by Zaitsev and Polyanin [17],

and finally, by Sachdev and Janna [18],

As such, even in an age dominated by computational studies, there remains a distinct role for the development and implementation of analytical treatments. As summarized by Sachdev and Janna [18],

Now-ubiquitous quantitative code verification efforts may thus be viewed as a core element within a broader program of study; this integrated picture should also feature complementary analytical studies to aid in the understanding of essential phenomenologies underpinning increasingly complicated computational science simulations.

As characterized by Saltelli et al. [19],

Sensitivity analysis is a prime example of a field that has in recent years been almost wholly subsumed by purely computational endeavors; the literature is voluminous surrounding site-specific efforts and implementations within physics, engineering, biology, earth science, population dynamics, economics, and many other areas. This state of affairs is also reflected in the expositions of both the canonical primers and various critiques surrounding the subject.

The evolution of sensitivity analysis from its historical and largely analytical roots to modern-day computational programs of study has largely proceeded in tandem with parallel developments in large-scale computational science, as described in Sec. 1. The reasons for this outcome include but are not necessarily limited to:

1. Sensitivity analysis has emerged as one of the cornerstone processes through which generalized mathematical models or codes may be assessed through an integrated program of verification, validation, and uncertainty quantification. Modern-day sensitivity studies are therefore most commonly encountered in the context of code evaluations.

2. Fully coupled, global sensitivity analyses (as opposed to their “one factor at a time” or purely local or derivative-based counterparts) are increasingly viewed as necessary components of fully rigorous code assessment strategies. These techniques typically demand repeated code execution under coupled sampling spanning the entire space of possible parameter realizations.

3. Tremendous advances in the power and widespread availability of high-performance computing resources has made conceptually and computationally simple “brute-force method” sensitivity analysis approaches more viable than ever before. In these studies, even large or otherwise complicated codes can be rapidly, repeatedly executed under conjoined parameter sampling.

In these often-necessary realizations of sensitivity analysis practices, complementary analytical studies assume the same historical relevance and underpinning basis, limitations, and advantages as outlined in Sec. 1 for the attendant mathematical models. A classical example of this phenomenology with broad relevance to the nuclear engineering community is adjoint analysis, as detailed among many others by Keepin [20], Henry [21], and Lewins [22]. As noted by Lewins [22] in the context of reactor fuel loading and Goertzel's theorem,

and, moreover, Lewins [22] also establishes the essential equivalence between the notions of statistical weight, importance and the adjoint function, and sensitivity according to its interpretation as a local perturbation δ:

where, in this discussion, N and $†N$ are the forward and adjoint neutron populations or fluxes, respectively, $ηΣfv$ is essentially a fission neutron multiplication factor, and $ηΣfvN$ is a fission neutron production rate. As such, the deep analytical connectivity between the concepts of importance, statistical weight, and sensitivity is well established, and in the modern day widely used in computational sensitivity analysis settings [23].

Like the scenarios discussed in Sec. 1, Stacey, Greenspan, and Lewins and Becker discuss the appropriateness of analytical adjoint based approaches in nuclear engineering, since adjoint sensitivity calculations are best suited for applications when one is interested in the change in response to many input parameters [24–26]. Nuclear engineering is not unique in that a desired response (e.g., the neutron flux) depends on many input parameters (e.g., nuclear data). However, if an analytical model has a closed-form solution, then solving the system from the forward direction is more appropriate since the process is considerably more straightforward and does not require knowledge of the adjoint solution.

Closely related to analytical adjoint methods is Cacuci's “forward sensitivity analysis procedure” (FSAP) based on the concept of the Gâteaux generalized directional derivative [27]. This formalism (if not its practical implementation) is entirely analytical and often features minimal computational overhead; in turn, however, it only provides decoupled sensitivity information on the local or first-order level. As such, when implemented in the context of analytical mathematical models, the FSAP shares the same general drawbacks and benefits associated with analytical modeling in general, and as outlined in Sec. 1.

In short, the potentially narrow scope of analytically computed sensitivity information is counterbalanced by a more complete and lucid functional representation that may prove informative of more general developments. This important idea forms the motivation of this work, to be initiated in a specific context with relevance to the nuclear engineering community.

The notions set forth in Secs. 1 and 1.1 are entirely general and may be applied in any number of contexts. Of particular relevance to this study is the implementation of the aforementioned techniques in the context of modeling and simulation of containers intended for the long-term storage of spent nuclear reactor fuel (hereafter referred to as “spent fuel casks”).

Indeed, over the past decade modeling and analysis of spent fuel casks have proceeded along almost entirely computational lines. Computational modeling efforts pertaining to spent fuel casks are extensive; these studies are motivated by a lack of a long-term nuclear waste storage site for spent power reactor fuel. As a result, many studies have investigated models of spent fuel casks using computer codes, includes those conducted by Harkness et al. [28], Miller et al. [29], Liao and Yang [30], and Greulich et al. [31] who used various codes to ensure the continued viability of spent fuel casks.

In tandem with increasing computational evaluations of spent fuel casks, some small amount of verification, validation, and uncertainty quantification has also been performed to corroborate and rigorize various simulation results. Sensitivity analysis pertaining to spent fuel casks has been mostly focused on transportation issues. Sanders and Westfall [32], Radulescu et al. [33], and Brady and Sanders [34] investigate how cask transportation limits are sensitive to spent fuel material compositions. However, as sensitivity analysis in this context is not limited to cask transportation limits, Gao conducted a computational sensitivity analysis study investigating the effects of multigroup and continuous energy cross section data and changes to geometric fidelity [35].

Against this broader backdrop, and in the spirit of the ideas presented in Secs. 1 and 1.1, complementary analytical sensitivity analysis modeling of spent fuel cask scenarios is intended to serve two essential purposes:

1. These studies can serve as guides for targeting future, entirely computational sensitivity analysis studies pertaining to spent fuel cask scenarios, thus potentially achieving computational cost savings where necessary.

2. As described throughout Sec. 1, analytical studies may also serve as guides for interpreting, understanding, and rigorizing certain results of existing and future computational studies pertaining to spent fuel casks.

The primary objective of this work is therefore to execute an analytical study along these lines. The target for analysis is a Holtec International HI-STORM 100 spent fuel cask as built off previous work described by Remedes et al. [36]. Figure 1 shows a Holtec HI-STORM 100 spent nuclear fuel cask with a partially inserted multipurpose canister (MPC) [37]. Spent nuclear fuel is held in the MPC where each rectangular prism or “fuel cell” holds 264 thin cylindrical fuel rods and 25 instrumentation rods. The overpack (the region consisting of concrete and outer shell) provides passive radiation attenuation and a simple neutron flux distribution. However, the region interior to the MPC, hereafter referred to as the “fuel region,” features a complex neutron flux distribution which suggests a need for deeper analysis. Further, the fuel region is the source of radiation in the system, making it the most important region for investigation and requiring the most analysis.

In support of this objective, Sec. 2 contains additional information on the computational details surrounding a representative HI-STORM 100 fuel region, and the associated derivation of a representative analytical mathematical model. Section 3 provides a discussion of the general concepts of sensitivity analysis as well as a detailed description of the analytical methods used in this work. Section 4 motivates the use of a reduced-fidelity model by comparing details from the high-fidelity and reduced-fidelity models. After corroborating the reduced-fidelity model, this section analyzes the results of this sensitivity investigation on the analytic models and on the HI-STORM 100 spent fuel cask. It is important to note that the analytic model used in Sec. 4 is chosen to be simple and easy to understand to demonstrate the process as well as the types of conclusions which can be found by undertaking an analysis of this kind. Finally, Sec. 5 reviews some conclusions and recommendations for future study.

The fuel region of the HI-STORM 100 spent fuel cask features various materials including spent UO₂ nuclear fuel, a stainless-steel basket, boron-containing neutron absorbing pads, and helium backfill. The geometric configuration of these materials is highly complex, as depicted in Fig. 1. As such, a single mathematical model capable of describing the neutron flux in the fuel region is not analytically tractable; therefore, it must be treated computationally. Then, informed by the results of the detailed computational study, a simplified representative model is developed using a variety of assumptions and approximations. In essence, analyzing properties of the detailed model (i.e., neutron energies and directions) leads to information that aids in determining an analytic model for sensitivity analysis investigations.

A high-fidelity computational model of the spent fuel cask appearing in Fig. 1 has been developed using the Monte Carlo N-particle (MCNP) code system, as detailed in Appendix  A. In order to identify appropriate modeling simplifications to inform an analytic model choice, the energy spectrum and angular distribution of the neutron flux (as simulated using the high-fidelity MCNP model) and associated cross section data for various materials are analyzed at various locations in the fuel region.

Figure 2 shows the energy spectrum of the neutron flux throughout the cask fuel region as calculated by MCNP; these plots show that the neutron energy spectrum has little variation throughout the fuel region. This outcome results from the uniform distribution of fuel rods through the fuel region. Further, the lack of thermalizing materials (which are used to lower the energy of neutrons) in the fuel region suggests that there is little change in the neutron spectrum. Therefore, it is assumed that energy dependence of the neutrons is uniform throughout the fuel region (for the purpose of analytic modeling). Further analysis of Fig. 2 also indicates how to analytically treat the neutron flux energy dependence throughout the fuel region. According to Figs. 2(a)–2(h), the percentage of neutrons above 10 keV varies between 78% at inner radius values to 71% at the edge of the fuel region. As such, a mono-energetic handling of the neutron flux energy dependence in the analytic model is assumed since the majority of neutrons have energies between 10 keV and 10 MeV, if an appropriate group weighting spectrum is used as described by Bell and Galsstone [38]. These results are also useful to corroborate the reduced-fidelity MCNP model developed in Sec. 4.

After choosing a method for analytically treating the energy dependence of the neutron flux, it is also necessary to determine a method for handling its directional dependence. Figure 3 shows the angular distribution, as calculated by MCNP, of the neutron flux 0.5 cm from the centerline (Fig. 3(a)) and at the edge of the fuel region (Fig. 3(b)); the angular distribution is tallied at these locations to capture the two extents of the angular neutron flux. Figure 3(a) shows the neutron flux is slightly inward-peaked 0.5 cm from the centerline, with 50.28% of all neutrons at that position traveling toward the centerline. Figure 3(a) further indicates the neutron flux at the outer edge of the fuel region is outward-peaked, with 57.29% of all neutrons at that position traveling away from the centerline. This phenomenon owes to significant neutron leakage from the fuel region to exterior regions featuring smaller neutron population densities. All together these results suggest that the neutron flux is nearly isotropic throughout the cask fuel region, with a maximum deviation from isotropy of about 7% near its outer edge.

Finally, Fig. 4 shows the neutron mean free path (MFP; or average distance between interactions) associated with each material present within the fuel region, as calculated using MCNP. The MFP within the helium components is approximately 1 km for neutron energies near 1 MeV, where the neutron source flux is most intense. As the thickest region of helium appearing in the fuel region is a 10 cm run between the fuel cells and edge of the fuel region, the overwhelmingly large MFP indicates that there will be a negligible number of neutrons interacting in helium. As such, the first material assumption is that the helium outside of the fuel cells is negligible for the purpose of constructing a counterpart analytic model.

The remaining fuel region materials feature neutron MFPs of approximately 1 cm for neutron energies near 1 MeV. The various thicknesses of these materials as appearing in the fuel region are also of the order of 1 cm; as such, their propensity to interact with neutrons cannot be neglected. However, since these materials are evenly distributed (i.e., the materials exist throughout the fuel region and not just at a single location) and feature similar MFPs, a homogenization technique is used to approximate the geometry in the fuel region. Consistent with this interpretation, a simplified geometric representation of the cask fuel region features a cylinder-shaped homogeneous fuel material preserving the weight ratio of each material in the fuel region. The volume of the homogenous cylinder of fuel material is determined to preserve the volume from the 32 original fuel cells; the resulting radius of the cylinder is approximately 75 cm. The volume around the cylinder of homogenous fuel is treated as a vacuum in the mathematical model. The radius of the homogenized fuel is about two orders of magnitude greater than the MFP of the materials used in the fuel region (e.g., 100 cm radius of fuel $≫$ 1 cm MFP).

In summary and in the interest of constructing a complementary analytical representation of the neutron population behavior within the cask fuel region (as calculated using MCNP), inspection of the various features appearing in Figs. 1–4 suggests several modeling simplifications.

1. Aside from various isotope production and depletion processes featuring characteristic time scales spanning weeks to years, the spent fuel cask is essentially a static object. It is therefore assumed that the analytical representation of the cask is entirely time-independent (hereafter referred to as “static”).

2. The neutron energy spectrum within the fuel region is uniform and essentially “fast;” that is, it principally exists at fission neutron energies (i.e., 1–2 MeV) with minimal thermalization. As such, the analytical model used to characterize the cask fuel region is taken be approximately monoenergetic.

As a consequence of these observations and associated simplifications, a static, mono-energetic balance law model is used to analytically characterize the neutron population information within the cask fuel region.

To derive this model, for a convex but otherwise arbitrary control volume V (with infinitesimal element $dV$⁠, surface s, and infinitesimal surface element $ds$ with an associated unit normal $n$⁠), and given the definitions

$r≡$ position vector within V (cm),

$ϕ(r)≡$ neutron scalar flux at point $r(cm−2 s−1)$⁠,

$J(r)≡$ net neutron current at $r(cm−2 s−1)$⁠,

and the most probable neutron reactions occurring within the spent fuel cask,

$Σa(r)≡$ probability per unit path length at point $r$ that a neutron is absorbed; otherwise known as the macroscopic total absorption cross section (cm⁻¹),

$Σf(r)≡$ probability per unit path length at point $r$ that a neutron induces a fission event; otherwise known as the macroscopic fission cross section (cm⁻¹),

$ν¯≡$ average number of neutrons born per fission event,

$S(r)≡$ neutron source density (e.g., spontaneous fission) at point $r(cm−3 s−1)$⁠,

the salient reaction rates within V are defined by

$∭Σa(r)ϕ(r)dV=$ number of neutrons absorbed per unit time in $V (s−1)$⁠,

$∭ν¯Σf(r)ϕ(r)dV=$ number of neutrons born via fission per unit time in $V (s−1)$⁠,

$∯J(r)·nds=$ number of neutrons per unit time that that leak from V through the surface $s (s−1)$⁠,

$∭S(r)dV=$ number of neutrons born from source per unit time in $V (s−1)$⁠.

As written and within the scope of its underlying assumptions, Eq. (4) is an exact neutron balance relationship describing the interaction between the neutron scalar flux $ϕ$ and net neutron current vector $J$⁠, both of which are regarded as unknowns. As such, the first-order differential equation given by Eq. (4) cannot be solved without the invocation of an additional closure relation between $ϕ$ and $J$⁠.

At this point, additional observations and assumptions resulting from inspection of Figs. 1–4 may be used to further simplify Eq. (4):

1. The neutron angular distribution within the cask fuel region is essentially “flat;” that is, there is no significant indication of neutron motion in some preferred direction. Consequently, both neutron scattering and neutron transport within the fuel region is taken to be isotropic.

2. The neutron mean free path at most locations within the fuel region is generally much smaller than the distance at that same location to the outer edge of the fuel region itself. This result implies that neutrons typically undergo numerous interactions within the fuel region prior to escape.

3. The angular distribution and mean free path trends in turn indicate that from a macroscopic standpoint, neutrons act diffusively throughout the majority of the fuel region. As such, an equivalent homogeneous fuel region is assumed. Moreover, the mass-averaged material properties within this homogeneous fuel region are assumed to be uniform.

where the constant D is called the “diffusion coefficient.” Equation (5), also known as Fick's Law, is introduced under the observation that neutrons act diffusively throughout the majority of the cask fuel region. Indeed, Eq. (5) is indicative of a diffusion process: it indicates that neutrons move “down the gradient” from regions of high neutron concentration to low.

The diffusion approximation has known limitations. Duderstadt and Hamilton [39] outline the conditions that must be met in order for the diffusion approximation to be valid as: (1) the neutron flux is several MFPs (⁠$1/Σt$⁠) from material boundaries or neutron sources, (2) the material is weakly absorbing, and (3) the neutron current is changing slowly as compared to time between interactions [39].

Equation (9) is intended to qualitatively reproduce the neutron flux behavior at the outer surface of a nonreentrant convex body, as otherwise observed from more general neutron transport scenarios (see, for example, Lamarsh and Baratta [41]). Since Eq. (8) is second-order, it requires one more boundary condition which will be introduced in due course.

Further simplification of Eq. (8) is afforded by the specification of a spatial coordinate system, so as to resolve the Laplacian operator appearing therein. Again, additional observations and assumptions resulting from Fig. 1 may be used to further simplify Eq. (8):

• The spent fuel cask is a right circular cylinder with a height much larger than the outer radius of the interior fuel region. From the perspective of a central slice along its height, the spent fuel cask may therefore be regarded as approximately infinite in the vertical direction. This observation in turn suggests a one-dimensional (1D) cylindrical representation of the fuel region.

With an acute understanding of the assumptions that give rise to them, the k = 0 and 1 instantiations of Eq. (11) (with boundary conditions given by Eq. (9)—with $rb→rb$—and Eq. (13)) will be used for the analytical sensitivity analysis studies conducted throughout the remainder of this work.

where $R$ is a vector of the system responses, $F$ is a vector of the model equations (e.g., a vector containing the diffusion equation), $y$ is the state vector (e.g., a vector of $ϕ$ values), $α$ is the vector of the system input parameters, where the vector $F$ can also represent nonlinear model equations; however, the following discussion is limited to linear equations for the purpose of this work.

Re-expressing Eq. (15), using Eq. (16), yields the sought after sensitivity information $∂R/∂α$⁠. However, this approach can be algebraically involved since it requires solving the set of equations $F$ for each parameter variation.

where $Sc,αi$ is the sensitivity coefficient for parameter α_i and R_j is the response value corresponding to the jth element of $δR$ [45].

The sensitivity coefficients are used to determine the “importance” of each parameter. Parameters with larger sensitivity coefficients have a larger impact on the system response. The signs of sensitivity coefficients are also relevant, as they indicate the direction of change in the response given a change in a parameter. Meaning, if the sensitivity coefficient has a negative value for a given parameter, increasing the value of that parameter will cause the value of the response to decrease. On the other hand, if the sensitivity coefficient has a positive value, increasing the associated parameter value will cause an increase in the response value.

Using this methodology, Sec. 4 investigates the effects of perturbing input parameters relating to nuclear data in the solution to Eq. (8), the 1D cylindrical mono-energetic neutron diffusion equation. The intention here is to demonstrate the utility of this procedure (which could be applied to even more complicated models) by executing it on a simple, well-known model.

Section 2 introduced the 1D cylindrical mono-energetic neutron diffusion equation. This model uses experimental data in the form of cross sections. These cross sections are input parameters, and this work demonstrates sensitivity analysis on the aforementioned diffusion equation with respect to them.

where S⁰ is the intrinsic neutron source, I₀ is the modified Bessel function of the first kind, and $r̃0$ is the extrapolated radius of the fuel region equivalent to $rb0+d0$⁠. Again, the superscript 0 is used to denote the nominal value of each input parameter or response function.

where $Σc0$ is the nominal capture cross section, $rb0$ is the outer radius of the nominal cask fuel region, and the nominal total absorption cross section is redefined using $Σa0≡Σc0+Σf0$⁠. Fundamental sensitivity coefficient results written in terms of the parameters Σ_s, Σ_c, $ν¯$⁠, Σ_f, and r_b are then determined by applying the G-derivative to each of Eqs. (40)–(42) and substituting the results into their respective places in Eqs. (32)–(35).

Redefining the sensitivity coefficients for B, D, and $r̃$ in terms of those for Σ_c, Σ_s, $ν¯$⁠, Σ_f, and r_b is a straightforward process similar to how the coefficients were found for B, D, and $r̃$ above. Taking the G-derivative of each of Eqs. (40)–(42), each equation is redefined to be expressible in the terms $δΣc, δΣs, δν¯, δΣf$⁠, and $δrb$⁠. These definitions are then used in the sensitivity coefficients summarized in Eqs. (36)–(39) to yield the final expressions.

These values are then substituted into Eq. (31).

Numerical evaluation of the analytical sensitivity coefficients calculated within Sec. 4 proceeds only under prescription of the nominal input parameters S⁰, $Σc0, Σs0, ν¯0, Σf0$⁠, and $rb0$ (and thus D⁰, B⁰, and $r̃0$⁠). However, central to the analytical developments in Secs. 2 and 4 is the definition of a homogeneous cask fuel material featuring these lumped input parameters. The relevance of the surrogate analytical results featured in Secs. 2 and 4 to any detailed (i.e., nonhomogeneous) cask fuel model thus hinges on the development of a set of S⁰, $Σc0, Σs0, ν¯0, Σf0$⁠, and $rb0$ values that are themselves somehow representative of the detailed cask model.

In turn, as S⁰, $Σc0, Σs0, ν¯0, Σf0$⁠, and $rb0$ are associated with a homogeneous fuel representation, their relevance to a high-fidelity cask representation (e.g., the detailed model as described in Sec. 2) may be assessed through the development and simulation of a complementary MCNP model featuring an equivalent homogeneous fuel region. Figure 5 depicts the top–down view of such a hypothetical cask as developed in MCNP; this instantiation preserves isotopic weight fractions within the fuel region as featured within the detailed model. Put another way, all materials contained in the 32 fuel bundles present in the detailed model are homogenized preserving isotopic mass fraction. Details of this model are further discussed in Appendix  B. For ease of discussion, this MCNP model will henceforth be referred to as the “helium model” due to the annulus of helium surrounding the otherwise homogenous fuel.

Comparison of MCNP simulation results between the detailed and helium models indicates, in turn, the potential appropriateness and limitations of employing surrogate analytical models featuring homogenized fuel (i.e., such as those employed throughout Secs. 2 and 4) for informing more detailed simulation studies. In particular, credibility of the homogeneous fuel assumption and the favorable analytical consequences thereof is associated with, for example, and in part, quantifiable “agreement” between the two computational models in the calculated neutron energy spectrum and angular distribution throughout the cask fuel region.

To this point, Fig. 6 compares the neutron energy spectrum at various locations in the fuel region between the detailed and helium models. The helium model consistently under-predicts the thermal component of the neutron energy spectrum, as seen in Figs. 6(a)–6(h). The helium model under-predicts the neutron flux at lower energies because the model features evenly distributed material throughout the entire fuel region (extending to 74.680 cm). The key materials are the evenly distributed fuel and neutron absorbing pad materials which readily absorb thermal neutrons. Since materials are evenly distributed in the helium model, rather than having discrete extents surrounded by free-streaming regions (i.e., helium backfill within a fuel pin lattice and structural materials, as in the detailed model), comparatively more thermal neutrons are absorbed on fuel and neutron absorbing pad materials. The final subfigure, Fig. 6(h), also shows an increase in neutron energy spectrum near 0.02 eV for both models, due to thermalized neutrons returning to the fuel region from the surrounding concrete.

Figure 7 depicts the neutron angular distribution within each computational model at 1.25 cm from the centerline (Fig. 7(a)) and at 84.30 cm from the centerline (Fig. 7(b)). In the helium model, 50.03% of the neutrons are traveling away from the cask centerline (outward) at 1.25 cm, compared to 49.71% of the neutrons traveling outward at the same location within the detailed model. In the helium model, 57.10% of the neutrons are traveling outward at 84.34 cm, compared to 57.28% of the neutrons traveling outward at the same location within the detailed model.

Together, Figs. 6 and 7 indicate that within the indicated metrics, the helium model is representative of the detailed model at approximately the 85% level (that is, the indicated metrics differ by maximum of 15% between the two models). Therefore, cross section and other input parameter data derived from the helium model can be expected to be representative of the detailed model to no more than the same degree of fidelity. Indeed, a variety of other simulation output metrics are available for quantitative comparison of the helium and detailed computational models (e.g., the neutron current evaluated throughout the fuel region), and the attendant assessment of the homogeneous fuel assumption underpinning the developments of Secs. 2 and 4. Moreover, the helium model itself is not unique among the possible computational surrogate models potentially associated with the analytical outcomes featured in Secs. 2 and 4. Accordingly, the comparisons provided in Figs. 6 and 7 are intended to serve principally, for example, purposes. Moreover, and in any event, in justifying the selection of surrogate analytical models, representative input parameters, and nominal input parameter values, it is most imperative to identify sources of error as opposed to outright elimination of each source of error. In turn, identification and quantification of any such discrepancies illuminates the strengths and weakness of the various included models, which then aids in adequately characterizing their underling physics and the overall appropriateness for a scenario of interest.

With these notions in mind, the representative homogeneous fuel composition employed in the helium model may therefore be used to determine an associated set of nominal input parameters S⁰, $Σc0, Σs0, ν¯0$⁠, and $Σf0$ for use with the analytical results appearing in Sec. 4, featuring an associated quantification of their relevance to the detailed model. Given fuel composition of the helium model, these input parameters are evaluated using the nuclear data processing code NJOY [46], where the necessary calculations are proceeded by weighting the cross section values against the neutron source energy spectrum. Otherwise, the nominal input parameter $rb0$ is the radius of the homogenized fuel material, 74.68 cm. Table 1 provides a summary of parameter values calculated for the homogeneous fuel associated with the helium model.

Figure 8 compares the scalar neutron flux as calculated using both the MCNP detailed model and the hypothetical homogeneous fuel material calculated using Eq. (23) with Table 1 data. The neutron flux decreases approximately by half over the fuel region, and for the analytical result, the shape of the curve is determined by the input parameter B⁰. At 74.68 cm, the analytical neutron flux becomes flat since the neutrons have entered the free-streaming region at a comparatively large cylindrical radius (i.e., the infinite region of vacuum surrounding the fuel). The neutron flux calculated from the detailed model features a similar response at the same location, since at that point, neutrons have entered the largely noninteracting helium backfill region surrounding the fuel cells. The inset graph in Fig. 8 shows the error between the scalar neutron flux predicted by the two models; in this metric, the two models again agree within 15% across the entire cask fuel region. Some of the principal details causing this variation include three small depressions present in the detailed model results; the stainless-steel lattice and neutron absorbing pads cause these depressions, the fine structure of which is effectively “smoothed out” in homogeneous fuel models.

Figure 9 depicts the sensitivity coefficients $Sc,i$ associated with the elemental parameters $i=$S, Σ_c, Σ_s, $ν¯$⁠, Σ_f, and r_b appearing within the analytical model given by Eq. (23), as calculated using Eqs. (32)–(35), (36)–(39), and (50)–(52) and the data appearing in Table 1. Several trends are immediately evident from Fig. 9:

1. The sensitivity coefficient associated with the intrinsic neutron source term S is unity since the source term itself appears simply as a scalar multiplier within Eq. (23).

2. The sensitivity coefficient associated with r_b exhibits the most dramatic change across the radius of the cask fuel region. In fact, the value of this sensitivity coefficient increases to 13.24 at 74.68 cm. Perturbing r_b is effectively perturbing the location at which the boundary condition given Eq. (9) is applied; and for this reason, $Sc,rb$ increases drastically from r = 40 cm to r = 74.78 cm since boundary conditions are imperative in constructing unique analytical solutions. This phenomenon also explains why $Sc,rb<0.04$ for the first 40 cm of the cask r = fuel region, as the flux at these values is less affected by the boundary condition at 74.68 cm, and more affected by the boundary condition at the centerline given by Eq. (13). Finally, $Sc,rb>0$ since increasing the outer radius value forces the flux to remain at higher values through the radius of the fuel.

3. The sensitivity coefficient associated with the capture cross section Σ_c is negative throughout the entire homogenous fuel region. This phenomenon indicates that as the capture cross section increases, the neutron flux decreases. This behavior is physically plausible since capture is a pure loss mechanism (i.e., as more neutrons are lost to capture, the value of the neutron flux becomes smaller). $Sc,Σc$ has an inflection point and increases in value near 73 cm from the centerline, since loss terms are forcing the flux to meet to the boundary condition given by Eq. (9).

4. Figure 9 shows that positive perturbations in $ν¯$ cause uniformly positive perturbations in the neutron flux. This trend is physically plausible since increasing the number of neutrons generated through fission events will increase the flux value throughout a multiplying material. Along these same lines, the sensitivity coefficient for the fission cross section Σ_f is also uniformly positive since increasing the likelihood of fission will in turn increase the number of neutrons in the homogeneous fuel material (i.e., as the number of neutrons available for transport increases, the flux increases). Moreover, while there appears to be a strong correlation between $Sc,ν¯$ and $Sc,Σf$ as appearing in Fig. 9, the two coefficients are not identical since Σ_f appears decoupled from $ν¯$ as part of its inclusion in the definition of D given by Eq. (23).

5. Otherwise, the sensitivity coefficients associated with Σ_f, $ν¯$⁠, Σ_s, and Σ_c have a similar shape: they are nearly flat for a majority of the cask's radial extent, before trending toward zero near the outer surface of the cask. This phenomenon is a consequence of all these terms appearing within the definition of B as given by Eq. (23), which in turn controls the shape of the analytical neutron flux. The relationship between these input parameters demonstrates how the structure of the neutron flux is related to the structure of the sensitivity coefficients, since the G-derivative is a linear operator.

6. The sensitivity coefficient associated with the scattering cross section Σ_s exhibits the most nontrivial behavior; it is positive and increasing for $r<$ 66.84 cm, positive and decreasing for 66.84 cm $<r<70.93$ cm, and negative for $r>$ 70.93 cm to the cask outer radius. In turn, these features are indicative of the relative importance of a variety of gain and loss mechanisms occurring within Eq. (23). In particular, for $r<$ 70.93 cm neutron scattering serves a gain mechanism: it acts to spatially redistribute but otherwise preserve the neutron flux within the monoenergetic diffusion model (i.e., in the absence of thermalization). For $r>$ 70.93 cm, neutron scattering is a loss mechanism: scattering in proximity to the outer boundary of the fuel region serves to increase leakage processes where neutrons are lost through the outer surface of the fuel material. The inflection point occurring at $r=$ 66.84 cm is then indicative of the spatial location where the role of neutron scattering begins to transition: its presence owes to the approximate nonreentrant boundary condition given by Eq. (9), which is intended to include leakage mechanics within the analytical diffusion model. That is, if the neutron flux was instead terminated at the physical extent of the fuel region, the analytical model would predict no neutron leakage and rather a zero neutron flux there. In this case, $Sc,Σs$ would then be uniformly positive, which is clearly a nonphysical result in the neighborhood of the cask outer boundary.

To further understand and better rank the importance of the various competing physical phenomenologies included in Eq. (23), Fig. 10 depicts the absolute value of each sensitivity coefficient plotted in Fig. 9. Several additional trends are immediately evident from Fig. 10:

1. For a majority of the cask radius, Σ_c is the most important input parameter; however, its importance drops near the cask outer radius as a result of the increase in $Sc,rb$⁠.

2. For a majority of the cask radius, S is the second most important input parameter; however, near r = 50 cm, $Sc,S$ is briefly the most important parameter before becoming the second most important parameter resulting from the increase in $Sc,rb$⁠.

3. For the inner portion of the cask radius, $ν¯$ and Σ_f are the third and fourth most sensitive parameters, respectively. However, the sharp increase in $Sc,rb$ relegates $ν¯$ and Σ_f to the fourth and fifth most important parameters near r = 40 and 30 cm for $ν¯$ and Σ_f, respectively.

4. For approximately the first 30 cm of the cask, r_b is the fifth most important parameter until approximately r = 30 cm, where $S,c,rb$ increases and overtakes all the other parameters to become the most important parameter in the system.

5. For a majority of the cask radius, Σ_s is the least important input parameter; however, it becomes the fifth most important parameter very near the cask outer radius.

The trends manifesting in Figs. 9 and 10 are principally due to the r-dependent interplay between the capture and leakage loss mechanisms present in Eq. (23). For example, capture is the dominant loss mechanism near the cask centerline, as shown in Fig. 9: a neutron initially located there is most likely to undergo many interactions before escaping from the cask outer surface. Conversely, leakage becomes an increasingly important loss mechanism near the cask outer radius, the importance of which is observed to eventually exceed that of capture. This physical interplay noticeably manifests in the behavior of $Sc,Σs$ and $Sc,Σc$ as depicted in Figs. 9 and 10, for example, at the point where $Sc,Σs$ changes sign, $Sc,Σc$ changes slope.

where $σx,0$ is the unperturbed cross section value, $ϕ(σx,0)$ is the neutron flux evaluated with respect to the nominal value cross sections, and $Δσx≡σx−σx,0$⁠.

$ϕ0$ is the unperturbed neutron flux. The values of the sensitivity coefficients calculated by Eq. (57) are then comparable to the definition provided in Eq. (22).

Figure 11 shows the comparison between the sensitivity coefficients calculated from the helium model in MCNP and the analytic sensitivity coefficients. The sensitivity coefficients from the detailed model are not included in the analysis of the fuel region, as limitations of MCNP's perturbation capabilities preclude sensitivity analysis in this region. Moreover, only the analytic sensitivity coefficients which have comparable computational values (Σ_f, Σ_s, and Σ_c) from the helium MCNP model are displayed.²

Figure 11 shows that the analytic and computational values of $Sc,Σf$ agree within 5% throughout the fuel region. Near r = 60 cm, the value of the analytic sensitivity coefficients pertaining to the fission cross section begins to decrease where the computational values of $Sc,Σf$ remain flat. Here is the first example of the effect of boundary values. The analytic model in the fuel region requires a boundary value at r = 78.79 cm (the extrapolated radius of the fuel). The flux is chosen to vanish at the extrapolated boundary value. Therefore, the value of the analytic $Sc,Σf$ decreases reflecting the decreasing flux value approaching the boundary value location. The computational model does not share this boundary value, and the computational values of $Sc,Σf$ do not decrease as a result.

Comparing the analytically calculated values of $Sc,Σs$ and their computational counterparts in Fig. 11 shows considerable disagreement. There are clear benefits when the two models agree; in this scenario, a code user understands the physics at the level of the analytic models. However, the scenario when the two models do not agree still provides insight into the problem leading to a more rigorous analysis of a simulation. The relative error between the two models is nearly constant over the first 50 cm of the fuel region for each of the parameters. The analytic values of $Sc,Σs$ are positive due to the choice of a mono-energetic analytic model. That is, choosing a mono-energetic model prevents the neutrons from slowing down which, in turn, does not capture the how the probability of absorption increases as neutrons thermalize. As a result, neutron scattering can only act as a loss term through leakage, which will not occur until a neutron is significantly close to a boundary. The computational models use continuous energy cross section data, which captures thermalization and indirectly leads to neutron loss through capture of thermal neutrons, in addition to the aforementioned leakage process, yielding a negative sensitivity coefficient value. Further, the two models have similar shapes, flat before breaking downward. This behavior occurs in the analytic models because of the boundary value at the extrapolated radius. The sensitivity coefficients pertaining to loss terms increase in value at the boundary value since the neutron flux is being forced to zero. However, the computational model does not have boundary values at this location. Instead, the neutron flux is decreasing because the flux is directed outward (from Fig. 3(b)) confirming that neutrons are leaking from the fuel region. Therefore, increasing the scattering cross section will increase the chance that a neutron leaks from the fuel region. Even though the two models do not agree, understanding the causes for the disagreement is as important in understanding the problem as having matching results.

Finally, the values of $Sc,Σc$ for analytic and computational models disagree as shown in Fig. 11. However, both curves are negative, since absorption is a loss term. The values of the analytically determined $Sc,Σc$ initially show a reduction around r = 60 cm, near the location where $Sc,Σs$ goes negative (r = 66.84 cm). From this, a relationship is again seen between the loss terms, Σ_c and Σ_s. The computationally computed values of $Sc,Σc$ do not show a drastic reduction at the same location, since scattering is always a loss term in the computational models. That is, there is no location where scattering physics changes from a preservation term to a loss term; therefore, there is no drastic change in $Sc,Σc$ in the computational model. Another difference between the two models occurs near the boundary of the fuel region, at r = 73 cm where the analytically calculated values of $Sc,Σc$ has an inflection point. This inflection point results from the boundary value forcing the flux to zero, resulting in larger negative values for $Sc,Σc$ and $Sc,Σs$⁠. The computational model does not have this inflection point since the flux is not forced to zero at this location in the helium model.

The inconsistency between the sensitivity coefficients from the helium model and analytic models is a result of oversimplifying the continuous energy cross section data when using only two energy groups. In an effort to investigate the effect of better representing the continuous cross section data by increasing the number of energy groups, a 30-group instantiation of the helium model is developed in MCNP, using the preloaded multigroup formulations included in that code [47]. While the cross section data is energy dependent in the 30-group MCNP models, the simulated neutron flux values are integrated over energy, similarly to the neutron flux values of the continuous energy MCNP models. Figure 11 shows the results from the 30-group sensitivity analysis. The sensitivity coefficients from the 30-group model for Σ_c and Σ_s better agree with the corresponding sensitivity coefficients from the helium model than the sensitivity coefficients from the analytic and detailed models. However, the sensitivity coefficients pertaining to Σ_f from the 30-group model show more disagreement than from the analytic model. This conclusion shows that 30 energy groups are insufficient to capture the sensitivity information from the detailed model in the fuel region. However, the 30-group sensitivity coefficients allude to the disagreement between the analytic and computational sensitivity coefficients being a result of using too few energy groups in the analytic model.

Analytical or otherwise reduced-fidelity models are useful tools for enhancing analysis of high-fidelity physics simulations based on the extensive computational modeling used in nuclear engineering. Furthermore, performing sensitivity analysis on these same models reveals the underlying mathematical structure inherent to a scenario, leading to a deeper understanding of the salient physics. Incorporating a study of appropriate analytical models into this paradigm acts as part of a broader program of study, which underpins the results from increasingly complicated computational science simulations. Further, the addition of analytically computed sensitivity information proves informative as a guide in interpreting, understanding, and rigorizing results of existing and future computational studies.

In the spirit of established analytical and computational model comparison techniques and outcomes as described in Sec. 1, the various analytical results, examples, and commentary provided in Secs. 2, 4, and 4.1 represent an example of how an analytical sensitivity analysis study can be used to setup, precondition, and eventually inform or compare against a complementary computational sensitivity analysis study. Within this conceptual strategy, and against the backdrop of the detailed MCNP computational model of a HI-STORM 100 spent nuclear fuel storage cask, the results appearing within Secs. 2, 4, and 4.1 exemplify a more general recipe justifying the development and execution of local sensitivity analysis formalisms within the context of surrogate analytical models:

1. Establish a high-fidelity computational model and extract key features of the simulation output.

2. Based on these key features, establish a reduced-fidelity model of the same underlying scenario; preferably this model is amenable to analytical or semi-analytical solution.

3. Execute a sensitivity analysis study on the reduced-fidelity model; again, preferably this study will be amenable to analytical or semi-analytical evaluation.

4. Scenario dependent evaluation of the analytical or semi-analytical sensitivity structure requires nominal input parameters; these must also be consistent with the key features extracted from the high-fidelity computational model.

5. Establish scenario-dependent sensitivity trends and input parameter importance ranking to precondition additional high-fidelity computational sensitivity analysis studies.

Sections 2, 4, and 4.1 exemplify this process in its application to a 1D cylindrical mono-energetic neutron diffusion equation model of a spent nuclear fuel cask fuel region. For example, in this case, the parasitic capture cross section and 1D cylindrical outer radius of the cask fuel region were found to be the parameters causing the most uncertainty in the neutron flux interior to the cask fuel region. Moreover, the inconsistencies present between the analytically computed sensitivity coefficients for Σ_s and their computationally calculated counterparts illuminate the importance of handling the energy dependence, as demonstrated in Sec. 4.2. Further, this shows that discrepancies between results provide value to the analysis process by identifying salient physics.

More broadly, results of this type are capable of guiding future research to reduce uncertainty in the most important input parameters inherent to a given scenario of interest. Further, by identifying the most important input parameters, a code user can identify if any simplifications were made when developing a simulation input which would affect the results. From these conclusions, a user could either change the simulation input to address any insufficiencies or explain the insufficiencies and identify pathways for improvement. Either decision results in a more thorough examination of the scenario of interest, which is ultimately the goal of any scientific study.

The purpose of the numerical example present in this work is to demonstrate the process of determining and comparing sensitivity coefficients between analytic and computational models. However, the concepts presented in this work extend beyond the models provided herein. In addition to this necessary program of study, there appears to be a nearly limitless sequence of higher-fidelity analytical fuel cask models in which the G-derivative formalism may be brought to bear. Candidate analytical models along these lines include but are not necessarily limited to multigroup neutron diffusion models, multigroup P_n or S_n neutron transport models, and multigroup integral or integro-differential neutron transport models. Depending on the physical processes of interest, each of these models may be formulated as static or time-dependent, in various representative geometries, and featuring any number of multimaterial regions. Again, the ultimate intent of analytical sensitivity analysis studies within any of these formalisms is to enable comparison to complementary computational results, using the formalism exemplified in Secs. 2 and 4.

Finally, and as indicated in Sec. 1.2, programs of sensitivity analysis as applied to computational models of spent nuclear fuel casks appears to be an area ripe for further advancement in research and development. This being the case, and in tandem with the aforementioned potential for new, analogous analytical treatments, there also appears to be ample opportunity for the computational evaluation of not only local sensitivity information as pertaining to spent fuel casks, but also the more complete global metrics as described by Saltelli et al. and many other authors.

Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy (Contract No. 89233218CNA000001).

The authors would also like to thank Joe Schmidt, Jeff Favorite, Tom Saller, C. J. Solomon and Joel Kulesza for their mentorship and knowledge in using the computational tools.

• US Department of Energy through the Los Alamos National Laboratory (Funder ID: 10.13039/100000015).

• National Nuclear Security Administration (Funder ID: 10.13039/100006168).

B =

material buckling parameter

d =

extrapolated distance

D =

diffusion coefficient

$e 0$ =

concatenation of nominal state vector and input parameter values

$F$ =

model equations

$hy$ =

vector of perturbed state values

$hα$ =

vector of perturbed input parameters values

$J(r)$ =

net neutron current at $r$

k =

geometric parameter where k = 0, 1, 2 for planar, cylindrical, and spherical geometries, respectively

M =

number of parameters in sensitivity study

MCNP =

Monte Carlo N-particle program

MFP =

mean free path

MPC =

multipurpose canister

$r$ =

position vector in cm

$R$ =

system responses

$rb$ =

radius vector to material boundary

S =

neutron source density

$Sc,αi$ =

sensitivity coefficient pertaining to parameter α_i

V =

control volume

$y$ =

state vector

$α$ =

vector of the system input parameters

$αi0$ =

nominal value of parameter α_i

$δαi$ =

perturbation of parameter α_i

ϵ =

infinitesimal deviation

η_i =

sensitivity parameter

$ν¯$ =

number of neutrons emitted per fission event

Σ_a =

neutron absorption interaction cross section (⁠$Σa=Σc+Σf$⁠)

Σ_c =

neutron capture cross section

Σ_f =

fission cross section

Σ_s =

scattering interaction cross section

Σ_t =

total interaction cross section

$ϕ(r)$ =

neutron flux, also the total track length traveled per unit time per unit volume by neutrons through a control volume

The foundational basis for this work is an MCNP simulation of the HI-STORM 100 spent fuel cask. Remedes et al. also use the results of this MCNP model as a foundational element for the results assessment methodology, and a description of the model is provided therein [36]. A brief description of what is called the “detailed model” is also provided here for completeness. This MCNP model is designed to account for many of the small geometric details of the cask. The purpose of this appendix is to briefly describe the development of the MCNP model.

The MPC and overpack are modeled in MCNP to predict the simulated interior neutron flux spatial distribution averaged over the height of the cask as a function of radial distance from the centerline. Figure 12 shows the side view of the cask, and the corresponding cross-sectional view is provided in Fig. 13. In Fig. 12, the innermost rectangles are the interior lattice of the MPC which holds the fuel cell. The rectangular lattice structure is easier to see in Fig. 13 as the square lattice of boxes that bound each fuel cell. The materials inside the MPC are surrounded by helium. The MPC itself encloses the fuel materials, steel lattice, and helium gas. The MPC is surrounded by air that can enter the cask through vents in the concrete (the thickest annulus). Figure 12 shows the locations of four of eight of the vents as the two gaps at the bottom of the concrete and two gaps at the top of the concrete. Air is used to transfer decay heat away from the MPC and into the environment. The MPC sits on a carbon steel plate and a concrete pedestal. Above the MPC is another carbon steel plate and concrete. The entire outer radius, as well as top and bottom, of the concrete is covered with a layer of carbon steel to protect the cask from the environment as well as from accidental damage. Each square in Fig. 13 shows the locations of a fuel cell. A single fuel cell is discussed next.

Figure 14 is a zoomed in view of a single fuel cell in the fuel region of the detailed model. The smaller circles are fuel rods, which are clad in zircalloy. The larger circles are water regions used to represent instrumentation. Each fuel cell contains two boron-carbide and aluminum pads which are called the “neutron absorbing pads” in this work. Finally, the square lattice is made of stainless steel 304 structures which hold the fuel cells in place. The fuel material is characterized using data from the Next Generation Safeguards Initiative spent fuel libraries [49]. These libraries contain data representative of used nuclear fuel rods from various fuel bundle types ranging in initial enrichment values and burn-up times. The fuel for this work is chosen to come from a Westinghouse 17 × 17 fuel bundle with an initial enrichment of 3 wt.% ²³⁵U and a burn-up value of 30 GWd/MTU. The probabilistic nature of fission means each fuel rod will have a unique isotopic composition which would induce small variations into the interior neutron flux and potentially hide features occurring in the neutron flux. Therefore, an average fuel rod composition is developed based on the total mass of each isotope present in all fuel rods of a single fuel cell. Handling the fuel composition in this manner makes salient physics more readily available for identification. Tables 2 and 3 summarize the isotope name, neutron source strength, and weight fraction of the $(α,n)$ reaction and spontaneous fission sources, respectively.