Development of a Coupled Depletion Perturbation Theory Methodology in Continuous-Energy Monte Carlo Depletion Simulations

Abstract

This work develops a depletion perturbation theory methodology coupled with generalized perturbation theory to produce fully coupled nuclide number density sensitivity coefficients for cross sections, initial nuclide number densities, and decay constants. The method is implemented into OpenMC, a high-fidelity Monte Carlo simulation code, utilizing continuous-energy cross sections to produced high resolution sensitivity coefficients over the energy spectrum.

Keywords:

• Depletion perturbation theory
• sensitivity coefficients
• OpenMC

I. INTRODUCTION

Depletion simulations allow Monte Carlo simulations to predict the time-dependent evolution of a system’s nuclide inventories. Significant uncertainty exists in nuclear data for short-lived or rare isotopes (as illustrated in Fig. 1) in fission product yields and in radioisotope decay branching ratios, yet few methods exist to resolve the impact of uncertainty in nuclear data in time-dependent problems. Furthermore, analyzing the impact of this uncertainty using stochastic uncertainty quantification methods can be computationally challenging, as Monte Carlo transport–driven depletion simulations are generally computationally intensive, even for simple cases. Schemes such as the predictor-corrector method require three Monte Carlo simulations (once for the predictor in the initial system state, once for the corrector, and once for the final system state). As a result, uncertainty quantification methods for depletion simulations are limited in resolution and inherently computationally expensive.

Fig. 1. Uncertainty in nuclear data for short-lived radioisotopes limits the accuracy of depletion simulations. ENDFB/VII.0 cross sections for ²³⁸Np (red) lack the nuclear data resolution present in data for ²³⁷Np.

Figure 1 demonstrates an example of the measurement limitations of nuclear data. A half-life of 2.117 days for ²³⁸Np makes it difficult to isolate the nuclide and accurately measure the cross section. It may be possible to improve the nuclear data for ²³⁸Np by observing the sensitivity of longer-lived or stable isotopes to ²³⁸Np and applying data assimilation methods to infer the ²³⁸Np values that maximize agreement between the simulation results and the experimental measurements for systems that contain ²³⁸Np. While this work does not achieve such lofty goals, the potential to infer nuclear data that cannot be otherwise measured directly has motivated the development of these depletion sensitivity analysis methods. This work combines depletion perturbation theory (DPT)^[1] with recent advances in reaction rate sensitivity methods^[2] to enable high-resolution, adjoint-based sensitivity analysis for continuous-energy Monte Carlo transport–based depletion simulations.

Generalized perturbation theory (GPT) methods estimate the sensitivity of reaction rates in the Bateman depletion matrix, which are coupled with adjoint functionals from these depletion equations to obtain depletion sensitivity coefficients.^[2] This methodology was initially implemented in a simple one-dimensional (1-D) test code^[3] and is demonstrated here for a simple pressurized water reactor (PWR) pin cell model using the OpenMC Monte Carlo transport code. Sensitivity coefficients are estimated for nuclide densities with respect to cross-section data, initial nuclide densities, and radioactive decay constants and are then evaluated in accuracy through comparison with reference direct perturbation sensitivity coefficients.

II. THEORY

Sensitivity coefficients are relative derivatives that quantify the fractional change in a response $\mathit{R}$ with respect to a fractional change in data $\mathit{\alpha }$,

(1) $S_{\mathit{R},\mathit{\alpha }}=\frac{\mathrm{\partial R}/R}{\mathrm{\partial \alpha }/\alpha }.$(1)

This work assumes $\mathit{R}$ to be the number density of nuclides and $\mathit{\alpha }$ is either a neutron cross section for a specific nuclide, the initial number density of a specific nuclide, or a decay constant for a specific nuclide.

The theory for basic depletion sensitivity coefficients was developed by Williams in 1979.^[1] The depletion sensitivity coefficient can be derived from first-order perturbations to the forward and backward depletion and buildup equations.

The forward depletion equation can be written as

(2) $\frac{d\vec{n}}{\mathit{dt}}=\mathit{A}\vec{n},$(2)

where $\mathit{A}$ is a transmutation matrix describing the pathways for nuclide transmutation. For the specific production of a nuclide $n_i$,

(3) $\frac{d{n}_i}{\mathit{dt}}=-{\lambda }_i{n}_i-\int \mathit{dE\varphi }\left(\mathit{E}\right){\sigma }_r\left(\mathit{E}\right)+\sum _{j{\ne }_i}{n}_j\left({\lambda }_{\mathit{j}\to \mathit{i}}+\int \mathit{dE\varphi }\left(\mathit{E}\right){\sigma }_{\mathit{j}\to \mathit{i}}\left(\mathit{E}\right)\right).$(3)

Equation (3)(3) $\frac{d{n}_i}{\mathit{dt}}=-{\lambda }_i{n}_i-\int \mathit{dE\varphi }\left(\mathit{E}\right){\sigma }_r\left(\mathit{E}\right)+\sum _{j{\ne }_i}{n}_j\left({\lambda }_{\mathit{j}\to \mathit{i}}+\int \mathit{dE\varphi }\left(\mathit{E}\right){\sigma }_{\mathit{j}\to \mathit{i}}\left(\mathit{E}\right)\right).$(3) can be understood as the sum of losses from the total decay constant ${\lambda }_i$ of $n_i$ and the total reaction rate of any reactions destroying $n_i$ with gains from all other nuclides $n_{\mathit{j}\ne \mathit{i}}$ decaying into $n_i$ based on the decay constant (adjusted for branching) ${\lambda }_{\mathit{j}\to \mathit{i}}$ and the reaction rate of any reactions that transmute $n_j$ into $n_i$ via ${\sigma }_{\mathit{j}\to \mathit{i}}.$ Candidates for ${\sigma }_r$ are essentially all nonscattering reactions, e.g., radiative capture/absorption of a neutron $\left(\mathit{n},\mathit{\gamma }\right)$, fission, and any other reaction that would cause the nuclide to change.

Similarly, ${\sigma }_{\mathit{j}\to \mathit{i}}$ is explicitly only reactions that cause a nuclide $n_j$ to transmute into $n_i$, which severely limits the actual nuclides that have such a valid reaction, e.g., ²³⁸U can undergo radiative capture to become ²³⁹U while ²³⁴U cannot become ²³⁹U in a single radiative capture reaction, as ²³⁴U would transmute to ²³⁵U.

The second equation of interest is the adjoint depletion equation

(4) $-\frac{d{\vec{n}}^{\ast }}{\mathit{dt}}={\mathit{A}}^{\ast }{\vec{n}}^{\ast }+{\vec{S}}^{\ast },$(4)

where $S^{\ast }$ is an arbitrary adjoint source that provides initial importance to the system, and $A^{\ast }$ is the adjoint transmutation matrix. When $\mathit{A}$ is quasi static, $A^{\ast }$ is the transpose of $\mathit{A}$.

To derive the depletion sensitivity coefficient, first-order perturbations can first be applied to Eqs. (2)(2) $\frac{d\vec{n}}{\mathit{dt}}=\mathit{A}\vec{n},$(2) and (4)(4) $-\frac{d{\vec{n}}^{\ast }}{\mathit{dt}}={\mathit{A}}^{\ast }{\vec{n}}^{\ast }+{\vec{S}}^{\ast },$(4) ,^[1]

(5) $\frac{\mathrm{\partial }}{\mathrm{\partial t}}\mathit{\delta }\vec{n}=\mathit{\delta A}\vec{n}+\mathit{A}\mathit{\delta }\vec{n}.$(5)

Multiply Eq. (5)(5) $\frac{\mathrm{\partial }}{\mathrm{\partial t}}\mathit{\delta }\vec{n}=\mathit{\delta A}\vec{n}+\mathit{A}\mathit{\delta }\vec{n}.$(5) by ${\vec{n}}^{\ast }$ and take the inner product to obtain

(6) $⟨{\vec{n}}^{\ast }\frac{\mathrm{\partial }}{\mathrm{\partial t}}\mathit{\delta }\vec{n}⟩=<{\vec{n}}^{\ast }\mathit{\delta A}>\vec{n}+<{\vec{n}}^{\mathbf{\ast }}\mathit{A}\mathit{\delta }\vec{n}>.$(6)

Multiply the adjoint depletion [Eq. (4)(4) $-\frac{d{\vec{n}}^{\ast }}{\mathit{dt}}={\mathit{A}}^{\ast }{\vec{n}}^{\ast }+{\vec{S}}^{\ast },$(4) ] by $\mathit{\delta }\vec{n}$ and take the inner product to obtain

(7) $⟨-\mathit{\delta }\vec{n}\frac{\mathrm{\partial }}{\mathrm{\partial t}}{\vec{n}}^{\ast }⟩=<\mathit{\delta }\vec{n}{\mathit{A}}^{\mathbf{\ast }}{\vec{n}}^{\ast }>+<\mathit{\delta }\vec{n}{\vec{S}}^{\ast }>.$(7)

Using the properties of adjoints and the inner product of a derivative produces an alternate form for Eqs. (6)(6) $⟨{\vec{n}}^{\ast }\frac{\mathrm{\partial }}{\mathrm{\partial t}}\mathit{\delta }\vec{n}⟩=<{\vec{n}}^{\ast }\mathit{\delta A}>\vec{n}+<{\vec{n}}^{\mathbf{\ast }}\mathit{A}\mathit{\delta }\vec{n}>.$(6) and (7)(7) $⟨-\mathit{\delta }\vec{n}\frac{\mathrm{\partial }}{\mathrm{\partial t}}{\vec{n}}^{\ast }⟩=<\mathit{\delta }\vec{n}{\mathit{A}}^{\mathbf{\ast }}{\vec{n}}^{\ast }>+<\mathit{\delta }\vec{n}{\vec{S}}^{\ast }>.$(7) , respectively,

(8) $⟨{\vec{n}}^{\ast }\frac{\mathrm{\partial }}{\mathrm{\partial t}}\mathit{\delta }\vec{n}⟩={\left[n^{\ast }\mathit{\delta n}\right]}_{\mathit{t}=t_i}^{\mathit{t}=t_f}-⟨\mathit{\delta }\vec{n}\frac{\mathrm{\partial }}{\mathrm{\partial t}}{\vec{n}}^{\ast }⟩$(8)

(9) $⟨-\mathit{\delta }\vec{n}\frac{\mathrm{\partial }}{\mathrm{\partial t}}{\vec{n}}^{\ast }⟩=<{\vec{n}}^{\ast }\mathit{A}\mathit{\delta }\vec{n}>+<\mathit{\delta }\vec{n}{\vec{S}}^{\ast }>.$(9)

There are a few ways of joining Eqs. (6)(6) $⟨{\vec{n}}^{\ast }\frac{\mathrm{\partial }}{\mathrm{\partial t}}\mathit{\delta }\vec{n}⟩=<{\vec{n}}^{\ast }\mathit{\delta A}>\vec{n}+<{\vec{n}}^{\mathbf{\ast }}\mathit{A}\mathit{\delta }\vec{n}>.$(6) to (9)(9) $⟨-\mathit{\delta }\vec{n}\frac{\mathrm{\partial }}{\mathrm{\partial t}}{\vec{n}}^{\ast }⟩=<{\vec{n}}^{\ast }\mathit{A}\mathit{\delta }\vec{n}>+<\mathit{\delta }\vec{n}{\vec{S}}^{\ast }>.$(9) together to obtain the depletion sensitivity equation:

1. Use Eqs. (6)(6) $⟨{\vec{n}}^{\ast }\frac{\mathrm{\partial }}{\mathrm{\partial t}}\mathit{\delta }\vec{n}⟩=<{\vec{n}}^{\ast }\mathit{\delta A}>\vec{n}+<{\vec{n}}^{\mathbf{\ast }}\mathit{A}\mathit{\delta }\vec{n}>.$(6) and (8)(8) $⟨{\vec{n}}^{\ast }\frac{\mathrm{\partial }}{\mathrm{\partial t}}\mathit{\delta }\vec{n}⟩={\left[n^{\ast }\mathit{\delta n}\right]}_{\mathit{t}=t_i}^{\mathit{t}=t_f}-⟨\mathit{\delta }\vec{n}\frac{\mathrm{\partial }}{\mathrm{\partial t}}{\vec{n}}^{\ast }⟩$(8) to obtain the left side of Eq. (7)(7) $⟨-\mathit{\delta }\vec{n}\frac{\mathrm{\partial }}{\mathrm{\partial t}}{\vec{n}}^{\ast }⟩=<\mathit{\delta }\vec{n}{\mathit{A}}^{\mathbf{\ast }}{\vec{n}}^{\ast }>+<\mathit{\delta }\vec{n}{\vec{S}}^{\ast }>.$(7) .

2. Equate the common term in Eqs. (6)(6) $⟨{\vec{n}}^{\ast }\frac{\mathrm{\partial }}{\mathrm{\partial t}}\mathit{\delta }\vec{n}⟩=<{\vec{n}}^{\ast }\mathit{\delta A}>\vec{n}+<{\vec{n}}^{\mathbf{\ast }}\mathit{A}\mathit{\delta }\vec{n}>.$(6) and (9)(9) $⟨-\mathit{\delta }\vec{n}\frac{\mathrm{\partial }}{\mathrm{\partial t}}{\vec{n}}^{\ast }⟩=<{\vec{n}}^{\ast }\mathit{A}\mathit{\delta }\vec{n}>+<\mathit{\delta }\vec{n}{\vec{S}}^{\ast }>.$(9) .

3. Apply the product rule to obtain a result that utilizes Eqs. (6)(6) $⟨{\vec{n}}^{\ast }\frac{\mathrm{\partial }}{\mathrm{\partial t}}\mathit{\delta }\vec{n}⟩=<{\vec{n}}^{\ast }\mathit{\delta A}>\vec{n}+<{\vec{n}}^{\mathbf{\ast }}\mathit{A}\mathit{\delta }\vec{n}>.$(6) and (7)(7) $⟨-\mathit{\delta }\vec{n}\frac{\mathrm{\partial }}{\mathrm{\partial t}}{\vec{n}}^{\ast }⟩=<\mathit{\delta }\vec{n}{\mathit{A}}^{\mathbf{\ast }}{\vec{n}}^{\ast }>+<\mathit{\delta }\vec{n}{\vec{S}}^{\ast }>.$(7) .

All three methods yield the same result, and method 3 is used here. Starting with the product rule,

(10) $⟨\frac{\mathrm{\partial }}{\mathrm{\partial t}}\left({\vec{n}}^{\ast }\mathit{\delta }\vec{n}\right)⟩=⟨{\vec{n}}^{\ast }\frac{\mathrm{\partial }}{\mathrm{\partial t}}\left(\mathit{\delta }\vec{n}\right)⟩+⟨\mathit{\delta }\vec{n}\frac{\mathrm{\partial }}{\mathrm{\partial t}}\left({\vec{n}}^{\ast }\right)⟩,$(10)

Eq. (10)(10) $⟨\frac{\mathrm{\partial }}{\mathrm{\partial t}}\left({\vec{n}}^{\ast }\mathit{\delta }\vec{n}\right)⟩=⟨{\vec{n}}^{\ast }\frac{\mathrm{\partial }}{\mathrm{\partial t}}\left(\mathit{\delta }\vec{n}\right)⟩+⟨\mathit{\delta }\vec{n}\frac{\mathrm{\partial }}{\mathrm{\partial t}}\left({\vec{n}}^{\ast }\right)⟩,$(10) is equal to Eq. (7)(7) $⟨-\mathit{\delta }\vec{n}\frac{\mathrm{\partial }}{\mathrm{\partial t}}{\vec{n}}^{\ast }⟩=<\mathit{\delta }\vec{n}{\mathit{A}}^{\mathbf{\ast }}{\vec{n}}^{\ast }>+<\mathit{\delta }\vec{n}{\vec{S}}^{\ast }>.$(7) subtracted from Eq. (6)(6) $⟨{\vec{n}}^{\ast }\frac{\mathrm{\partial }}{\mathrm{\partial t}}\mathit{\delta }\vec{n}⟩=<{\vec{n}}^{\ast }\mathit{\delta A}>\vec{n}+<{\vec{n}}^{\mathbf{\ast }}\mathit{A}\mathit{\delta }\vec{n}>.$(6) ,

(11) $⟨\frac{\mathrm{\partial }}{\mathrm{\partial t}}\left({\vec{n}}^{\ast }\mathit{\delta }\vec{n}\right)⟩=<{\vec{n}}^{\ast }\mathit{\delta A}\vec{n}>-<\mathit{\delta }\vec{n}{\vec{S}}^{\ast }>.$(11)

Using the fundamental theorem of calculus to evaluate the time derivative in the inner product yields

(12) $⟨\frac{\mathrm{\partial }}{\mathrm{\partial t}}\left({\vec{n}}^{\ast }\mathit{\delta }\vec{n}\right)⟩=⟨{\left[{\vec{n}}^{\ast }\mathit{\delta }\vec{n}\right]}_{\mathit{t}=t_i}^{\mathit{t}=t_f}⟩.$(12)

The inner product on the right side of Eq. (12)(12) $⟨\frac{\mathrm{\partial }}{\mathrm{\partial t}}\left({\vec{n}}^{\ast }\mathit{\delta }\vec{n}\right)⟩=⟨{\left[{\vec{n}}^{\ast }\mathit{\delta }\vec{n}\right]}_{\mathit{t}=t_i}^{\mathit{t}=t_f}⟩.$(12) is no longer a formal inner product, as the time integration has been removed, and the inner product notation is preserved in Eq. (12)(12) $⟨\frac{\mathrm{\partial }}{\mathrm{\partial t}}\left({\vec{n}}^{\ast }\mathit{\delta }\vec{n}\right)⟩=⟨{\left[{\vec{n}}^{\ast }\mathit{\delta }\vec{n}\right]}_{\mathit{t}=t_i}^{\mathit{t}=t_f}⟩.$(12) for convenience.

Using Eqs. (11)(11) $⟨\frac{\mathrm{\partial }}{\mathrm{\partial t}}\left({\vec{n}}^{\ast }\mathit{\delta }\vec{n}\right)⟩=<{\vec{n}}^{\ast }\mathit{\delta A}\vec{n}>-<\mathit{\delta }\vec{n}{\vec{S}}^{\ast }>.$(11) and (12)(12) $⟨\frac{\mathrm{\partial }}{\mathrm{\partial t}}\left({\vec{n}}^{\ast }\mathit{\delta }\vec{n}\right)⟩=⟨{\left[{\vec{n}}^{\ast }\mathit{\delta }\vec{n}\right]}_{\mathit{t}=t_i}^{\mathit{t}=t_f}⟩.$(12) , and assuming the perturbations are relative to $\mathit{\delta \alpha }/\mathit{\alpha }$, yields the number density sensitivity coefficient, where $\mathit{\alpha }$ represents an uncertain nuclear data parameter (cross section or initial nuclide density), as^[1]

(13) $⟨\frac{\delta \vec{n}}{\mathit{\delta \alpha }/\mathit{\alpha }}{\vec{S}}^{\ast }⟩=⟨{\vec{n}}^{\ast }\frac{\mathit{\delta A}}{\mathit{\delta \alpha }/\alpha }\vec{n}⟩-⟨{\left[{\vec{n}}^{\ast }\frac{\delta \vec{n}}{\mathit{\delta \alpha }/\mathit{\alpha }}\right]}_{\mathit{t}=t_i}^{\mathit{t}=t_f}⟩.$(13)

This yields the depletion sensitivity equation present in Ref. [1]. The $\mathit{\delta A}$ term inherently depends on the flux and reaction rates present in the depletion operator. GPT will next account for this dependency.

The GEAR-MC/CLUTCH GPT method estimates the sensitivity of the reaction rates used in the Bateman depletion matrix $\mathit{A}$ using the quasi-static Linear Boltzmann Transport Equation B^[2],

(14) $⟨\frac{\mathit{\delta A}}{\mathit{\delta \alpha }/\alpha }\vec{n}⟩=⟨{\mathrm{\Gamma }}^{\ast }\frac{\mathrm{\partial }\mathit{B}}{\mathrm{\partial \alpha }/\alpha }\mathit{\varphi }⟩.$(14)

Equation (14)(14) $⟨\frac{\mathit{\delta A}}{\mathit{\delta \alpha }/\alpha }\vec{n}⟩=⟨{\mathrm{\Gamma }}^{\ast }\frac{\mathrm{\partial }\mathit{B}}{\mathrm{\partial \alpha }/\alpha }\mathit{\varphi }⟩.$(14) is extended using total derivatives to account for sensitivity effects from perturbations to nuclide density,

(15) $⟨\frac{\mathit{\delta A}}{\mathit{\delta \alpha }/\alpha }\vec{n}⟩=⟨{\mathrm{\Gamma }}^{\ast }\frac{\mathrm{\partial }\mathit{B}}{\mathrm{\partial \alpha }/\alpha }\mathit{\varphi }⟩+⟨{\mathrm{\Gamma }}^{\ast }\frac{\mathrm{\partial }\mathit{B}}{\mathrm{\partial }\vec{n}/\vec{n}}\frac{\mathrm{\partial }\vec{n}/\vec{n}}{\mathrm{\partial \alpha }/\mathit{\alpha }}\mathit{\varphi }⟩.$(15)

This adjustment suggests that nuclear data perturbations introduce a direct-effect perturbation to reaction rates, the first term, and an indirect-effect through perturbing the nuclide densities in these reaction rates, the second term. The $\mathit{\alpha }$ in Eq. (15)(15) $⟨\frac{\mathit{\delta A}}{\mathit{\delta \alpha }/\alpha }\vec{n}⟩=⟨{\mathrm{\Gamma }}^{\ast }\frac{\mathrm{\partial }\mathit{B}}{\mathrm{\partial \alpha }/\alpha }\mathit{\varphi }⟩+⟨{\mathrm{\Gamma }}^{\ast }\frac{\mathrm{\partial }\mathit{B}}{\mathrm{\partial }\vec{n}/\vec{n}}\frac{\mathrm{\partial }\vec{n}/\vec{n}}{\mathrm{\partial \alpha }/\mathit{\alpha }}\mathit{\varphi }⟩.$(15) is similarly any nuclear data parameter of interest, as described previously.

The Linear Boltzmann Transport Equation can be defined as

(16) $\begin{array}{rl}& \frac{1}{v}\frac{\mathrm{\partial }\psi (\vec{r},\mathit{E},\mathrm{\Omega },\mathrm{t})}{\mathrm{\partial }\mathit{t}}=\begin{array}{c}\mathit{Q}\left(\vec{r},\mathit{E},\mathrm{\Omega },\mathrm{t}\right)\\ \mathit{External}\mathit{source}\end{array}+\begin{array}{c}\int \int \mathit{\psi }\left(\vec{r},E^{\prime },\mathrm{\Omega }{ }^{\mathrm{\prime }},\mathrm{t}\right){\mathrm{\Sigma }}_s\left(\vec{r},\mathit{E}{ }^{\mathrm{\prime }}\to \mathit{E},\mathrm{\Omega }{ }^{\mathrm{\prime }}\cdot \mathrm{\Omega }\right)\mathit{d}\mathrm{\Omega }{ }^{\mathrm{\prime }}\mathrm{d}{\mathrm{E}}^{\prime }\\ \mathit{Scattering}\end{array}\\ & +\begin{array}{c}\frac{{\displaystyle \chi \left(\vec{r},\mathit{E}\right)}}{{\displaystyle 4\mathit{\pi }}}\int \int \mathit{\nu }{\mathrm{\Sigma }}_f\left(\vec{r},\mathit{E}{ }^{\mathrm{\prime }}\right)\mathit{\psi }\left(\vec{r},E^{\prime },\mathrm{\Omega }{ }^{\mathrm{\prime }},\mathrm{t}\right)\mathit{d}\mathrm{\Omega }{ }^{\mathrm{\prime }}\mathrm{dE}{ }^{\mathrm{\prime }}\\ \mathit{Fission}/\mathit{Multiplication}\end{array}-\left[\begin{array}{c}\mathrm{\Omega }\cdot \mathrm{\nabla }\\ \mathit{Leakage}\end{array}+\begin{array}{c}{\mathrm{\Sigma }}_t\left(\vec{r},\mathit{E}\right)\\ \mathit{Collisions}\end{array}\right]\mathit{\psi }(\vec{r},\mathit{E},\mathrm{\Omega },\mathrm{t})& (16)\end{array}$(16)

where $\mathit{r}$ represents position, $\mathit{E}$ represents energy, $\mathrm{\Omega }$ represents direction, and $\mathit{t}$ represents time. The quasi-static assumption relaxes the time dependence into a piecewise function, with the time dependence being established in ${\mathrm{\Sigma }}_s,{\mathrm{\Sigma }}_f,{\mathrm{\Sigma }}_t$ as the nuclide field (and thus cross sections) evolve with time, consequentially causing the angular flux $\mathit{\psi }$ to change. The left side of Eq. (16)(16) $\begin{array}{rl}& \frac{1}{v}\frac{\mathrm{\partial }\psi (\vec{r},\mathit{E},\mathrm{\Omega },\mathrm{t})}{\mathrm{\partial }\mathit{t}}=\begin{array}{c}\mathit{Q}\left(\vec{r},\mathit{E},\mathrm{\Omega },\mathrm{t}\right)\\ \mathit{External}\mathit{source}\end{array}+\begin{array}{c}\int \int \mathit{\psi }\left(\vec{r},E^{\prime },\mathrm{\Omega }{ }^{\mathrm{\prime }},\mathrm{t}\right){\mathrm{\Sigma }}_s\left(\vec{r},\mathit{E}{ }^{\mathrm{\prime }}\to \mathit{E},\mathrm{\Omega }{ }^{\mathrm{\prime }}\cdot \mathrm{\Omega }\right)\mathit{d}\mathrm{\Omega }{ }^{\mathrm{\prime }}\mathrm{d}{\mathrm{E}}^{\prime }\\ \mathit{Scattering}\end{array}\\ & +\begin{array}{c}\frac{{\displaystyle \chi \left(\vec{r},\mathit{E}\right)}}{{\displaystyle 4\mathit{\pi }}}\int \int \mathit{\nu }{\mathrm{\Sigma }}_f\left(\vec{r},\mathit{E}{ }^{\mathrm{\prime }}\right)\mathit{\psi }\left(\vec{r},E^{\prime },\mathrm{\Omega }{ }^{\mathrm{\prime }},\mathrm{t}\right)\mathit{d}\mathrm{\Omega }{ }^{\mathrm{\prime }}\mathrm{dE}{ }^{\mathrm{\prime }}\\ \mathit{Fission}/\mathit{Multiplication}\end{array}-\left[\begin{array}{c}\mathrm{\Omega }\cdot \mathrm{\nabla }\\ \mathit{Leakage}\end{array}+\begin{array}{c}{\mathrm{\Sigma }}_t\left(\vec{r},\mathit{E}\right)\\ \mathit{Collisions}\end{array}\right]\mathit{\psi }(\vec{r},\mathit{E},\mathrm{\Omega },\mathrm{t})& (16)\end{array}$(16) is treated as zero in the quasi-static assumption, as the system is treated as being in a steady state (and thus in balance). The angular flux $\mathit{\psi }$ is integrated over all angles to yield flux $\mathit{\varphi }$. $\mathit{Q}$, the external source, is treated as zero. To ensure balance, a factor of $\frac{1}{{\mathit{k}}_{\mathit{eff}}}$ is applied to the fission/multiplication term.

The GEAR-MC/CLUTCH methodology is used to compute the GPT sensitivity coefficients from Eq. (16)(16) $\begin{array}{rl}& \frac{1}{v}\frac{\mathrm{\partial }\psi (\vec{r},\mathit{E},\mathrm{\Omega },\mathrm{t})}{\mathrm{\partial }\mathit{t}}=\begin{array}{c}\mathit{Q}\left(\vec{r},\mathit{E},\mathrm{\Omega },\mathrm{t}\right)\\ \mathit{External}\mathit{source}\end{array}+\begin{array}{c}\int \int \mathit{\psi }\left(\vec{r},E^{\prime },\mathrm{\Omega }{ }^{\mathrm{\prime }},\mathrm{t}\right){\mathrm{\Sigma }}_s\left(\vec{r},\mathit{E}{ }^{\mathrm{\prime }}\to \mathit{E},\mathrm{\Omega }{ }^{\mathrm{\prime }}\cdot \mathrm{\Omega }\right)\mathit{d}\mathrm{\Omega }{ }^{\mathrm{\prime }}\mathrm{d}{\mathrm{E}}^{\prime }\\ \mathit{Scattering}\end{array}\\ & +\begin{array}{c}\frac{{\displaystyle \chi \left(\vec{r},\mathit{E}\right)}}{{\displaystyle 4\mathit{\pi }}}\int \int \mathit{\nu }{\mathrm{\Sigma }}_f\left(\vec{r},\mathit{E}{ }^{\mathrm{\prime }}\right)\mathit{\psi }\left(\vec{r},E^{\prime },\mathrm{\Omega }{ }^{\mathrm{\prime }},\mathrm{t}\right)\mathit{d}\mathrm{\Omega }{ }^{\mathrm{\prime }}\mathrm{dE}{ }^{\mathrm{\prime }}\\ \mathit{Fission}/\mathit{Multiplication}\end{array}-\left[\begin{array}{c}\mathrm{\Omega }\cdot \mathrm{\nabla }\\ \mathit{Leakage}\end{array}+\begin{array}{c}{\mathrm{\Sigma }}_t\left(\vec{r},\mathit{E}\right)\\ \mathit{Collisions}\end{array}\right]\mathit{\psi }(\vec{r},\mathit{E},\mathrm{\Omega },\mathrm{t})& (16)\end{array}$(16) and is used to obtain the $\frac{\mathrm{\partial B}}{\mathrm{\partial \alpha }/\alpha }$ terms in Eqs. (14)(14) $⟨\frac{\mathit{\delta A}}{\mathit{\delta \alpha }/\alpha }\vec{n}⟩=⟨{\mathrm{\Gamma }}^{\ast }\frac{\mathrm{\partial }\mathit{B}}{\mathrm{\partial \alpha }/\alpha }\mathit{\varphi }⟩.$(14) and (15)(15) $⟨\frac{\mathit{\delta A}}{\mathit{\delta \alpha }/\alpha }\vec{n}⟩=⟨{\mathrm{\Gamma }}^{\ast }\frac{\mathrm{\partial }\mathit{B}}{\mathrm{\partial \alpha }/\alpha }\mathit{\varphi }⟩+⟨{\mathrm{\Gamma }}^{\ast }\frac{\mathrm{\partial }\mathit{B}}{\mathrm{\partial }\vec{n}/\vec{n}}\frac{\mathrm{\partial }\vec{n}/\vec{n}}{\mathrm{\partial \alpha }/\mathit{\alpha }}\mathit{\varphi }⟩.$(15) . These methods require estimating intergenerational and intragenerational importance terms, which describe, respectively, how the flight path and interactions of neutrons in the current generation generate a response of interest, and how the current generation of neutrons birth fission neutrons that themselves generate a response of interest. This is a memory intensive process that allows for the generalized flux importance ${\mathrm{\Gamma }}^{\ast }$ to be computed, although these algorithms are shown to scale well to systems with a large number of responses.^[4]

II.A. Power Normalization Sensitivity Effects

Finally, a power normalization correction is introduced based on Ref. [1] that accounts for how perturbations affect the flux power normalization during each depletion step. Power normalization is important to address, as solving for the flux profile (or the reaction rate profile) only specifies the shape of the flux and does not constrain the peak value of the flux. Reactor power is used as an outside constraint to determine the peak value of the flux. In this case, the fission reaction rate was utilized to determine the power of the reactor, as each fission yields some amount of energy. Using the notation from the OpenMC documentation^[5]

(17) $H^{\prime }=1.602\mathit{x}{10}^{-19}\left[\frac{\mathrm{J}}{\mathrm{eV}}\right]\times \mathit{H}\left[\frac{\mathrm{eV}}{\mathit{source}}\right],$(17)

where $\mathit{H}$ is the heating rate and $H^{\prime }$ is a normalization factor of the heating rate in units of joules per source particle. Energy deposition from absorption and scattering reactions or from the decay of nuclides is neglected when computing reactor power to simplify the system. The normalization factor is then combined with a desired power target $\mathit{P}$ to find what factor $\mathit{f}$ all flux scores will be scaled by

(18) $\mathit{f}=\mathit{\hspace{0.17em}}\frac{P}{H^{\prime }}=\frac{\left[\mathrm{J}/\mathrm{s}\right]}{\left[\mathrm{J}/s\mathrm{ource}\right]}=\left[\frac{\mathit{source}}{\mathrm{s}}\right].$(18)

The normalized flux can then be obtained in typical units

(19) $\mathit{\varphi }\mathrm{\prime }=\frac{\mathit{f\varphi }}{V}=\frac{\left[\mathit{source}/\mathrm{s}\right]\left[\mathit{particle}\ast \mathrm{cm}/s\mathrm{ource}\right]}{\left[\mathrm{c}{\mathrm{m}}^3\right]}=\left[\frac{\mathit{particle}}{\mathrm{c}{\mathrm{m}}^2\ast \mathrm{s}}\right].$(19)

In the case of a fixed-source simulation, $\mathit{f}$ is prescribed or already known (as was the case with the implementation in Ref. [3]). Equations (17)(17) $H^{\prime }=1.602\mathit{x}{10}^{-19}\left[\frac{\mathrm{J}}{\mathrm{eV}}\right]\times \mathit{H}\left[\frac{\mathrm{eV}}{\mathit{source}}\right],$(17) , (18)(18) $\mathit{f}=\mathit{\hspace{0.17em}}\frac{P}{H^{\prime }}=\frac{\left[\mathrm{J}/\mathrm{s}\right]}{\left[\mathrm{J}/s\mathrm{ource}\right]}=\left[\frac{\mathit{source}}{\mathrm{s}}\right].$(18) , and (19)(19) $\mathit{\varphi }\mathrm{\prime }=\frac{\mathit{f\varphi }}{V}=\frac{\left[\mathit{source}/\mathrm{s}\right]\left[\mathit{particle}\ast \mathrm{cm}/s\mathrm{ource}\right]}{\left[\mathrm{c}{\mathrm{m}}^3\right]}=\left[\frac{\mathit{particle}}{\mathrm{c}{\mathrm{m}}^2\ast \mathrm{s}}\right].$(19) can be understood as finding a simple factor to multiply the flux by to achieve the desired power target.

Williams proposed a method for estimating the power normalization term by first estimating an adjoint function for the power norm during time step $\mathit{i}$,

(20) $P_i^{\ast }\left(\mathrm{Williams}\right)=\frac{{\int }_{t_i}^{t_{\mathit{i}+1}}<N^{\ast }<{\varphi }_i\mathit{R}{>}_{\mathrm{\Omega },\mathit{E}}\mathit{N}{>}_V\mathit{dt}+\frac{\mathrm{\partial R}}{\mathrm{\partial }{\varphi }_i}}{<\mathit{M}{\varphi }_i{\sigma }_f{N_i>}_{\mathrm{\Omega },\mathit{E},\mathit{V}}},$(20)

where the $\frac{\mathrm{\partial R}}{\mathrm{\partial }{\varphi }_i}$ term equals zero when estimating nuclide field sensitivities because flux sensitivity effects are captured in the GPT indirect-effect term. The $P_i^{\ast }$ term can estimate the power normalization sensitivity component $S_{R,\alpha }^{po\mathrm{wer}\mathit{norm}}$,

(21) $S_{R,\alpha }^{po\mathrm{wer}\mathit{norm}}\left(\mathit{Williams}\right)=\frac{\alpha }{R}\sum _i{P}_i^{\ast }\left[{\varphi }_i\frac{\mathrm{\partial }}{\mathrm{\partial }\mathit{\alpha }}\left({\sigma }_f{N}_i\right)\right].$(21)

Williams provided a sample calculation for the power normalization sensitivity effects in Ref. [1], but this example both implies that activation or destruction of nonfissioning isotopes affects the fission power normalization and also does not discriminate between the importance effects from separate isotopes.

Initial prototyping in OpenMC employed Williams’s derivation and found it to produce discrepant results. Therefore, this work develops an adjustment to Williams’s derivation for power normalization sensitivity effects.

The adjoint of the power normalization term $P_i^{\ast }$ during time step $\mathit{i}$ must be separated into sensitivity components for each isotope $\mathit{j}$,

(22) $P_{\mathit{i},\mathit{j}}=\frac{\mathrm{\Delta }{\mathit{N}}_{\mathit{i},\mathit{j}}{N}_{\mathit{i},\mathit{j}\to \mathit{j}}}{<\mathit{M}{\varphi }_i{\sigma }_{\mathit{f},\mathit{j}}{N}_{\mathit{i},\mathit{j}}{>}_{\mathrm{\Omega },\mathit{E},\mathit{V}}}.$(22)

This term accounts for the sensitivity of the production/destruction rate for isotope $\mathit{j}$ to the system’s overall fission rate. For example, isotopes that have reached equilibrium would not be impacted by changes in a system’s power normalization, while isotopes that are transmuting rapidly would experience slower transmutation rates if the power normalization constant decreased due to increases in fission cross sections.

This definition of $P_{i,j}^{\ast }$ differs from Williams’s $P_i^{\ast }$ both in that it includes a $N_{i,j\to j}^{\ast }$ term, which was found to produce more accurate sensitivity estimates in a semi-empirical study, and that $P_{i,j}^{\ast }$ depends on the isotope $\mathit{j}$. Williams may have intended for his $P_i^{\ast }$ definition to incorporate this nuclide dependency, but it was not clear to the authors whether this was the case. Williams included a sample depletion sensitivity calculation in Ref. [1], but this simple case did not require an isotope-dependent $P_i^{\ast }$ to yield correct sensitivity results.

From here, the power norm sensitivity component of isotope $\mathit{j}$ during time step $\mathit{i}$ can be estimated by

(23) $S_{{\mathit{N}}_{\mathit{i},\mathit{j}},\mathit{\alpha }}^{\mathit{power}\mathit{norm}}=\sum _{\mathit{k}=\mathit{fissile}\mathit{isotopes}}\frac{P_{i,j}^{\ast }}{N_{\mathit{i},\mathit{j}}}\left[{\varphi }_i{\sigma }_{\mathit{f},\mathit{k}}{N}_{\mathit{i},\mathit{k}}{S}_{{\mathit{N}}_{\mathit{i}-1,\mathit{k}},\mathit{\alpha }}\right].$(23)

In essence, this term weights the importance of changes to the system’s overall fission rate relative to an isotope’s production/destruction by the sensitivity of the overall fission rate to the parameter $\mathit{\alpha }$. The fission product yields utilized are independent of each other in that they can be understood as a single fission produces an average a fractional amount of all nuclides described by the parent nuclide’s fission yield.

Williams presented a similar definition in Eq. (21)(21) $S_{R,\alpha }^{po\mathrm{wer}\mathit{norm}}\left(\mathit{Williams}\right)=\frac{\alpha }{R}\sum _i{P}_i^{\ast }\left[{\varphi }_i\frac{\mathrm{\partial }}{\mathrm{\partial }\mathit{\alpha }}\left({\sigma }_f{N}_i\right)\right].$(21) that included the fission cross-section term ${\sigma }_f$ in the innermost parenthesis of the expression. This introduced an unnecessary direct-effect term when $\mathit{\alpha }$ represents a fission cross section, similar to how the direct-effect term is accounted for in GPT sensitivity calculations. Williams’s example depletion sensitivity calculations were performed only for changes in initial number densities, meaning that they could not have detected this error in Eq. (21)(21) $S_{R,\alpha }^{po\mathrm{wer}\mathit{norm}}\left(\mathit{Williams}\right)=\frac{\alpha }{R}\sum _i{P}_i^{\ast }\left[{\varphi }_i\frac{\mathrm{\partial }}{\mathrm{\partial }\mathit{\alpha }}\left({\sigma }_f{N}_i\right)\right].$(21) .

III. PROOF OF PRINCIPLE

A simple model based off a PWR pin cell (shown in Fig. 2) was utilized to demonstrate the depletion sensitivity analysis methodology. The 2.4% enrichment UO₂ fuel pin was depleted for 31 days over eight burn steps at a typical PWR power density of 90 kW/L. Prior to implementing these methods in OpenMC, a simple 1-D, three-group, three-nuclide fixed-source slab model was used to demonstrate this methodology.^[3] This initial test case showed good agreement, thus prompting the implementation into OpenMC.

Fig. 2. UO₂ pin cell with a helium air gap and zircalloy cladding surrounded by water.

The burnable nuclide inventory consisted of 12 nuclides, ${ }^{234}\mathrm{U}$, ${ }^{235}\mathrm{U}$, ${ }^{238}\mathrm{U}$, ${ }^{239}\mathrm{U}$, ${ }^{239}\mathrm{Np}$, and ${ }^{239}\mathrm{Pu}$, and the following fission products, ${ }^{135}\mathrm{I}$, ${ }^{135}\mathrm{Cs}$, ${ }^{135}\mathrm{Xe}$, ${ }^{136}\mathrm{Xe}$, ${ }^{156}\mathrm{Gd}$, and ${ }^{157}\mathrm{Gd}$. This abbreviated list of nuclides was selected to limit computational resources when implementing and testing this methodology, and includes nuclides with short half-lives (${ }^{239}\mathrm{U}$ and ${ }^{135}\mathrm{I}$), significant transmutation cross sections (${ }^{135}\mathrm{Xe}{\to }^{136}\mathrm{Xe}),$ and stable fission products (${ }^{156}\mathrm{Gd}$ and ${ }^{157}\mathrm{Gd}$). The ENDFB-VII.1 nuclear data library in OpenMC was used for all simulations.^[6] The fission products and ${ }^{239}\mathrm{U}$, ${ }^{239}\mathrm{Np}$, and ${ }^{239}\mathrm{Pu}$ existed as only trace quantities at the start of the simulation for the purpose of tallying nonzero GPT responses.

A few changes were made to the depletion/transmutation physics for simplicity. As ${ }^{136}\mathrm{I}$ has a half-life of 83.6 s, this decay was treated as if it occurs instantaneously. This means that any neutron absorption of ${ }^{135}\mathrm{I}$ was presumed to directly transmute to ${ }^{136}\mathrm{Xe}$, instead of producing ${ }^{136}\mathrm{I}$ and then decaying (or undergoing any transmutation reactions). Alpha decay was neglected, as the half-lives of ${ }^{234}\mathrm{U}$, ${ }^{235}\mathrm{U}$, ${ }^{238}\mathrm{U}$, and ${ }^{239}\mathrm{Pu}$ are long in comparison to the total simulation time. Similarly, the fission yields were truncated to match the abbreviated list of nuclides. If a nuclide transmutes or decays into a nuclide that is not in the burnable nuclide inventory, it was removed from the system.

The method was implemented using a simple predictor depletion scheme and utilizing the Chebyshev Rational Approximation Method to perform the forward depletion.^[7] Because the depletion operator is quasi-static for the predictor scheme, solving for the adjoint depletion operator ${\mathit{A}}^{\mathbf{\ast }}$ simply required transposing the original depletion matrix $\mathit{A}$. While this method could be applied to predictor-corrector simulations, it has not yet been tested but could be naively applied to the corrector step.The method results were compared against sensitivity coefficients produced from reference direct perturbation calculations, which were examined to confirm linearity.

III.A. Burn History

Table I describes the 31-day burn history of the pin cell. During each depletion step, a Monte Carlo transport calculation was performed to calculate the GPT sensitivity coefficients and to tally the reaction rates. The numerator for the GPT reaction rate sensitivities was the reaction rate of interest (e.g., fission, radiative capture), while the denominator was the total fission rate for the system. A total of 50 million particles were simulated per Monte Carlo depletion step.

TABLE I Depletion Time-Step Profile of the Test Problem

Depletion step	1	2	3	4	5	6	7	8
Burn duration (days)	0.1	0.3	0.6	1.0	3.0	6.0	10.0	10.0
Total burn time (days)	0.1	0.4	1.0	2.0	5.0	11.0	21.0	31.0

III.B. Depletion Sensitivity Results

In the following figures (Figs. 3-10 and Figs.,13-15), the results of various depletion sensitivity coefficients are presented, grouped by the type of sensitivity coefficient. The DPT sensitivity coefficients were compared against the sensitivity coefficients produced by direct perturbation of $\pm 1\mathrm{\%}$, $3\mathrm{\%}$, $5\mathrm{\%}$, and $10\mathrm{\%}$. The DPT sensitivity coefficient is presented in green, while the direct perturbation result is in blue, with red lines representing the standard error in the direct perturbation linear regression. While a 10% perturbation is large, the direct perturbations were generally found be linear. A comment regarding the quality of the DPT results based on a graph of $\left(\frac{\mathit{\delta R}}{R}\right)$ versus $\left(\frac{\mathit{\delta \alpha }}{\alpha }\right)$ was made to address this concern on a per result basis.

Fig. 3. Sensitivity of nuclide inventories to the ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ fission cross section.

Fig. 4. Sensitivity of nuclide inventories to the ${\mathit{\hspace{0.17em}}}^{238}\mathit{U}\left(\mathit{n},\mathit{\gamma }\right)$ cross section.

Fig. 5. Sensitivity of nuclide inventories to the ${\mathit{\hspace{0.17em}}}^{234}\mathrm{U}\left(\mathit{n},\mathit{\gamma }\right)$ cross section.

Fig. 6. Energy-dependent number density sensitivity of ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$ to the (left) ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ fission cross section and (right) the${\mathit{\hspace{0.17em}}}^{238}\mathrm{U}\left(\mathit{n},\mathit{\gamma }\right)$ cross section.

Fig. 7. Sensitivity of nuclide inventories to the initial ${\mathit{\hspace{0.17em}}}^{238}\mathrm{U}$ number density.

Fig. 8. Sensitivity of the nuclide inventories to the initial ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ number density.

Fig. 9. Sensitivity of the nuclide inventories to the initial ${\mathit{\hspace{0.17em}}}^{234}\mathrm{U}$ number density.

Fig. 10. Sensitivity of the nuclide inventories to the ${\mathit{\hspace{0.17em}}}^{239}\mathrm{U}$ decay constant.

The initial inventory of depleting nuclides was dominated by uranium isotopes: ²³⁸U comprised 97.55 at. % of the initial burnable inventory, ²³⁵U 2.43%, and ²³⁴U 0.02%. Thus, the initial number density perturbations were based on these three nuclides.

Three primary mechanisms affecting the sensitivity coefficients can be observed in the graphs: direct effects, indirect effects, and power normalization. Direct effects can be understood as the direct production (or destruction) of nuclides through the reaction of interest, while indirect effects affect the production (or destruction) of nuclides by inducing perturbations to the neutron flux. The power normalization describes how changes in the fission profile impacts the normalization of the flux and thus overall rate of transmutation occurring in the system.

III.B.1. Nuclide Inventory Sensitivity to Cross Sections

This section discusses the sensitivity of the nuclide inventory at various time steps to specific cross sections: the ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ fission cross section, the ${\mathit{\hspace{0.17em}}}^{238}\mathrm{U}$ $\left(\mathit{n},\mathit{\gamma }\right)$ cross section, and the ${\mathit{\hspace{0.17em}}}^{234}\mathrm{U}$ $\left(\mathit{n},\mathit{\gamma }\right)$ cross section. DPT sensitivity coefficients are presented both with and without power normalization to quantify the impact power normalization has on the produced sensitivity coefficients. For sensitivity coefficients that are not directly tied to the fission rate (such as absorption), power normalization generally shows a minimal impact. On the other hand, fission cross sections and fissile nuclide number densities are greatly impacted by power normalization.

Most of the sensitivity coefficients exhibit strong agreement with the direct perturbation results, with the exception of ${\mathit{\hspace{0.17em}}}^{239}\mathrm{U}$, ${\mathit{\hspace{0.17em}}}^{135}\mathrm{I}$, and ${\mathit{\hspace{0.17em}}}^{156}\mathrm{Gd}$ (Fig. 3). The observed disagreement for ${\mathit{\hspace{0.17em}}}^{156}\mathrm{Gd}$ is possibly stochastic noise, while the disagreement observed for ${\mathit{\hspace{0.17em}}}^{239}\mathrm{U}$ and ${\mathit{\hspace{0.17em}}}^{135}\mathrm{I}$ is thought to be an artifact of the DPT sensitivity coefficient integration scheme for these relatively short-lived isotopes. Studies during the initial OpenMC DPT method development encountered similar discrepancies and succeeded in resolving this discrepancy by significantly increasing the resolution of the DPT sensitivity integration scheme. The power correction term provided a minor correcting factor but was not particularly prominent.

The ${\mathit{\hspace{0.17em}}}^{234}\mathrm{U}$ number density is only indirectly impacted by the ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ fission cross section absorbing the flux. The ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ number density is strongly proportional to direct effect of ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ fission cross sections. Increasing ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ cross sections has a similar direct effect on the fission product densities. Conversely, increasing ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ fission cross sections reduces the amount of flux seen by ${\mathit{\hspace{0.17em}}}^{238}\mathrm{U}$ and significantly hinders the buildup of isotopes in the ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$ activation chain.

Good agreement was again observed for all isotopes except ${\mathit{\hspace{0.17em}}}^{239}\mathrm{U}$, ${\mathit{\hspace{0.17em}}}^{135}\mathrm{I}$, and ${\mathit{\hspace{0.17em}}}^{156}\mathrm{Gd}$ (Fig. 4). Significant discrepancy existed for the ${\mathit{\hspace{0.17em}}}^{239}\mathrm{U}$ sensitivity estimates, but ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Np}$ and ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$ both agreed well with the direct perturbation result. The fission products were partially directly derived from the ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$ buildup from the ${\mathit{\hspace{0.17em}}}^{238}\mathrm{U}\left(\mathit{n},\mathit{\gamma }\right)$ cross section, but were mainly generated by the fission of ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$, which was indirectly impacted by flux self-shielding from ${\mathit{\hspace{0.17em}}}^{238}\mathrm{U}$ capture reactions.

The direct perturbation results in Fig. 5 were strongly contaminated by stochastic noise except for the ${\mathit{\hspace{0.17em}}}^{234}\mathrm{U}$ and ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ sensitivities; however, all DPT sensitivity estimates agreed well with the direct perturbation results within this noise. This stochastic noise was not surprising, as ${\mathit{\hspace{0.17em}}}^{234}\mathrm{U}$ comprised less than 1% of the total uranium, and therefore produced minor sensitivity effects. That said, the fact that the DPT sensitivity coefficients agreed well with these noisy, but small in magnitude, results speaks highly to the DPT accuracy. These sensitivity coefficients span a range in magnitude from 10⁻³ to 10⁻⁷, and any small bias or error in the DPT sensitivity coefficient methodology would likely be readily apparent.

The adjoint-weighted DPT methodology also allowed for the estimation of high-resolution, energy-dependent sensitivity coefficients, as shown in Fig. 6. The ability to estimate DPT sensitivity coefficients to this resolution is a significant strength of the DPT methodology; obtaining sensitivity coefficients to this degree of detailed energy resolution allows for improved understanding of which evaluated nuclear data contribute the highest degree of uncertainty to depletion responses and can also be used to develop neutron energy filters to optimize the production of specific radioisotopes. The curse of dimensionality makes estimating energy-dependent sensitivity coefficients very computationally expensive for stochastic uncertainty analysis and direct perturbation methods.

Figure 6 shows that the fission of ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ suppresses the production of ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$, while ${\mathit{\hspace{0.17em}}}^{238}\mathrm{U}$ absorption facilitates the production of ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$, particularly at resonance energies. These results were unsurprising, as ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$ clearly relies on ${\mathit{\hspace{0.17em}}}^{238}\mathrm{U}$ neutron absorption and power normalization due to ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ fission lowering the effective ${\mathit{\hspace{0.17em}}}^{238}\mathrm{U}$ absorption rate. However, this methodology will provide valuable insight into less intuitive radioisotope production mechanisms (e.g., ${\mathit{\hspace{0.17em}}}^{252}\mathit{Cf}$ production). By verifying that integral quantities (e.g., energy-independent cross-section sensitivity) are correct, it is possible to have confidence in the accuracy of the components.

III.B.2. Nuclide Inventory Sensitivity to Initial Nuclide Density

This section details the sensitivity of time-dependent nuclide inventories to initial nuclide densities. This section examines the sensitivity of all nuclide number densities to the initial inventory of ${\mathit{\hspace{0.17em}}}^{238}\mathrm{U}$, ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$, and ${\mathit{\hspace{0.17em}}}^{234}\mathrm{U}$, the only nuclides with significant initial concentrations. These number density sensitivities are unconstrained and thus describe the effect of perturbing a nuclide’s number density without holding the material density constant. The sensitivity coefficients start at zero for nuclides that are not the perturbation target, as a nuclide’s initial inventory is presumed to be independent of other nuclides.

The number density sensitivity coefficients demonstrated strong agreement with the direct perturbation results. The observed discrepancy for the cross-section sensitivities was not present except for the ${\mathit{\hspace{0.17em}}}^{157}\mathrm{Gd}$ sensitivities, which may perhaps be contaminated by stochastic noise.

Shown in Fig. 7, the direct effects for the ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$ production chain exhibited good agreement, although the power normalization term provided a minor adjustment that provided noticeable correction for the sensitivity coefficients. The indirect effect impact on ${\mathit{\hspace{0.17em}}}^{234}\mathrm{U}$ and ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ also showed good agreement, as did most of the fission products. The power normalization term played a pivotal role in the ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ number density sensitivity, and in the sensitivity of many of the ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ fission products; several of these sensitivities changed signs due to the power normalization term. Here, the increased production of ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$, resulting from the increased ${\mathit{\hspace{0.17em}}}^{238}\mathrm{U}$ number density, resulted in more fissions from nuclides that were not ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$. As a result, the consumption of ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ was reduced to maintain the desired power in the simulation. This highlights the importance of considering how the power normalization affects the scaling of the flux.

Shown in Fig. 8, the nuclide number density sensitivities to the ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ initial number density are very strongly impacted by the power normalization term. The power normalization term essentially corrects all sensitivity coefficients, often changing their sign to properly match with the direct perturbation result. Because the power in the reactor is set, the increased number density of ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ does not result in more ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ being consumed and transmuted (since the amount of power produced by the system must remain the same). The increase in ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ number density suppresses the ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$ production chain, as more of the flux is occupied by fissions in ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ resulting in less flux being available for ${\mathit{\hspace{0.17em}}}^{238}\mathrm{U}$ to absorb and produce ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$. In this case, the power normalization can be understood as the ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$ acting as a shield for other nuclides, suppressing the amount of flux they receive.

Some small discrepancies were observed in the ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$ activation chain sensitivities, but it was unclear if these discrepancies were due to discrepancies in the ${\mathit{\hspace{0.17em}}}^{239}\mathrm{U}$ sensitivities or perhaps unresolved uncertainty in the GPT reaction rate sensitivity tally estimates.

As was observed for the ${\mathit{\hspace{0.17em}}}^{234}\mathrm{U}$ cross-section sensitivities, the ${\mathit{\hspace{0.17em}}}^{234}\mathrm{U}$ initial number density sensitivities in Fig. 9 were strongly contaminated by noise. Upon inspection, the direct perturbations, outside of ${\mathit{\hspace{0.17em}}}^{234}\mathrm{U}$ and ${\mathit{\hspace{0.17em}}}^{235}\mathrm{U}$, were washed out completely by stochastic noise. For the most part, the sensitivity coefficients oscillated around zero. Because ${\mathit{\hspace{0.17em}}}^{234}\mathrm{U}$ makes up a small fraction of the depletable inventory, the low sensitivity coefficients for the nuclide inventory were expected and understandable; however, the DPT and direct perturbation sensitivities seemed to agree well within the bounds of this stochastic noise.

III.B.3. Nuclide Inventory Sensitivity to Radioisotope Decay Constants

This section details the sensitivity of the time-dependent nuclide inventory to specific radionuclide decay constants. In this specific instance, the decay constants do not allow for any branching (i.e., there is only one decay path). This study examined the nuclide inventory sensitivities to the ${\mathit{\hspace{0.17em}}}^{239}\mathrm{U}$, ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Np}$, ${\mathit{\hspace{0.17em}}}^{135}\mathrm{I}$, and ${\mathit{\hspace{0.17em}}}^{135}\mathrm{Xe}$ decay constants, which are presented in this section. All other nuclides were presumed to not decay, as their half-life is long compared to the length of the burn time.

The DPT ${\mathit{\hspace{0.17em}}}^{239}\mathrm{U}$ decay constant sensitivities generally exhibited generally good agreement with the reference direct perturbation results. There was a minor discrepancy in the sensitivity coefficient for ${\mathit{\hspace{0.17em}}}^{239}\mathrm{U}$, which is believed to be due to the aforementioned DPT sensitivity coefficient integration scheme; in contrast, ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Np}$ and ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$ both exhibited fantastic agreement with the direct perturbation results.

The ${\mathit{\hspace{0.17em}}}^{239}\mathrm{U}$ direct perturbation results were linear for all nuclides at the start of the time steps (Fig. 10). The direct perturbation results were large for ${\mathit{\hspace{0.17em}}}^{239}\mathrm{U}$, ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Np}$ and ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$, but were very small for the rest of the nuclide inventory. As the simulation progressed, all sensitivity coefficients (aside from ${\mathit{\hspace{0.17em}}}^{239}\mathrm{U}$) degraded into a noise regime by the end of the simulation. By the second time step, 0.4 days, the rest of the nuclide inventory had degraded into a noise regime.

At the 5-day mark, ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Np}$ and ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$ were still somewhat linear. At the 10-day mark, ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Np}$ was nonlinear while ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$ was significantly degraded but still had a potential trend, but by the 20-day mark ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$ was also nonlinear. This was despite the fact that the fluctuation order was ${10}^{-4}$. This was understandable upon inspecting Fig. 10, as the ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Np}$ and ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$ sensitivity coefficients trended toward zero as time went on, while the rest of the inventory appeared to be independent of the ${\mathit{\hspace{0.17em}}}^{239}\mathrm{U}$ decay constant.

The DPT ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Np}$ decay constant sensitivities in Fig. 11 generally agreed well with the reference direct perturbation results; however, there was a trend in the fission product sensitivities that did not appear to be well captured.

Fig. 11. Sensitivity of the nuclide inventories to the ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Np}$ decay constant.

Similar to ${\mathit{\hspace{0.17em}}}^{239}\mathrm{U}$, the ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Np}$ decay constant sensitivities started out linear as well. For most of the nuclide inventory, the results were small ($\frac{\mathit{\delta R}}{R}{}_{}\pm {10}^{-8}$), while as expected, ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Np}$ and ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$ were linear and large. By the second time step, 0.4 days, the rest of the nuclide inventory had degraded into a noise regime. By the end of the simulation, only ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Np}$ and ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$ were linear while the rest of the nuclide inventory had degraded into the noise regime. In comparison to ${\mathit{\hspace{0.17em}}}^{239}\mathrm{U}$, there did not appear to be a well-defined point where the sensitivities were expected to degrade, though extrapolating based on the trend in ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$ gives a prediction that perhaps by 60 days that sensitivity coefficient would also tend toward zero and become noise.

In Fig. 12 ${\mathit{\hspace{0.17em}}}^{135}\mathrm{I}$ decay constant sensitivities showed generally good agreement; however, there was an apparent issue that resulted in null sensitivity coefficients for the uranium isotopes and the ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$ buildup chain.

Fig. 12. Sensitivity of the nuclide inventories to the ${\mathit{\hspace{0.17em}}}^{135}\mathrm{I}$ decay constant.

The direct perturbation sensitivity estimates were well resolved for ${\mathit{\hspace{0.17em}}}^{135}\mathrm{I}$, ${\mathit{\hspace{0.17em}}}^{135}\mathrm{Cs}$, ${\mathit{\hspace{0.17em}}}^{135}\mathrm{Xe}$, and ${\mathit{\hspace{0.17em}}}^{136}\mathrm{Xe}$ during the first time step. The rest of the nuclide inventory was very insensitive ($\pm {10}^{-16}$ magnitude) to the ${\mathit{\hspace{0.17em}}}^{135}\mathrm{I}$ decay constant. In the final time step, the ${\mathit{\hspace{0.17em}}}^{135}\mathrm{I}$ and ${\mathit{\hspace{0.17em}}}^{135}\mathrm{Cs}$ direct perturbation sensitivities remained well resolved, but the ${\mathit{\hspace{0.17em}}}^{136}\mathrm{Xe}$ sensitivity estimates were contaminated by stochastic noise, as were the ${\mathit{\hspace{0.17em}}}^{135}\mathrm{Xe}$ sensitivities, which oscillated around zero. However, Fig. 12 indicates that the prevalence of stochastic noise for xenon was unsurprising, as the sensitivity coefficients were tending toward zero as time progressed. The rest of the nuclide inventory remained noisy as before, with fluctuations up to $\pm {10}^{-4}$.

As observed in Fig. 13 for ${\mathit{\hspace{0.17em}}}^{135}\mathrm{I}$, the ${\mathit{\hspace{0.17em}}}^{135}\mathrm{Xe}$ decay constant sensitivities generally exhibited good agreement; however, there was an apparent issue that resulted in null sensitivity coefficients for the uranium isotopes and the ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$ buildup chain.

Fig. 13. Sensitivity of the nuclide inventories to the ${\mathit{\hspace{0.17em}}}^{135}\mathrm{Xe}$ decay constant.

At the first time step, the direct perturbations were linear for ${\mathit{\hspace{0.17em}}}^{135}\mathrm{Cs}$, ${\mathit{\hspace{0.17em}}}^{135}\mathrm{Xe}$, and ${\mathit{\hspace{0.17em}}}^{136}\mathrm{Xe}$. The rest of the nuclide inventory was noise ($\pm {10}^{-16}$ magnitude). The ${\mathit{\hspace{0.17em}}}^{135}\mathrm{Cs}$, ${\mathit{\hspace{0.17em}}}^{135}\mathrm{Xe}$, and ${\mathit{\hspace{0.17em}}}^{136}\mathrm{Xe}$ remained linear throughout the problem, and the rest of the noisy inventory did not improve, though the noise oscillated at magnitudes up to $\pm {10}^{-4}$.

Figures 14 and 15 demonstrate the lack of sufficient resolution in the direct perturbations for decay constant sensitivities that indirectly impact decay and the activation chains. This appeared to be a trend for the decay constant sensitivities presented in this paper. This limited resolution is believed to be a product of stochastic noise, as the decay constants had no impact on the initial flux profile and only indirectly impacted the flux profile for the later time steps. In comparison to cross sections and number densities, the decay constants did not directly participate in establishing the flux profile or contribute to any power normalization. Thus, these sensitivity coefficients were expected to be significantly affected by stochastic noise.

Fig. 14. Direct perturbation data for ${\mathit{\hspace{0.17em}}}^{239}\mathrm{U}$ after the first two depletion cycles.

Fig. 15. Direct perturbation data for ${\mathit{\hspace{0.17em}}}^{239}\mathrm{U}$ at the later time steps.

IV. CONCLUSIONS

This paper presented proof of principle for a depletion sensitivity coefficient methodology in continuous-energy OpenMC Monte Carlo simulations. Recent advances in GPT reaction rate sensitivity analysis for Monte Carlo simulations were combined with previous DPT to develop fully coupled depletion sensitivity coefficients in Monte Carlo depletion simulations. The DPT methodology developed in this work produced sensitivity coefficient estimates that were consistent with direct perturbation sensitivity coefficients for cross section and number density sensitivities with few notable exceptions. The ${\mathit{\hspace{0.17em}}}^{239}\mathrm{Pu}$ buildup chain tended to have issues with ${\mathit{\hspace{0.17em}}}^{239}\mathrm{U}$ (see Figures 3 and 4), perhaps due to numerical convergence, which caused prominent disagreement with the reference direct perturbation results.

Through use of a sample problem containing various transmutation paths and fission products with a range of half-lives and absorption cross sections, this work successfully verified the calculation of direct effect, indirect effect, and power normalization sensitivity components. The power normalization term was found to play a crucial role in calculating several sensitivity coefficients, most notably the sensitivity of nuclide inventories to the initial number density. Discrepancies with the decay constant sensitivities were attributed to stochastic noise and are currently being investigated for further systems.

The adjoint-weighted depletion sensitivity analysis methods developed by this work offer the potential to estimate depletion sensitivity coefficients to a degree of resolution that is generally not achievable using stochastic uncertainty analysis methods. In future work, these high-resolution sensitivity coefficients may be combined with covariance data via the Sandwich Equation to quantify the impact of nuclear data uncertainty on depletion simulation results.

Recent improvements to nuclear data covariance data libraries, such as the fission yield covariance data library recently developed by Lovell et al., will prove essential in enabling this Sandwich Equation–based uncertainty analysis.^[8] The availability of high-resolution depletion sensitivity coefficients and increasingly complete covariance data libraries enables a variety of powerful and rapid uncertainty analysis techniques, including similarity assessment for designing optimal benchmark experiment measurements^[9] and data assimilation to calibrate nuclear data, which are difficult to measure directly.^[10]

Previous work indicates that CLUTCH-based GPT sensitivity methods scale well as the number of responses increases; the memory footprint and run time of those simulations have been observed to increase approximately by 1% to 2% per additional response.^[4] Thus, while computing a large array of responses for a detailed depletion simulation may incur a significant computational cost, the DPT methodology may be preferable to the very large array of direct perturbation calculations required for a stochastic uncertainty analysis.

Disclosure Statement

The authors report there are no competing interests to declare.

Additional information

Funding

This research was performed using funding received from the U.S. Department of Energy Office of Nuclear Energy’s Nuclear Energy University Programs under contract number DE-NE0008890.