## Abstract

Nuclear science and engineering is a field increasingly dominated by computational studies resulting from increasingly powerful computational tools. As a result, analytical studies, which previously pioneered nuclear engineering, are increasingly viewed as secondary or unnecessary. However, analytical solutions to reduced-fidelity models can provide important information concerning the underlying physics of a problem and aid in guiding computational studies. Similarly, there is increased interest in sensitivity analysis studies. These studies commonly use computational tools. However, providing a complementary sensitivity study of relevant analytical models can lead to a deeper analysis of a problem. This work provides the analytical sensitivity analysis of the one-dimensional (1D) cylindrical mono-energetic neutron diffusion equation using the forward sensitivity analysis procedure (FSAP) developed by Cacuci. Further, these results are applied to a reduced-fidelity model of a spent nuclear fuel cask, demonstrating how computational analysis might be improved with a complementary analytic sensitivity analysis.

## 1 Introduction

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],

“Given a differential equation to analyze, most people spend only a small amount of time using analytical tools and then use a computer to see what the solution ‘looks like.’”

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],

“… for a long time… [these] solutions were treated by most researchers as though they were merely isolated ‘exact’ solutions to special problems: elegant, sometimes useful, but extremely limited in significance. It was only gradually realized that these solutions were actually of much broader significance…”

and moreover, by Zaitsev and Polyanin [17],

“… exact solutions of differential equations play an important role in the proper understanding of qualitative features of many phenomena and processes in various areas of natural science…”

and finally, by Sachdev and Janna [18],

“… the search for exact solutions is now motivated by the desire to understand the mathematical structure of the solutions, and hence, a deeper understanding of the physical phenomena described by them. Analysis, computation, and not insignificantly, intuition all pave the way to their discovery…”

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],

“… understanding the validity and place of exact/approximate analytical solution(s) in the general context can be greatly enhanced by numerical simulation. In short, there must be a continuous interplay of analysis and computation if a… problem is to be successfully tackled.”

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.

### 1.1 Sensitivity Analysis.

As characterized by Saltelli et al. [19],

“… sensitivity analysis is aimed… at priority setting, to determine what factor most needs better determination, and to identify the weak links of the assessment chain (those that propagate most variance in the output)…”

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:

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.

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.

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,

“… if instead of adding fuel we move a small amount of the existing fuel from one place to another, there will be an increase in reactivity if the fuel is moved from a region of low [importance] to a region of high [importance]: [importance] is a statistical weight for fuel…”

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 *δ*:

“… there are several forms of the statistical weight. The appropriate weighting for neutrons in a system is $\u2020N$. When there is a perturbation in properties, $\delta [\eta \Sigma fv]$ for example, the source or rate of production of new neutrons is $\delta [\eta \Sigma fv]N$. Then the weighting of the perturbation as an effective source is by $\u2020NN$ where

Ncauses the source and $\u2020N$ describes its effect. Suitably normalized, such a bilinear weighting is known as the statistical weight of the region…”

where, in this discussion, *N* and $\u2020N$ are the forward and adjoint neutron populations or fluxes, respectively, $\eta \Sigma fv$ is essentially a fission neutron multiplication factor, and $\eta \Sigma 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.

### 1.2 Motivation and Charter.

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:

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.

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.

## 2 Mathematical Model

The fuel region of the HI-STORM 100 spent fuel cask features various materials including spent UO_{2} 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 $\u226b$ 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.

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

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\u2261$ position vector within *V* (cm),

$\varphi (r)\u2261$ neutron scalar flux at point $r(cm\u22122\u2009s\u22121)$,

$J(r)\u2261$ net neutron current at $r(cm\u22122\u2009s\u22121)$,

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

$\Sigma a(r)\u2261$ probability per unit path length at point $r$ that a neutron is absorbed; otherwise known as the macroscopic total absorption cross section (cm^{−1}),

$\Sigma f(r)\u2261$ probability per unit path length at point $r$ that a neutron induces a fission event; otherwise known as the macroscopic fission cross section (cm^{−1}),

$\nu \xaf\u2261$ average number of neutrons born per fission event,

$S(r)\u2261$ neutron source density (e.g., spontaneous fission) at point $r(cm\u22123\u2009s\u22121)$,

the salient reaction rates within *V* are defined by

$\u222d\Sigma a(r)\varphi (r)dV=$ number of neutrons absorbed per unit time in $V\u2009(s\u22121)$,

$\u222d\nu \xaf\Sigma f(r)\varphi (r)dV=$ number of neutrons born via fission per unit time in $V\u2009(s\u22121)$,

$\u222fJ(r)\xb7nds=$ number of neutrons per unit time that that leak from *V* through the surface $s\u2009(s\u22121)$,

$\u222dS(r)dV=$ number of neutrons born from source per unit time in $V\u2009(s\u22121)$.

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 $\varphi $ 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 $\varphi $ and $J$.

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

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.

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.

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.

*or S*

_{n}*expansions of the far more general Boltzmann neutron transport equation (see, for example, Duderstadt and Hamilton or Lewis and Miller [39,40]), from which it also follows:*

_{n}*and Σ*

_{t}*are, respectively, the macroscopic total interaction and isotropic scattering cross sections at point $r$, such that*

_{s}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/\Sigma 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].

*s*of

*V*and

*d*is an “extrapolation distance” (suggested via notation to be uniformly added to each element of $rb$ in Eq. (9)) given by

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.

*r*=

*0, as*

## 3 Local Sensitivity Analysis Primer

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 $\varphi $ values), $\alpha $ 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 $\u2202R/\u2202\alpha $. However, this approach can be algebraically involved since it requires solving the set of equations $F$ for each parameter variation.

*M*number of parameters as

*ϵ*is interpreted as an infinitesimal deviation from the nominal value of a given parameter, and the rightmost expression is the definition of the G-derivative. In general, the evaluated result of Eq. (20) can be written for each element of $\delta R$ as

*j*th element of $\delta R$,

*η*contains sensitivity information for the parameter

_{i}*α*, and

_{i}*M*is the number of parameters being studied. The values of

*η*are used to calculate the sought-after sensitivity coefficients, which provide a relative comparison between parameters. The sensitivity coefficients are thus calculated using $\delta Rj$ as

_{i}where $Sc,\alpha i$ is the sensitivity coefficient for parameter *α _{i}* and

*R*is the response value corresponding to the

_{j}*j*th element of $\delta 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.

## 4 Local Sensitivity Analysis of Representative Spent Fuel Cask 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*^{0} is the intrinsic neutron source, *I*_{0} is the modified Bessel function of the first kind, and $r\u03030$ 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.

*j*subscript is dropped as there is only one response function in this demonstration. Finally, determining the sensitivities for each input parameter using Eqs. (24)–(28) in Eq. (20) is equivalent to replacing each input parameter in Eq. (23) with

*r*-dependent functions appearing in Eq. (31) are defined by

*B*

^{0}as well as in

*D*

^{0}. In practice, the values for

*D*

^{0},

*B*

^{0}, and $r\u03030$ are calculated from experimental data or geometry (in the case of $r\u03030$). Therefore, it is necessary to express each of the above input parameters according to their individual definitions using the following equations:

where $\Sigma 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 $\Sigma a0\u2261\Sigma c0+\Sigma f0$. Fundamental sensitivity coefficient results written in terms of the parameters Σ* _{s}*, Σ

*, $\nu \xaf$, Σ*

_{c}*, and*

_{f}*r*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).

_{b}Redefining the sensitivity coefficients for *B*, *D*, and $r\u0303$ in terms of those for Σ* _{c}*, Σ

*, $\nu \xaf$, Σ*

_{s}*, and*

_{f}*r*is a straightforward process similar to how the coefficients were found for

_{b}*B*,

*D*, and $r\u0303$ above. Taking the G-derivative of each of Eqs. (40)–(42), each equation is redefined to be expressible in the terms $\delta \Sigma c,\u2009\delta \Sigma s,\u2009\delta \nu \xaf,\u2009\delta \Sigma f$, and $\delta rb$. These definitions are then used in the sensitivity coefficients summarized in Eqs. (36)–(39) to yield the final expressions.

*δB*,

*δD*, and $\delta r\u0303$ as

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

### 4.1 Numerical Example.

Numerical evaluation of the analytical sensitivity coefficients calculated within Sec. 4 proceeds only under prescription of the nominal input parameters *S*^{0}, $\Sigma c0,\u2009\Sigma s0,\u2009\nu \xaf0,\u2009\Sigma f0$, and $rb0$ (and thus *D*^{0}, *B*^{0}, and $r\u03030$). 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*^{0}, $\Sigma c0,\u2009\Sigma s0,\u2009\nu \xaf0,\u2009\Sigma f0$, and $rb0$ values that are themselves somehow representative of the detailed cask model.

In turn, as *S*^{0}, $\Sigma c0,\u2009\Sigma s0,\u2009\nu \xaf0,\u2009\Sigma 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*^{0}, $\Sigma c0,\u2009\Sigma s0,\u2009\nu \xaf0$, and $\Sigma 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.

Parameter | Value |
---|---|

S^{0} | 20.14 neutrons/cm^{3} s |

$\Sigma c0$ | 0.061 1/cm |

$\Sigma f0$ | 0.0090 1/cm |

$\Sigma s0$ | 0.10 1/cm |

$\nu \xaf0$ | 2.65 neutrons |

r_{b} | 74.68 cm |

Parameter | Value |
---|---|

S^{0} | 20.14 neutrons/cm^{3} s |

$\Sigma c0$ | 0.061 1/cm |

$\Sigma f0$ | 0.0090 1/cm |

$\Sigma s0$ | 0.10 1/cm |

$\nu \xaf0$ | 2.65 neutrons |

r_{b} | 74.68 cm |

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*^{0}. 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}*, Σ

*, $\nu \xaf$, Σ*

_{s}*, and*

_{f}*r*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:

_{b}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).The sensitivity coefficient associated with

*r*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_{b}*r*is effectively perturbing the location at which the boundary condition given Eq. (9) is applied; and for this reason, $Sc,rb$ increases drastically from_{b}*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.The sensitivity coefficient associated with the capture cross section Σ

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,\Sigma 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)._{c}Figure 9 shows that positive perturbations in $\nu \xaf$ 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 Σ

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,\nu \xaf$ and $Sc,\Sigma f$ as appearing in Fig. 9, the two coefficients are not identical since Σ_{f}appears decoupled from $\nu \xaf$ as part of its inclusion in the definition of_{f}*D*given by Eq. (23).Otherwise, the sensitivity coefficients associated with Σ

, $\nu \xaf$, Σ_{f}, and Σ_{s}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_{c}*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.The sensitivity coefficient associated with the scattering cross section Σ

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,\Sigma s$ would then be uniformly positive, which is clearly a nonphysical result in the neighborhood of the cask outer boundary._{s}

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:

For a majority of the cask radius, Σ

is the most important input parameter; however, its importance drops near the cask outer radius as a result of the increase in $Sc,rb$._{c}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$.For the inner portion of the cask radius, $\nu \xaf$ and Σ

are the third and fourth most sensitive parameters, respectively. However, the sharp increase in $Sc,rb$ relegates $\nu \xaf$ and Σ_{f}to the fourth and fifth most important parameters near_{f}*r*= 40 and 30 cm for $\nu \xaf$ and Σ, respectively._{f}For approximately the first 30 cm of the cask,

*r*is the fifth most important parameter until approximately_{b}*r*= 30 cm, where $S,c,rb$ increases and overtakes all the other parameters to become the most important parameter in the system.For a majority of the cask radius, Σ

is the least important input parameter; however, it becomes the fifth most important parameter very near the cask outer radius._{s}

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,\Sigma s$ and $Sc,\Sigma c$ as depicted in Figs. 9 and 10, for example, at the point where $Sc,\Sigma s$ changes sign, $Sc,\Sigma c$ changes slope.

### 4.2 Comparison With Computational Model.

where $\sigma x,0$ is the unperturbed cross section value, $\varphi (\sigma x,0)$ is the neutron flux evaluated with respect to the nominal value cross sections, and $\Delta \sigma x\u2261\sigma x\u2212\sigma x,0$.

*p*, as

_{x}*p*is used to normalize the amount of response change to the amount of perturbation. Further, Eq. (54) can be rewritten in terms of

_{x}*p*after applying the chain rule as

_{x}*σ*is defined as

_{x}$\varphi 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}*, Σ

*, and Σ*

_{s}*) from the helium MCNP model are displayed.*

_{c}^{2}

Figure 11 shows that the analytic and computational values of $Sc,\Sigma 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,\Sigma 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,\Sigma 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,\Sigma f$ do not decrease as a result.

Comparing the analytically calculated values of $Sc,\Sigma 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,\Sigma 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,\Sigma 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,\Sigma c$ initially show a reduction around *r *=* *60 cm, near the location where $Sc,\Sigma s$ goes negative (*r *=* *66.84 cm). From this, a relationship is again seen between the loss terms, Σ* _{c}* and Σ

*. The computationally computed values of $Sc,\Sigma 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,\Sigma c$ in the computational model. Another difference between the two models occurs near the boundary of the fuel region, at*

_{s}*r*=

*73 cm where the analytically calculated values of $Sc,\Sigma 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,\Sigma c$ and $Sc,\Sigma 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 Σ

*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 Σ*

_{s}*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.*

_{f}## 5 Conclusion

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:

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

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.

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

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.

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.

### 5.1 Recommendations for Future Work.

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

_{n}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.

## Acknowledgment

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.

## Funding Data

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

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

## Nomenclature

*B*=material buckling parameter

*d*=extrapolated distance

*D*=diffusion coefficient

- $e\u20090$ =
concatenation of nominal state vector and input parameter values

- $F$ =
model equations

- $hy$ =
vector of perturbed state values

- $h\alpha $ =
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,\alpha i$ =
sensitivity coefficient pertaining to parameter

*α*_{i}- V =
control volume

- $y$ =
state vector

- $\alpha $ =
vector of the system input parameters

- $\alpha i0$ =
nominal value of parameter

*α*_{i}- $\delta \alpha i$ =
perturbation of parameter

*α*_{i}*ϵ*=infinitesimal deviation

*η*=_{i}sensitivity parameter

- $\nu \xaf$ =
number of neutrons emitted per fission event

- Σ
=_{a} neutron absorption interaction cross section ($\Sigma a=\Sigma c+\Sigma f$)

- Σ
=_{c} neutron capture cross section

- Σ
=_{f} fission cross section

- Σ
=_{s} scattering interaction cross section

- Σ
=_{t} total interaction cross section

- $\varphi (r)$ =
neutron flux, also the total track length traveled per unit time per unit volume by neutrons through a control volume

## Footnotes

The sensitivity coefficients for *S* and *r _{b}* are not directly computable from MCNP using the PERT card. MCNP results are given in units of

*per source neutron*, therefore the value of $Sc,S$ is likely 1 as there is a linear relationship between the internal source strength and MCNP simulated neutron flux. However, discussion of the computational sensitivity coefficients is limited to values that are entirely computationally attainable, and $Sc,S$ is not. Determining the values of $Sc,rb$ through computational means would require running multiple simulations where

*r*, the fuel radius, is changed in each simulation, as geometry perturbations are precluded from MCNP's perturbation capabilities. Finally, perturbations in $\nu \xaf$ are not compatible with MCNP perturbation capabilities, and therefore, cannot be computationally determined.

_{b}### Appendix A: Detailed Model Summary

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.% ^{235}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 $(\alpha ,n)$ reaction and spontaneous fission sources, respectively.

Isotope | Source strength neutrons/cm^{3} s | Weight fraction |
---|---|---|

$234U$ | 5.307 × 10^{−5} | 1.087 × 10^{−3} |

$238Pu$ | 1.743 × 10^{−1} | 8.338 × 10^{−5} |

$239Pu$ | 2.512 × 10^{−2} | 0.004 |

$240Pu$ | 4.072 × 10^{−2} | 0.002 |

$241Pu$ | 1.222 × 10^{−4} | 0.001 |

$242Pu$ | 1.201 × 10^{−4} | 3.829 × 10^{−3} |

$241Am$ | 1.797 × 10^{−1} | 2.081 × 10^{−5} |

$243Am$ | 1.400 × 10^{−3} | 6.823 × 10^{−5} |

$242Cm$ | 1.671 × 10^{−7} | 7.585 × 10^{−6} |

$243Cm$ | 7.315 × 10^{−4} | 1.28 × 10^{−7} |

$244Cm$ | 1.350 × 10^{−1} | 1.738 × 10^{−5} |

Isotope | Source strength neutrons/cm^{3} s | Weight fraction |
---|---|---|

$234U$ | 5.307 × 10^{−5} | 1.087 × 10^{−3} |

$238Pu$ | 1.743 × 10^{−1} | 8.338 × 10^{−5} |

$239Pu$ | 2.512 × 10^{−2} | 0.004 |

$240Pu$ | 4.072 × 10^{−2} | 0.002 |

$241Pu$ | 1.222 × 10^{−4} | 0.001 |

$242Pu$ | 1.201 × 10^{−4} | 3.829 × 10^{−3} |

$241Am$ | 1.797 × 10^{−1} | 2.081 × 10^{−5} |

$243Am$ | 1.400 × 10^{−3} | 6.823 × 10^{−5} |

$242Cm$ | 1.671 × 10^{−7} | 7.585 × 10^{−6} |

$243Cm$ | 7.315 × 10^{−4} | 1.28 × 10^{−7} |

$244Cm$ | 1.350 × 10^{−1} | 1.738 × 10^{−5} |

Isotope | Source strength neutrons/cm^{3} s | Weight fraction |
---|---|---|

$233U$ | 1.623 × 10^{−13} | 1.791 × 10^{−9} |

$234U$ | 1.201 × 10^{−7} | 1.087 × 10^{−3} |

$235U$ | 1.123 × 10^{−8} | 0.007 |

$236U$ | 2.040 × 10^{−6} | 0.003 |

$238U$ | 1.687 × 10^{−3} | 0.819 |

$237Np$ | 5.239 × 10^{−9} | 2.961 × 10^{−3} |

$238Pu$ | 3.326 × 10^{−2} | 8.338 × 10^{−5} |

$239Pu$ | 9.692 × 10^{−6} | 0.004 |

$240Pu$ | 2.985 × 10^{−1} | 0.002 |

$241Pu$ | 1.882 × 10^{−7} | 0.001 |

$242Pu$ | 1.005 × 10^{−1} | 3.829 × 10^{−3} |

$241Am$ | 8.223 × 10^{−5} | 2.081 × 10^{−5} |

$243Am$ | 6.827 × 10^{−6} | 6.823 × 10^{−5} |

$242Cm$ | 8.743 × 10^{−7} | 7.585 × 10^{−6} |

$243Cm$ | 3.961 × 10^{−6} | 1.281 × 10^{−7} |

$244Cm$ | 1.906 × 10^{1} | 1.738 × 10^{−5} |

$245Cm$ | 1.430 × 10^{5} | 8.515 × 10^{−7} |

$246Cm$ | 9.711 × 10^{2} | 6.809 × 10^{−8} |

Isotope | Source strength neutrons/cm^{3} s | Weight fraction |
---|---|---|

$233U$ | 1.623 × 10^{−13} | 1.791 × 10^{−9} |

$234U$ | 1.201 × 10^{−7} | 1.087 × 10^{−3} |

$235U$ | 1.123 × 10^{−8} | 0.007 |

$236U$ | 2.040 × 10^{−6} | 0.003 |

$238U$ | 1.687 × 10^{−3} | 0.819 |

$237Np$ | 5.239 × 10^{−9} | 2.961 × 10^{−3} |

$238Pu$ | 3.326 × 10^{−2} | 8.338 × 10^{−5} |

$239Pu$ | 9.692 × 10^{−6} | 0.004 |

$240Pu$ | 2.985 × 10^{−1} | 0.002 |

$241Pu$ | 1.882 × 10^{−7} | 0.001 |

$242Pu$ | 1.005 × 10^{−1} | 3.829 × 10^{−3} |

$241Am$ | 8.223 × 10^{−5} | 2.081 × 10^{−5} |

$243Am$ | 6.827 × 10^{−6} | 6.823 × 10^{−5} |

$242Cm$ | 8.743 × 10^{−7} | 7.585 × 10^{−6} |

$243Cm$ | 3.961 × 10^{−6} | 1.281 × 10^{−7} |

$244Cm$ | 1.906 × 10^{1} | 1.738 × 10^{−5} |

$245Cm$ | 1.430 × 10^{5} | 8.515 × 10^{−7} |

$246Cm$ | 9.711 × 10^{2} | 6.809 × 10^{−8} |

Monte Carlo N-particle requires input of a neutron source energy spectrum. The source spectrum is found using the Oak Ridge National Laboratory SCALE code named ORIGEN-S [50]. ORIGEN-S is an irradiation and decay code used, in this work, to calculate the energy spectrum of a material based on intrinsic neutron sources in the material. Figure 15 shows the neutron energy spectrum found using ORIGEN-S, which is caused by spontaneous fission of nuclei in the fuel (such as ^{252}Cf) and (*α*, *n*) reactions generating neutrons in the fuel.

Figure 16 shows the scalar neutron flux integrated over the height of the HI-STORM 100 spent fuel cask as a function of radial position. The fuel region is the innermost region, followed by the MPC, the air region, the concrete annulus, and finally the carbon steel shell. The vertical lines denote the interface between each of the material regions. The neutron flux is reduced by 54% in the fuel region and there is a further reduction of 39% occurring in the concrete annulus as shown in Fig. 16. These reductions in the neutron flux occur as a result of multiple processes: (1) the comparatively dense materials in the fuel region, (2) the inclusion of neutron-absorbing materials (e.g., boron), and (3) the large amount of thermalizing isotopes (e.g., hydrogen) present in the concrete annulus.

### Appendix B: Helium Model Summary

The helium model is a hypothetical model developed to investigate the effects of reducing geometric details in the fuel region and for comparison of sensitivity analysis results with the analytic model. Remedes et al. also feature this model for use in the results assessment methodology [36]. This MCNP model homogenizes the materials inside of the 32 fuel cells by weight. The homogenized weight fraction of each isotope is found by taking the total weight of each isotope in the 32 fuel cells divided by the total mass of all isotopes present in the 32 fuel cells. The homogenized fuel material is represented as a cylinder, shown in the top-view image (Fig. 17) and the side-view image (Fig. 18). The volume of the cylinder is determined to match the volume of the original 32 fuel cells. The region exterior to the homogenized fuel cylinder and interior to the MPC is filled with helium.

Figure 19 shows the neutron flux through the fuel region of the helium model as compared to the neutron flux simulated using the detailed model. The inset graph shows the relative error between the neutron flux simulated in the detailed model and the helium model. Remedes et al. provide an in-depth analysis of the results shown in Fig. 19 [36].