• $\epsilon_{\text{eff}}$: specific energy per fission (MeV).

The $P_{dE}(t, T)$, the decay heat power from neutron capture in fission products except for $^{133}$Cs is determined from the decay heat of the fission products $P'_d$ and reactor independent factors tabulated as a function of decay time $H(t)$:

$$P_{dE}(t, T) = P'_{d}(t, T) \cdot H(t). \quad (41)$$

The contribution to decay heat power from actinides is computed by two separate sources. The decay heat power from the decay of $^{239}$U and $^{239}$Np, $P_{dHE}(t, T)$, and the decay heat power from actinides other than $^{239}$U and $^{239}$Np, $P_{dA}(t, T)$:

$$P_{dHE}(t, T) = \frac{P_{\max}}{Q_{\text{eff}}} [F_{U}(t, T) + F_{N_{p}}(t, T)] \quad (42)$$

with

$$F_{U}(t, T) = E_{U} \cdot R \cdot (1 - e^{-\lambda_{U}T}) \cdot e^{-\lambda_{U}t} \quad (43)$$

$$F_{N_{p}}(t, T) = E_{N_{p}} \cdot R \left[ \frac{\lambda_{U}}{\lambda_{U} - \lambda_{N_{p}}} (1 - e^{-\lambda_{N_{p}}T}) e^{-\lambda_{N_{p}}t} - \frac{\lambda_{N_{p}}}{\lambda_{U} - \lambda_{N_{p}}} (1 - e^{-\lambda_{U}T} e^{-\lambda_{U}t}) \right] \quad (44)$$

and

• $E_U$: average energy from decay of $^{239}$U (MeV),

• $E_{N_p}$: average energy from decay of $^{239}$Np (MeV),

• R: $^{239}$U production from capture in $^{238}$U at the end of irradiation (at/s),

• $\lambda_U$: decay constant of $^{239}$U ($s^{-1}$), and

• $\lambda_{N_p}$: decay constant of $^{239}$Np ($s^{-1}$).

The production of $^{239}$U (R-factor) must be computed and provided by the user. The contribution to decay heat of all other actinides is accounted for by a tabulated correction factor, and with dependence on burnup and initial enrichment, applied to the uncorrected contribution for neutron capture decay heat of fission products by the following equation:

$$P_{dA}(t, T) = P'_{d}(t, T) \cdot \left( 10^{-4} \frac{A(t)B^{\alpha(t)}}{\epsilon_{\text{eff}}} \right) \quad (45)$$

with $A(t)$ a tabulated correction factor, $B$ the burnup, $\epsilon_{\text{eff}}$ the effective $^{235}$U enrichment (wt.% of initial $^{235}$U), and

$$\alpha(t) = 1.5 \cdot e^{-\frac{t}{10^9}}. \quad (46)$$

The total uncertainty is determined from the uncertainty in the fission product decay heat power without neutron capture in fission products and the uncertainty in the operating power. The other terms in the total decay heat, neutron capture in fission products and actinides, are defined to provide conservative overestimates of their contributions to the decay heat power, and the uncertainties in these terms are therefore not included in the total uncertainty:

$$\left[ \frac{\Delta P_d}{P_d} \right]^2 = \left[ \frac{\Delta P'_d}{P'_d} \right]^2 + \left[ \frac{\Delta P}{P} \right]^2 \quad (47)$$

with $P, \Delta P$ the reactor power and reactor power uncertainty. The decay heat uncertainty of fission products $\Delta P'_d$ is defined by:

$$\left(\Delta P'_{d}(t, T)\right)^{2} = \sum_{i=1}^{4} \left[ \left(\frac{\Delta Q_{i}}{Q}\right)^{2} + \frac{1}{Q_{i}P'_{di}} \left( \sum_{k=1}^{N} P_{i,k} \sum_{j=1}^{23} \frac{\beta_{i,j}}{\lambda_{i,j}} (1 - e^{-\lambda_{i,j}T_{k}}) \right)^{2} \left(\cdot e^{-\lambda_{i,j}T_{k}}\right)^{2} \right]. \quad (48)$$

The four sets of 23 coefficients $\beta_{i,j}$ used in the calculations of the uncertainty are provided by the standard for the thermal fissions of $^{235}$U, 239Pu, and $^{241}$Pu, and the fast fissions of $^{238}$U.

### 6.3.6 ISO 10645 (2022)

The standard provides procedures for calculating the decay heat power of non-recycled nuclear fuel irradiated in light water reactors. The methods apply to light water reactors (pressurized water and boiling water reactors) loaded with a nuclear fuel mixture consisting of $^{235}$U and $^{238}$U. By its implementation, the methods applied in this standard are very similar to those in ANSI/ANS-5.1-2014.

#### History

The first released in 1992 is almost identical to DIN 25463 1990:

• applicable to PWR and BWR,

• for fission products: 24 exponential terms for $^{235}$U, $^{238}$U, $^{239}$Pu and $^{241}$Pu,

• with explicit formula with approximated neutron flux for $^{134}$Cs production,

• based on a G-factor for neutron capture in fission products other than $^{133}$Cs,

• explicit description for both actinides $^{239}$U and $^{239}$Np,

• correction factor to account for actinides other than $^{239}$U and $^{239}$Np,

• initial $^{235}$U enrichment up to 4.1 wt.%,

• and a burnup limit of 52 MWd/kgU.

The 2022 implementation has the following characteristics:

• Fission products: 23 exponential terms for $^{235}$U, $^{238}$U, $^{239}$Pu and $^{241}$Pu from ANSI/ANS-5.1-2014,

• initial $^{235}$U enrichment increased up to 5 wt.%,

• Burnup range is extended to 62 MWd/kgU,

• Decay time less than $10^9$ seconds,

• New procedure for neutron capture in $^{133}$Cs, and

• Updated formulations and nuclear data.

#### Applicability and limitations

The following limitations and range of applications need to be considered:

• The methods in the standard are applicable to UO$_2$ depleted in either PWR or BWR.
* The standard is not applicable to MOX.

* The standard is applicable for cooling times from 0 to $10^9$ seconds.

* The standard yields conservative results for:

- Initial enrichment $a_0$, between 1.9 wt.% $\le a_0 \le 5$ wt.%,

- Burnup $B_f \le 12.5$ a$_0$ MWd/kgU,

- Power density $S \ge 5$ a$_0$ kW/kgU.

Decay heat power from delayed neutron-induced fission and activation of structural materials are not included. The user is advised to evaluate those and appropriately include them in any analysis of decay heat.

## Methodology

The total decay heat $P_N$ after a decay time $t$ which follows a power period $T$, is computed as a sum of the following separate sources:

$$P_N(t, T) = P_S(t, T) + P_B(t, T) + P_A(t, T) + P_{Cs}(t, T) + P_E(t, T), \quad (49)$$

where

* $P_S$: the decay heat power from fission products,

* $P_B$: the contribution of $^{239}$U and $^{239}$Np,

* $P_A$: the contribution of other nuclides,

* $P_{Cs}$: the contribution of $^{134}$Cs, and

* $P_E$: the decay heat power from neutron capture in other fission products.

The decay heat from fission products is computed by approximating a pulse power function by a sum of 23 exponential terms multiplied by the reaction rates of the four main fissile nuclides and integrated over the power intervals (k). The power function, the innermost sum in the equation below, represents the decay heat power per fission following an instantaneous pulse of a significant number of fission events:

$$P_S(t, T) = \sum_{i=1}^{4} \sum_{k=1}^{m} \frac{P_{ik}}{Q_i} \sum_{j=1}^{23} \frac{\alpha_{ij}}{\lambda_{ij}} (1 - \exp^{-\lambda_{ij}T_k}) \times \exp^{-\lambda_{ij}t_k}, \quad (50)$$

with the following factors:

* $\lambda_{ij}$, $\alpha_{ij}$: the coefficients in the 23-term exponential fit for the fission products decay heat power,

* $P_{ij}$: the fission rate fraction of nuclide i ($^{235}$U, $^{238}$U, $^{239}$Pu and $^{241}$Pu) in step k,

* $Q_i$: the specific energy per nuclide i ($^{235}$U, $^{238}$U, $^{239}$Pu and $^{241}$Pu), and

* t, T: the time on power and cooling time for step k.

The method and the data are those developed for ANSI/ANS-5.1-2014. The $^{239}$U and $^{239}$Np contributions to decay heat are treated separately from the rest of the actinides:

$$P_B(t, T) = \sum_{k=1}^{m} \frac{P_k}{Q} [F_U(t_k, T_k) + F_{Np}(t_k, T_k)]. \quad (51)$$

The term $P_k/Q$ is the total fission rate in time interval $k$, and the terms $F_U(t_k, T_k)$ and $F_{Np}(t_k, T_k)$ are computed by:

$$F_U(t_k, T_k) = E_U \times R (1 - \exp^{-\lambda_U T_k}) \times \exp^{-\lambda_U t_k} \quad (52)$$

$$F_{Np}(t_k, T_k) = E_{Np} \times R \left[ \frac{\lambda_U}{\lambda_U - \lambda_{Np}} (1 - \exp^{-\lambda_{Np} T_k}) \times \exp^{-\lambda_{Np} t_k} - \frac{\lambda_{Np}}{\lambda_U - \lambda_{Np}} (1 - \exp^{-\lambda_U T_k}) \times \exp^{-\lambda_U t_k} \right] \quad (53)$$

with

* $E_U$, $\lambda_U$: the mean decay energy and decay constant of $^{239}$U,

* $E_{Np}$, $\lambda_{Np}$: the mean decay energy and decay constant of $^{239}$Np, and

* R: the ratio of capture in $^{239}$U to total fission rate at the end of operating time.

An approximation for R is provided as a function of initial enrichment $a_0$ and final burnup $B_f$:

$$R = 0.974 a_0^{-0.504} + B_f (0.00883 - a_0 0.000726). \quad (54)$$

For initial enrichments between 1.9 wt.% and 5 wt.%, the approximation yields conservative results. The decay heat from decay of actinides other than $^{239}$U and $^{239}$Np, is determined from the decay heat power from fission products multiplied by tabulated vs decay time factor $A(t)$:

$$P_A(t, T) = A(t) \times P_S(t, T). \quad (55)$$

It stated that it will “yield conservatively high results”, provided the following conditions are fulfilled:

* Initial enrichment, 1.9 wt.% $\le a_0 \le 5.0$ wt.%,

* Burnup $\le 12.5 \times a_0$ MWd/kgU, and

* Power density $\ge 5.0 \times a_0$ kW/kgU.

The decay heat due to radioactive decay of $^{134}$Cs generated by neutron capture in $^{133}$Cs is given by the following equation:

$$P_{Cs}(t, T) = \frac{P}{Q} \lambda_4 E_{Cs} y \left[ \frac{1 - \exp^{-(\lambda_4 + \sigma_4 \Phi_{eff}) T_{eff}}}{\lambda_4 + \sigma_4 \Phi_{eff}} + \frac{\exp^{-(\sigma_3 \Phi_{eff} T_{eff})} - \exp^{-(\lambda_4 + \sigma_4 \Phi_{eff}) T_{eff}}}{\sigma_3 \Phi_{eff} - (\lambda_4 + \sigma_4 \Phi_{eff})} \right] \times \exp^{-\lambda_4 t} \quad (56)$$

with

$$\frac{P}{Q} = \frac{1}{T_{eff}} \sum_{k=1}^{m} \sum_{i=1}^{4} \frac{P_{ik}}{Q_i} T_k \quad (57)$$

and

* y: the mean cumulative $^{133}$Cs yield per fission,

* $\lambda_4$ the decay time constant for $^{134}$Cs (in sec$^{-1}$),

* $\sigma_3$ the spectrum averaged (n,γ) cross section in $^{133}$Cs (in barn),
<table><caption>**Table 16.** Measured and calculated decay heat value for the Ringhals-3 assembly 0E2.</caption><thead><tr><th></th><th>5823</th><th>6389</th><th>6390</th><th>7826</th><th>7837</th><th>7970</th></tr></thead><tbody><tr><td>Cooling time (days)</td><td>587.9 ± 7.0</td><td>566.0 ± 6.9</td><td>567.7 ± 6.9</td><td>522.4 ± 6.6</td><td>525.6 ± 6.6</td><td>520.1 ± 6.6</td></tr><tr><td>Measured DH (Watts)</td><td>587.9</td><td>568.4</td><td>568.3</td><td>535.7</td><td>525.4</td><td>521.8</td></tr><tr><td>Code 1 (Watts)</td><td>590.7</td><td>571.6</td><td>571.5</td><td>528.9</td><td>528.6</td><td>525.0</td></tr><tr><td>Code 2</td><td>579.8</td><td>561.0</td><td>561.0</td><td>519.4</td><td>556.6</td><td>515.6</td></tr><tr><td>Code 3</td><td>591.1</td><td>572.1</td><td>572.1</td><td>529.9</td><td>529.9</td><td>526.2</td></tr><tr><td>Code 4</td><td>584.2</td><td>564.0</td><td>564.0</td><td>520.1</td><td>519.8</td><td>516.0</td></tr><tr><td colspan="7" style="text-align: center;">History 1</td></tr><tr><td>ANSI/ANS-5.1-2014</td><td>738.7 ± 22.2</td><td>716.4</td><td>716.4</td><td>671.2</td><td>670.9</td><td>667.2</td></tr><tr><td>DIN 25463 (2014)</td><td>573.1 ± 22.5</td><td>554.7</td><td>554.6</td><td>514.2</td><td>513.9</td><td>510.5</td></tr><tr><td>ISO 10645:2022</td><td>685.7 ± 20.6</td><td>662.5</td><td>662.4</td><td>610.9</td><td>610.5</td><td>606.1</td></tr><tr><td>RG2 3.54 Rev2 (2018)</td><td>625.9</td><td>602.6</td><td>602.6</td><td>553.7</td><td>553.4</td><td>549.4</td></tr><tr><td colspan="7" style="text-align: center;">History 2</td></tr><tr><td>ANSI/ANS-5.1-2014</td><td>685.9 ± 23.8</td><td>665.2</td><td>665.2</td><td>623.3</td><td>623.0</td><td>619.6</td></tr><tr><td>DIN 25463 (2014)</td><td>576.5 ± 22.8</td><td>557.9</td><td>557.8</td><td>517.0</td><td>516.7</td><td>513.2</td></tr><tr><td>ISO 10645:2022</td><td>635.8 ± 22.7</td><td>614.2</td><td>614.2</td><td>566.6</td><td>566.2</td><td>562.2</td></tr><tr><td>RG2 3.54 Rev2 (2018)</td><td>605.8</td><td>583.3</td><td>583.2</td><td>536.1</td><td>535.8</td><td>532.0</td></tr></tbody></table>

* $\sigma_4$ the spectrum averaged absorption cross section in $^{134}$Cs (in barn),

* $E_{Cs}$ the recoverable energy per decay for $^{134}$Cs (in MeV),

* $T_{eff} = \sum_{k=1}^{m} T_k$ the effective irradiation time excluding outage periods,

* $Q_i$ the specific energy per fission for the four fissionable nuclides (in MeV),

* $\Phi_{eff} = \sum_{k=1}^{m} \Phi_k T_k$ the effective neutron flux (given in n/cm$^2$/s).

The value of $\Phi_k$ is given by $\Phi_k = S_k/a_{eff} \times 2.58 \times 10^{13}$ with $S_k$ the power density (in kW/kgU) and $a_{eff} = a_0/2 + 1.0$, being the effective enrichment (in wt.%).

The approximate neutron flux yields conservative $P_{Cs}$ for typical LWR enrichments and burnups by up to 5%. For burnup values less than 25 MWd/kgU, the approximate solution overestimates $P_{Cs}$ by up to 15%.

The $P_A$, the decay heat power from neutron capture in fission products except for $^{133}$Cs is determined from the decay heat of the fission products $P_s$ and reactor independent factors tabulated as a function of decay time $H(t)$:

$$P_E(t, T) = H(t) \times P_s(t, T). \qquad (58)$$

The uncertainty of calculated decay heat power is defined by:

$$\Delta P_N(t, T) = n \sqrt{(P_s(t, T))^2 + \left(P_N(t, T) \frac{\Delta P}{P}\right)^2} \quad (59)$$

with

* $n$ is the multiple of standard deviation for the chosen confidence interval,

* $\Delta P_s$ the uncertainty in decay heat from fission products, and

* $\Delta P/P$ is the uncertainty of relative thermal power.

“The other contributions to the decay heat power PB, PA, PCs, and PE shall be determined conservatively and therefore do not enter into the calculation of the error bandwidth. Using the approximate methods of this document for these contributions results in conservative estimates of the total decay heat power. Alternative methods, such as those based on comparisons of code predictions and isotopic measurements for the main nuclides contributing to each of these decay heat terms, shall be specified, and justified by the user”.

### 6.3.7 Example of application

The standards described in the previous sections can be applied to a large number of cases (see for instance Ref. [55]). They can nevertheless lead to different estimations of decay heat and it is of interest to present an example for a specific SNF case. To facilitate the comparison of decay heat values from standards, some simulation codes automatically include standard decay heat values during depletion calculations (e.g. the SNF code). Additionally, to help compare values within this community, Studsvik made available for all SG12 participants a simplified version of the SNF code, allowing them to calculate, with the same input definition, decay heat values for four standards, as presented in Table 16 for the example of the assembly 0E2. At the time of writing, this code is available from the SG12 webpage [8].

A specific example is presented in the following, extracted from the Clab data presented in the 2006 report [76]. The assembly named 0E2 corresponds to a PWR 17×17 assembly with 3.10% enriched $^{235}$U fuel, and its burnup is reported at 41.6 MWd/kgU. The full description can be found in the report, together with the value of one measurement, at a cooling of 5823 days (almost 16 years). Five other measurements were later reported in reference [79],
<figure><img src="image_30.png" /><figcaption>Fig. 27. Top: comparison of the data presented in Table 16 (dashed lines are the values from the code 1 to 5). Bottom: Comparison for the standard only, with both irradiation histories.</figcaption></figure>

for cooling time between 17 to 22 years. Measured and calculated values from the standards are reported in Table 16.

For these cooling periods, the important contributors to the decay heat are fission products, i.e. ¹³⁷mBa/¹³⁷Cs and ⁹⁰Y/⁹⁰Sr (almost 65% in total), and also ²⁴⁴Cm (10%), ²³⁸Pu (10%), and ²⁴¹Am (another 10%). The remaining parts originate in ²⁴⁰Pu and ¹⁵⁴Eu. Therefore, given the similar cooling periods, these 6 measurements are helping to validate the same nuclides.

Table 16 presents values from dedicated best estimate calculations (called “code 1” to “code 5” to keep anonymity), partly performed by different laboratories using different nuclear data libraries. Such calculations are very close to the measured values, and well within the reported experimental one-standard deviation. Regarding the application of the standards, two sets of histories are used: one simplified and one more detailed. Both sets can be accepted as corresponding to a good representation of the realistic irradiation history. If one first considers a single history (number 1 or 2), it can be observed that the values provided by the four standards present a large spread, larger than the one from the best-estimate calculations. As expected, three standards are conservative, and one (the DIN standard) is closer to the experimental values (see Fig. 27 top). Differences between standards are still noticeable, especially given the experimental and standard uncertainties. This example illustrates that standard values are not systematically coherent between themselves and are generally conservative (except for the DIN standard, but the conservatism is not systematic, as presented in Ref. [55]). The impact of the irradiation history (and the underlying assumptions) is also not negligible, changing the calculated values by more than 5%,

as presented in Figure 27 bottom. A dedicated discussion is out of the scope of this work, but the following remarks can be made for this example:

1. Standard values can be sensibly different compared to the ones from best-estimate codes.

2. Standard values are in general more conservative.

3. Standard values can disagree between themselves.

4. As for best-estimate calculations, standard values are sensitive to assumptions made in simplified irradiation histories.

## 6.4 Summation and best estimate methods

### 6.4.1 Summation method

Decay heat (rate) is a unique characteristic for any given nuclide composition of nuclear fuel and in general of radioactive matter. The decay rate (or activity) of a radioactive nuclide i is the product of its number density $N_i$ and its decay constant $\lambda_i$. Multiplying the decay rate $\lambda_i \cdot N_i$ by the recoverable energy $e_i$ released as radiation during the decay process gives the decay heat rate contributed by that nuclide,

$$H_i = e_i \cdot \lambda_i \cdot N_i. \qquad (60)$$

Then, the total decay heat rate H is given by the sum of each individual contribution

$$H = \sum_{i} H_{i}, \qquad (61)$$

hence the name summation method. Among the different methods to model the decay heat of radioactive nuclear fuel after shutdown, the summation method is considered the most general and comprehensive. Despite being widely accepted and adopted in the earliest nuclear applications especially, for crude predictions of the order of magnitudes in short-term decay heat analysis, simple semi-empirical formulations such as those proposed in the previous section have been generally nowadays superseded by the summation method [144]. The completeness of this method stems from explicitly considering the contributions of all nuclides involved in the decay process, including actinides in the nuclear fuel, fission, and activation products. Despite its simplicity, an accurate knowledge of a large number of nuclear parameters (decay constants, beta-electron, and gamma-ray energies) is required.

Given the radioactive nature of nuclear fuel, its decay heat rate inevitably varies with time according to the nuclide content evolution in the system. For reactor spent fuel it decreases following the time after shutdown, as the majority of the short-lived fission products decay into longer-lived ones. It follows that the accuracy of the summation method strictly depends on the correct solution of the Bateman equations. This step is generally considered more crucial than the summation calculation, as the complexity of the process increases when the fuel is placed in an irradiating environment.
The radioactive decay law formulated in 1902 by Ernest Rutherford and Frederick Soddy [145] indicates that the density variation $dN(t)/dt$ of radioactive matter is proportional to its activity. Here, $N(t)$ is the nuclide field vector containing each $i^{th}$ nuclide number density. In his publication at the beginning of the last century [146], Henry Bateman reported the analytical solution of this set of equations considering only natural decay and under the assumption that only the first nuclide in the decay chain had non-zero density.

Still, the concept can be generalized for radioactive matter placed in a particle-irradiating environment by adopting effective reaction rates (or simply reaction rates) that take into account the transmutation caused by particle-nucleus interactions. Each $i^{th}$ ordinary differential equation (ODE) describes the time-balance of a nuclide $i$ as

$$ \frac{dN_i(t)}{dt} = \sum_j \left[ (b_{j \to i} \lambda_j + \sum_p \phi_p \langle \sigma_{p,j \to i} \rangle) \cdot N_j(t) - (\lambda_i + \sum_p \phi_p \langle \sigma_{p,i} \rangle) \cdot N_i(t) \right] \quad (62) $$

where

* $N_i, N_j$ are the number densities of nuclides $i$ and $j$, respectively,

* $\lambda_i, \lambda_j$ are the decay constants of nuclides $i$ and $j$, respectively,

* $b_{j \to i}$ is the decay branching ratio for a nuclide $j$ to produce nuclide $i$,

* $\langle \sigma_{p,i} \rangle$ is the spectrum-averaged Doppler-broadened disappearance cross-section for nuclide $i$ induced by particle $p$,

* $\langle \sigma_{p,j \to i} \rangle$ is the spectrum-averaged Doppler-broadened cross section for nuclide $i$ induced by particle $p$ and that generates nuclide $j$,

* $\phi_p$ is the energy-integral of the flux for particle $p$, being either neutrons, protons, or other particles$^5$.

The spectrum-averaged cross-sections are

$$ \langle \sigma_{p,i} \rangle = \frac{\int \sigma_{p,i}(E) \varphi_p(E) dE}{\int \varphi_p(E) dE} \quad (63) $$

$$ \langle \sigma_{p,j \to i} \rangle = \frac{\sum_r \int \gamma_{p,j \to i,r} \sigma_{p,j,r}(E) \varphi_p(E) dE}{\int \varphi_p(E) dE} \quad (64) $$

where the energy distribution of the particle flux, i.e. the particle spectrum, is denoted by $\varphi_p(E)$. The production rate includes all reactions $r$ with cross-section $\sigma_{p,j,r}(E)$ with a non-zero yield $\gamma_{p,j \to i,r}(E)$ that generates the nuclide $i$ of interest. If nuclide $i$ is a direct fission product and nuclide $j$ is an actinide, $r$ represents the fission reaction and $\gamma_{p,j \to i,r}(E)$ corresponds to the independent fission yield. An independent fission yield is defined as the fraction of a specific nuclide produced directly from a single neutron-induced fission event of a given parent nuclide after the emission of prompt neutrons and gamma rays and before the radioactive decay of the fission fragments.

The particle flux field, normalized to the requested power production, follows the time-dependent neutron transport equation (see Ref. [147]) and is nonlinearly coupled to the calculation for the nuclide vector.

From the above description, the accuracy of the summation can depend on various approximations performed during the irradiation (see next sections), and also on the quality of nuclear data (see Sect. 6.4.3).

### 6.4.2 Case specific approximations

The transmutation or decay of fission products with large absorption cross sections or in general of strong neutron absorbers can significantly affect the spatial and energy distribution of the neutron flux on a time scale on the order of hours to days [148–150]. The buildup and depletion of actinides have an analogous effect on a longer time scale. Modifications of the irradiation conditions (e.g. reactor power, control rod position, soluble boron concentration) also alter the neutron flux. The dependence of the flux distribution on the nuclide composition is handled in burnup codes by decoupling the neutron transport calculation from the nuclide depletion problem, exploiting the slow variations of the neutron flux distribution (energy distribution as well as spatial distribution) in short time-scales. A quasi-static version of the space and energy-dependent burnup equations assumes static flux distribution that can be computed at the beginning of the time interval of interest by a deterministic or Monte Carlo-based particle transport solver. Then, this flux distribution and the resulting spectrum-averaged cross sections are used to integrate the depletion equations over the time interval for which the transition matrix can be assumed constant. This methodology reduces a nonlinear process into a series of equations that appear to be linearly coupled at given time instances.

There is a number of approximations, both physical and computational, that need to be applied when the burnup calculation is run. Computational approximations include discretization applied to space and time, i.e. depletion zone division and the number of depletion steps, as well as the burnup algorithm used for obtaining the coupled solution by iteration.

Physics approximations are related to the models used for the interaction physics. For deterministic transport codes this includes, for example, the methods applied for self-shielding, and for Monte Carlo codes the treatment of cross sections in the unresolved resonance range and the elastic scattering kernel. Some approximations, such as the definition of energy released per fission and the handling of fuel temperature distributions, result from simplifications in the computational model. The most significant uncertainties in burnup calculation, however, typically originate from the nuclear data libraries used in the calculation (see Sect. 6.4.3).

#### Spatial discretization

Because of the large costs in terms of computational power and memory necessary to run SNF depletion calculations for LWRs at the core level, a number of customary prac-

<footnote n="5">⁵ For most SNF applications nuclear reactions induced by particles other than neutrons are generally negligible.</footnote>
tices are nowadays accepted by the nuclear industry to reduce the computational running time, such as:

* Scaling down the core geometrical model to the assembly level with imposed boundary conditions,

* using a two-dimensional (2D) model,

* using a discrete number of depletion zones within the assembly,

* and using a discrete number of depletion zones within the pin.

The two last assumptions include grouping the depleted materials into homogeneous sub-regions, or depletion zones, inside which the reaction rates are assumed space-independent. If the spatial discretization is not carefully defined the local reaction rates may be under/over-estimated in regions where the neutron flux gradients are steep. This is a condition that typically occurs in presence of high-absorbing materials or at the core boundaries. Refining the spatial sub-division reduces the discretization bias, but in practice the total number of depletion zones is limited by the available computer capacity.

## Assembly model

To simulate the depletion of nuclear fuel in a fuel pellet, rod or fuel assembly (FA), computational models focus/extend up to the assembly boundaries. The rest of the core is defined by imposed boundary conditions. The first reason is computational time memory and resources: because of the large costs in terms of computational power and memory necessary to run SNF depletion calculations for LWRs at the core level. A second reason is that in many SNF-related benchmarks operational data are only provided on the reference assembly selected for the experimental campaign, whereas details on the adjacent assemblies are not documented [151,152]. The assembly shuffling is often not simulated when analyzed samples are taken from interior fuel rods as the impacts of neighboring assemblies on the nuclide contents in such rods were shown to be small after $^{148}$Nd normalization, see for instance reference [153]. On the contrary, it has been shown that for fuel rods at the periphery of an FA or if the reference FA is located at one edge of the core, modeling the FA surrounding environment has a significant impact on the accurate prediction of nuclide compositions as a function of burnup [154]. Under these conditions a more accurate representation of the rod environment is achieved by including the adjacent assemblies in the simulation model [155,156].

The down-scaling of the core model to a single assembly (or group of assemblies) is an acceptable approximation as long as the spatial boundary conditions are representative of the assembly position in the reactor core and evolve according to the fuel in-core reshuffling scheme. For instance, the use of reflective boundary conditions along all directions corresponds to an “infinite system”, which is correct if the neutron flux does not have a spatial gradient for several neutron mean free paths from the domain boundaries. Reflective conditions corresponding to a unitary albedo $\beta = 1$ (i.e. the ratio between the neutron currents entering and leaving the FA) preserve the neutronics environment of the reference assembly and

they are commonly used to represent a FA surrounded by fuel with the same level of burnup [152]. Other albedo boundary conditions are needed when a change of the surrounding media is expected to take place (i.e. around reflectors or strong absorbers such as control rod). Boundary conditions with albedo $\beta > 1$ can be representative of an FA surrounded by fuel with lower burnup or higher enrichment. On the other hand, $\beta < 1$ is found for assemblies adjacent to moderator channels or reflectors. It is actually common to find this strategy along in-core fuel management depletion calculations of commercial reactors in well established codes like CASMO-SIMULATE and PHOENIX-POLCA. The constraint in defining a certain type of boundary condition might have an impact on decay heat predictions and it is tightly connected to the computational approximations introduced with a spatial discretization.

## Radial depletion zones

The spatial dependence of the depletion in equation (62) is typically approximated by adopting a discrete number of volumes (or depletion zones) in which the model is homogenized. Then, for each depletion zone, a set of ODEs is defined using spatially-averaged reaction rates. If the spatial discretization is not carefully defined the local reaction rates may be under/over-estimated in regions where the neutron flux gradients are steep.

The volume subdivision should be physically driven by the system heterogeneities and by the spatial distribution of the neutron flux, e.g. at the model boundaries or in presence of high-absorbing materials. In an FA, depletion zones could follow the pin cell boundaries or include several adjacent pins. In general, a trade-off is sought between the computational time and the representativeness of the mode, knowing that refining the spatial subdivision reduces the discretization bias. In particular, in Monte Carlo burnup codes the smaller the depletion zones are the longer it takes a transport calculation to achieve a reliable estimate of the reaction rates.

In a FA model the flux distribution is affected by the type of applied boundary conditions. Reflective boundaries along all directions (i.e. $\beta = 1$) define an infinite medium with the nonphysical suppression of neutron leakage. Under these circumstances, it is foreseeable that a limited number of depletion zones could correctly reproduce the nuclide buildup in the system, as the neutron flux distribution is approximately flat. Conversely, white boundary conditions in all directions (i.e. $\beta = 0$) impose a steep gradient to the flux and more burnable zones should be used.

A volume discretization based on physical metrics such as the neutron flux or burnup spatial gradient can produce optimal results when it comes to accurately predict the nuclide content in SNF. This is particularly true for nuclides such as $^{134}$Cs, $^{241}$Am, and the curium isotopes whose concentrations do not evolve linearly with burnup. The same conclusions hold for decay heat calculations. Adopting a single depletion zone appears to be a reasonable option only in case of reflective boundary conditions, or for medium-term cooling times, where the most of the
decay heat is released by the decay chains of $^{137}$Cs and $^{90}$Sr, two nuclides whose concentration in SNF is almost linearly dependent on burnup in the case of LWR.

## 2D vs. 3D

At the beginning of life, the axial distribution of the neutron flux in a LWR FA approaches a cosine, with the neutron leakage being responsible of the flux depression at the vertical extremities of the assembly. Spectral effects generally occur in the same regions and depend on the core configuration above and below the fuel, e.g. presence of cladding, structural components, fertile fuel, axial blanket, coolant. An uneven axial flux distribution depletes the fuel at a higher rate near the axial midsection than at the ends. During operation, the burnup buildup decelerates at the center because of the accumulation of fission products and the depletion of the fissile material [157–160]. The burnup axial profile in a fuel assembly at discharge is relatively flat in the axial center with burnup peaks 1.1-fold the assembly average burnup, and considerably low at the extremities, down to about half the assembly average burnup [161]. In addition, the higher moderator density in the FA bottom half typically generates an axial offset with a higher burnup towards the bottom of the fuel assembly compared to the top. Some assemblies might display a different offset due to the proximity of control rods during their operational history, although this effect should in principle disappear at the time of discharge [158].

Bounding burnup axial profiles have been established for the criticality calculations of commercial LWR assemblies to maximize the effects that the low burnup regions near the ends of spent fuel have on the system reactivity, i.e. so-called “end effect”. However, such profiles do not ensure a conservative generation of the spent fuel nuclide concentrations with respect to decay heat applications. For example, using non-irradiated fuel is a conservative assumption for reactivity calculations, but not for decay heat predictions [162].

An accurate determination of the burnup accumulated axially in the fuel removes biases in the decay heat estimates that are otherwise introduced with simplified modeling assumptions. A typical approach to predict SNF integral quantities such as the decay heat rate of a FA implies using two-dimensional (2D) lattice models, with corrections to account for the axial flux variation [163]. On the other hand, despite burnup codes coupled to 3D transport solvers can accurately model the flux profile, often simplifications are introduced, either to reduce the complexity of the system or for lack of information on the operating conditions [164]. Parameters such as the material composition, temperature, and density in the FA are often homogenized to an average value, under the assumption that they have a linear dependence with the observable of interest.

The parametric study reported in reference [165] for a reference PWR FA assembly model (with an initial enrichment of 3.1%, and an average burnup of 36.7 GWd/tHM) shows that the decay heat rate does not scale linearly with modeling parameters such as the moderator temperature at long cooling times and the cycle power at shorter cooling times (see for instance Ref. [166]).

## Time discretization

Despite the latest enormous development in computer technology, depletion codes coupled to Monte-Carlo transport solvers are still considered computationally expensive, primarily because of their more and more frequent application to complex realistic problems such as a whole fuel assembly or a full-core problem. The most common solution adopted to reduce the computational burden is to introduce approximations in the coupled problem. In Monte Carlo burnup codes the solution for the neutron flux field tends to be more computationally expensive than the integration of the depletion equations. Hence, the time steps for complex problems are generally long in order to minimize the number of steady-state neutron transport calculations. Although the neutron flux field generally varies in time at a lower rate compared to many nuclide number densities, a coarse time-discretization might lead to an under- or over-estimation of the reaction rates at which individual nuclide concentrations are built up or depleted. This can result in significant discretization errors. Time-dependent spectral variations within a time step can be accounted for by using predictor-corrector algorithms, with a considerable improvement in the accuracy of the burnup solution [167,168]. It was demonstrated that the use of implicit predictor-corrector algorithms such as the modified Euler method can also prevent stability issues, e.g. spatial oscillations of the neutron flux driven by the neutron capture in $^{135}$Xe or, for longer time steps, by the fuel burnup [148,169]. Adopting advanced burnup algorithms with reaction rate interpolation/extrapolation capabilities can overcome the bias introduced by longer time-steps. Then, the time-dependence of the transition matrix in a depletion time-step can either be inherently handled by the ODE solver [170] or by further dividing the time-step into a number of sub-intervals [171,172].

Shortening the time steps reduces the bias on the burnup solutions, but since the total number of depletion steps is increased, so is the overall computational cost. The optimal time-step length is problem-dependent and, as a rule of thumb, it increases as the irradiated material is depleted. Numerical tests suggested that the computational efficiency of Monte Carlo burnup calculations might be actually improved by shortening the time steps and running transport calculations with a larger statistical uncertainty [173]. Following these results, algorithms are being developed to maximize the efficiency of the Monte Carlo burnup simulation by optimizing the time-step length [174].

## DBRC

The scattering of neutrons in the epithermal energy range can be significantly affected by thermal motion of target nuclides [175]. The probability of change of a neutron's initial energy and direction as a consequence of an interaction with a target nucleus (i.e. sampling the scattering kernel) is traditionally simulated in MC codes in reference [176]: (i) using a free gas model or pre-processed thermal neutron scattering kernel tables $S(\alpha, \beta)$ considering internal changes in the molecule (crystal effect), in
<figure><img src="image_31.png" /><figcaption>**Fig. 28.** Total and electromagnetic decay heat component for pulse irradiation of <sup>235</sup>U (top) and <sup>239</sup>Pu (bottom) with thermal neutrons.</figcaption></figure>

the thermal energy range; (ii) sampling the target nuclide velocity from a Maxwellian distribution in the epithermal energy range; and (iii) assuming the target to be stationary at sufficiently high neutron energies (400 kT).

While the traditional approach in the epithermal energy range is valid for light nuclides, heavy nuclides, in particular <sup>238</sup>U, should include the fact that the cross section is definitely not constant in the vicinity of the resonances [177]. To take this effect into account, resonant dependent probabilities (i.e. S(α, β) tables) for heavy nuclides were introduced [178]. However, this method is limited since the anisotropy of the azimuth angle between the interacting plane and the scattered neutron plane is not considered, and this is needed for the accurate handling of the scattering kernel [179]. Alternatively, Rothenstein [180] proposed a re-write of the elastic scattering probability density function, which is suitable for numerical calculations, and that it is implemented in state-of-the-art neutron transport codes such as MCNP [181], KENO [182] or SERPENT [183] to account for non-constant scattering cross sections.

### 6.4.3 Nuclear data

Nuclear data provide important quantities for the calculation of decay heat following the summation method. Precise information on cross sections (e.g. fission and capture), decay data (such as half-lives, released energies, branching ratios) or fission yields are needed if one

<figure><img src="image_32.png" /><figcaption>**Fig. 29.** Example of weighted cumulative fission yields, before and at the end of irradiation, for two UO₂ and MOX assemblies.</figcaption></figure>

needs to rely on calculated decay heat. Over the past decades, there have been numerous improvements from the nuclear data libraries (improved cross sections, fission yields, decay data), supported by a large range of validation, leading to more accurate summation calculations. Comparisons between various current nuclear data libraries are presented in Figure 28, for the total and decay heat (electromagnetic component) of <sup>235</sup>U and <sup>239</sup>Pu thermal fission pulses.

The experimental values and calculations presented in this figure are obtained from references [58,107,122,124]. Such better agreement between calculations, based on the improvement of nuclear decay data along the years, and experimental values provide confidence in the results from the summation method for SNF.

A detailed study on the need for nuclear data related to decay heat can be found in references [184] and [7]. From equation (62), one can see that the quality of the nuclear data directly affects the estimation of the nuclide concentration, and therefore of the calculated decay heat. As an example, the improvement of the decay heat calculations thanks to the TAGS data (see Sect. 6.4.4) was highlighted in many references (see for instance Refs. [185–187]). The quality of fission yields (cumulative as well as independent) also plays a prime role in the estimation of nuclide concentrations. Figure 29 presents two examples of cumulative fission yields, before and after irradiation of two UO₂ and MOX assemblies in a realistic core environment.

These examples highlight the changes of the so-called “double-hump curve” with burnup, and indicate the variety of nuclides to take into account to obtain an adequate SNF characterization.
<figure><img src="image_33.png" /><figcaption>Fig. 30. Schematic picture of the Pandemonium effect.</figcaption></figure>

### 6.4.4 TAGS data

#### Pandemonium effect

Decay heat calculations are based on summation calculations performed with evaluated decay data, as explained in Section 6.4.1. The average $\beta^-$ and $\gamma$-ray energies released per decay ($\bar{E}_{\beta^-}$, $\bar{E}_{\gamma}$) are obtained from the individual beta decay schemes of the decaying fission products using the $I_{\beta^-}$ feeding to each individual state in the daughter nucleus and their discrete $\gamma$-rays energies. Usually, the $I_{\beta^-}$ intensity distributions are obtained from high-resolution spectroscopy measurements with Ge detectors based on the $\gamma$-ray intensity balance technique which relies on the detection of all individual $\gamma$-rays coming from the levels fed in the $\beta^-$ process and their placement on the beta decay level scheme. This technique suffers from a well known systematic error called *Pandemonium effect* pointed out by Hardy et al. in 1977 [188].

Using a fictional nucleus called Pandemonium with a complex decay scheme, Hardy et al. showed that 14% of the $\gamma$-ray intensity above 1.7 MeV remains undetected under normal experimental conditions using a Ge(Li) detector. Due to the limited peak efficiency of Ge detectors, some $\gamma$-ray transitions from levels at high excitation energy are missed and consequently extra $\beta^-$ feeding is incorrectly assigned to low lying states, shifting the apparent $I_{\beta^-}$ distribution to lower excitation energies, see Figure 30. This leads respectively to an under-estimation of the $\gamma$-ray and over-estimation of the $\beta^-$ energies released in the decay. The largest effect is observed in the decays with large $Q_{\beta^-}$ values because the level density increases with the excitation energy (providing many possible weak $\gamma$ decay paths that can remain undetected) and those where the $\beta^-$ strength concentrates close to the end of the Q-value energy window [189].

#### Total Absorption Gamma-ray Spectroscopy technique

A different experimental technique, the Total Absorption $\gamma$-ray Spectroscopy technique (TAGS), based on the detection of the full $\gamma$-ray cascades which follow the $\beta^-$ decays instead of detecting the individual $\gamma$-rays was developed to

overcome the Pandemonium effect. This technique relies on a large detection efficiency of $\gamma$-rays, almost 100%, which is usually obtained with large scintillator crystals with a geometry as close as possible to $4\pi$. In the 1990s, a total absorption $\gamma$-ray spectrometer coupled to a $^{252}$Cf-based ISOL facility was used in a program of systematic study of the distributions of $\beta^-$ decay intensities of fission products by Greenwood et al. [190] aiming to improve decay data needed for decay heat calculations. Results were obtained for 48 nuclei on Rb, Sr, and Y isotopes and the region of $Z \ge 55$, with half-lives ranging, generally from more than 5 seconds to less than a few hours. Twenty-nine nuclei had their mean decay energies updated in the JEFF-3.1.1 library using these measurements. The nuclei updated, see Tables D.1 and D.2 (in Appendix D), were only those for which a detailed evaluation was not available [191]. No details were given on the analysis and a conservative 10% uncertainty was applied to all mean energies. The evaluation of these nuclei was also taken into account for the release of the ENDF/B-VII.1 library in 2011. Since then, the analysis methodology of the TAGS technique was considerably improved, especially by the Valencia group (see Refs. [192-195]).

In 2005, a WPEC Subgroup 25 was focused on the assessment and improvement of the evaluated decay data sub-libraries in order to obtain more accurate estimations of decay heat [7]. Recommendations were made to focus on a list of 37 particular fission products that merit measurement by TAGS in order to improve decay heat calculations. Consultant's meetings were also organized afterward by the Nuclear Data Section of IAEA to improve/complete the initial list of proposed TAGS measurements for decay heat both for U/Pu and Th/U fuel cycles but also for other possible applications associated with the determination of antineutrino spectra in reactors and the study of beta-delayed neutron emitters [196-198]. The meetings also monitored the progress of the different experimental campaigns and their impact. Four experimental campaigns were performed in 2007, 2009, 2014, and 2022 by the European TAS collaboration (CSIC-University of Valencia) Spain and the University of Rey, UK, joined by the SUBATECH team, Nantes, France in 2009) at the IGISOL facility in Jyväskylä, Finland (see Refs. [199-206]). TAGS experiments were also performed on nuclei from the NEA list at the Holifield Radioactive Ion Beam Facility at the Oak Ridge National Laboratory by the MTAS collaboration (USA, Poland) [207-209]. Tables D.3 and D.4 (Appendix D) list the nuclei that were measured and published by the different teams up to now from the NEA list. The gray cases are the TAGS $\bar{E}_{\beta^-}$ and $\bar{E}_{\gamma}$ already evaluated and included in the JEFF, JENDL, and ENDF libraries as ELP and EEM values. It should be noted that $^{102}$Tc and $^{101}$Nb were not included in JEFF-3.3 from the European TAGS measurements because the mean decay energies were almost identical to the values already in JEFF-3.1.1 for the first one and for the second one there was almost 50% uncertainty on the mean gamma-ray energy [191]. These cases were shown to not suffer from the Pandemonium effect (see Refs. [199,200]).
<figure><img src="image_34.png" /><figcaption>Fig. 31. Left: light particles decay heat component for pulse irradiation of $^{235}$U with thermal neutrons. Right: same for the electromagnetic decay heat component.</figcaption></figure>

## Impact on fission pulse decay calculations

Recently, some decay heat calculations of the Light Particle (LP) and Electromagnetic (EEM) components both for $^{235}$U and $^{239}$Pu thermal fission pulses were performed using the Serpent code [183] to calculate the impact of the TAGS measurements already included in the JEFF, JENDL, and ENDF evaluated libraries [210]. Decay heat calculations based on different releases of the JEFF library are compared in Figures 31 and 32 with the Tobias evaluation, Dickens, and Lowell decay heat measurements [107,122,124]. The decay heat calculations for the thermal fission pulses of $^{235}$U and $^{239}$Pu obtained with different versions of the series of JEFF libraries in their chronological order illustrate clearly the growing impact of the TAGS measurements with a decrease in the LP and an increase in the EEM decay heat components. This is a direct consequence of the adoption of Pandemonium-free decay data.

Decay heat calculations were also performed on thermal fission pulses in order to compare the performance of the latest evaluated libraries: ENDF/B-VIII.0 [211], JEFF-3.3 [58], and JENDL-5 [212]. The decay heat curves for the Light Particle component, both for $^{235}$U and $^{239}$Pu thermal fissions, are in good agreement with the decay heat measurements and their experimental uncertainties. The electromagnetic decay heat component of $^{235}$U thermal fission from 2 to 100 seconds cooling times is underestimated per the JEFF-3.3 calculation whereas both ENDF/B-VIII.0 and JENDL-5 calculations agree. At cooling times above 400 seconds, it is rather difficult to conclude as the three sets of measurements exhibit discrepancies that remain unresolved at this stage. For the electromagnetic decay heat component of $^{239}$Pu, the peak

at cooling times from 30 to 200 seconds is overestimated per the JENDL-5 calculation. Given the importance of $^{239}$Pu thermal fission for decay heat, further investigations were performed by repeating the same calculation using the JEFF-3.3 decay data library combined with three different fission yield libraries: ENDF/B-VIII.0, JEFF-3.3 and JENDL-5. As shown in Figure 33, both calculations performed with the ENDF/B-VIII.0 and JEFF-3.3 fission yield libraries reproduce the experimental data whereas the adoption of the JENDL-5 fission yield library leads similarly to an overestimation of the decay heat measurements in the 30–200 seconds cooling range. This overestimation would appear to be associated with the fission yield library and is now studied in detail per the JENDL team with respect to a release of fission yield data [210].

## 6.5 Methods for depletion solvers

The analytical solution of the Bateman equation is the matrix exponential function:

$$n(t) = \exp(At) \cdot n(t_0) \qquad (65)$$

where $n(t_0)$ is the nuclide density vector at time $t_0$. The part of the transition matrix A representative influencing the natural decay or the neutron-induced activation of a given nuclide usually contains reaction rates for only a few dozen of nuclides. On the other hand, the time evolution of the fuel composition in a nuclear reactor requires already 30–40 fission products and actinides that characterize the major SNF observables. In addition, several of these nuclides are generated via the decay or neutron-
<figure><img src="image_35.png" /><figcaption>Fig. 32. Left: light particles decay heat component for pulse irradiation of $^{239}$Pu with thermal neutrons. Right: same for the electromagnetic decay heat component.</figcaption></figure>

capture processes involving other nuclides, thus extending the list to several hundreds. Eventually, if one wanted to consider all nuclides for which nuclear data are provided in the evaluated radioactive decay data libraries, the number would increase to almost 4000, e.g. the JEFF-3.3 radioactive decay data library provides data for 3852 nuclides.

The exponential of a matrix can be computed in many ways involving different approximations [213], still in a transmutation-rich environment their implementation is not straightforward and requires to take into account the following considerations:

1. half-lives ($T_{1/2} = log(2)/\lambda$) can vary from billions of years to small fractions of seconds, thus generating a large spread in the transition matrix eigenvalues and making the system numerically stiff;

2. the large decay constants associated with highly unstable nuclides (e.g. $^{7}$B with a half-life $T_{1/2}$ of about $10^{-22}$ seconds) make the problem ill-conditioned, with a matrix norm up to $||A|| \simeq 10^{21}$;

3. long time steps also drastically increase the conditioning of the problem.

The main challenge in solving the depletion problem lies not in the mathematical form of the Bateman equations governing the changes in the nuclide concentrations, but rather in the large dimensions and the numerical characteristics of the system. In practice, however, there exists a wide range of state-of-the-art solvers that are capable of handling the depletion problem with a sufficient level of accuracy, and errors resulting from their limitations can be considered negligible compared to other factors. Comparisons of depletion algorithms in terms of performance and numerical accuracy are also reported in references [213,214]. A number of examples are presented below.

### Transmutation trajectory analysis (TTA)

The TTA method [215], also called the linear chain method, implies the analytical solution of one Bateman set of equations per each individual transmutation chain. The time evolution of the nuclide vector is broken down into chains that represent every possible reaction path. This method is used also for non-linear chains after they have been broken down into sets of linear chains. The total nuclide number densities are obtained from the sum of the results of each sub-chain for all the nuclides involved in the problem. Examples of codes that use this method are BISON [216], CINDER [217], MCB [218,219], and Serpent [183].

The depletion of nuclear fuel involves an enormous network of complex sub-chains that often happen to generate loops that cannot be linearized. As a consequence, assumptions are made to terminate unimportant chains based on multiple criteria, such as ignoring cyclic chains or terminating them after the first few loops. The choice of the termination criterion can play a role in the trade-off between computational speed and accuracy.

### Matrix exponential methods

The solution of the Bateman equation can be formally expanded as the Taylor series of the matrix exponential:

$$ \exp(At) = I + A \cdot t + \frac{(A \cdot t)^2}{2!} + \dots \quad (66) $$
<figure><img src="image_36.png" /><figcaption>Fig. 33. Impact of fission yield libraries on electromagnetic decay heat component for pulse irradiation of $^{239}$Pu with thermal neutrons.</figcaption></figure>

The matrix exponential method includes an approximation of the matrix exponential by a truncated Taylor series. This is implemented in the depletion code ORIGEN, available in the SCALE suite [220, 221], and adopted in many coupled neutron transport and fuel depletion codes (see for instance Monteburns [222], MOCUP [223], MCODE, ACAB [224], EVOLVOCODE [225] or VESTA [226]).

The ORIGEN method involves solving the matrix exponential in three phases, in which the contribution of short- and long-lived nuclides is separated. The short-lived nuclides initially present at time $t_0$ in the nuclide vector are treated independently by constructing for each nuclide all the populating chains consisting of only short-lived nuclides and solving them in a way analogous to the TTA method. In the second phase, the contributions of short-lived nuclides to long-lived nuclides are added to the initial densities of the latter. The transition matrix A is then re-constructed by assuming that the short-lived nuclides decay instantly, thus practically removing them from the system. By having removed the instabilities associated with short-lived nuclides, the matrix exponential can now be solved by using a power series. Eventually, the contributions of long-lived nuclides to short-lived nuclides are obtained iteratively using the results of the second phase as an initial guess and assuming a secular equilibrium

$$ \frac{dn(t)}{dt} = A \cdot n(t) = 0, \qquad (67) $$

valid for short-lived nuclides.

## Ordinary Differential Equation Methods

The Bateman equation adopted to describe the time evolution of nuclear fuel composition in a reactor is a system of ordinary differential equations (ODE), therefore it is natural to consider its resolution via numerical integration methods [227]. Ready-to-use packages of ODE solvers [228] are nowadays distributed for the most common computer programming languages.

Among the codes that rely on such methods, there is ALEPH-2 [229], which uses the RADAU IIA implicit Runge-Kutta method of order 5. This method, as well as

others involving numerical integration, is easy to use and requires only minimal programming for its implementation. RADAU IIA showed excellent stability for systems containing 3000 and more known nuclides, with the primary disadvantage being the relatively high cost in computer time compared to other solvers like ORIGEN. Still, these methods have been often developed for general nonlinear differential equations, and could potentially achieve much shorter running times if optimized for specific ODE problems.

Other numerical integration schemes (Runge-Kutta, BDF) could be considered for depletion calculations, but are limited by the nuclides with short half-lives that make the problem stiff. Such methods could be used if the short-lived nuclides were handled separately as done by ORIGEN. This is the case of the newest depletion module PHOENIX included in VESTA [226], which solves the Bateman equation either with a fourth-order Runge-Kutta method, a Taylor series development of the matrix exponential or the method of Arnoldi [230], after the separate treatment of the short-lived nuclides.

## Chebyshev rational approximation method (CRAM)

CRAM is a method to solve the matrix exponential equation that is implemented e.g. in Serpent [183]. It involves approximating the exponential function $e^z$ as a ratio between two polynomials $P_k(z)$ and $Q_k(z)$ of order $k$:

$$ e^z \approx \frac{P_k(z)}{Q_k(z)}. \qquad (68) $$

Eventually, the solution of the Bateman equation can be retrieved from solving $k/2$ linear systems using LU decomposition. The main advantage of this method is that for a chosen polynomial order $k$ and range of values $z$, the coefficients of the polynomials $P_k$ and $Q_k$ need to be determined only once. This method has been proven particularly effective for sparse matrices and, with order $k = 16$ polynomials, it can provide very accurate solutions to burnup equations without excluding any nuclides. More details can be found in reference [231].

## 6.6 Sources of uncertainties

As for many other calculated quantities based on reactor simulations involving neutron transport and thermal-hydraulics modules, the assembly decay heat can be affected by various parameters. Additionally, such parameters may not affect the calculated decay heat in the same way as a function of cooling time. For instance, quantities of relevance during (or just after) reactor transients may not be the same as for long-term cooling. As it is challenging to be exhaustive within this study, we will provide a number of examples, mostly for cooling time longer than months after shutdown.

Quantities of relevance, being significant sources of uncertainties for the assembly decay heat are generally
<table><caption>**Table 17.** Uncertainties for operational parameters, manufacturing tolerances, and burnup-induced technological changes.</caption><thead><tr><th rowspan="2">Parameters<br>Probability</th><th rowspan="2">Uncertainty [232]<br>–</th><th colspan="6">Design and operational parameters</th></tr><tr><th>1σ [233]<br>Normal</th><th>1σ [33,234]<br>Uniform</th><th>1σ [235]<br>Uniform</th><th>1σ [236]<br>See superscript</th><th>1σ [237]<br>See superscript</th><th>1σ [238]<br>Normal</th><th>1σ [239]<br>–</th></tr></thead><tbody><tr><td>Cladding thickness</td><td>± 40–50 µm</td><td>16.7 µm</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Cladding diameter (PWR)</td><td>± 200 mm</td><td>67 µm</td><td>0.5%</td><td>0.8%</td><td></td><td>0.1<sup>U</sup>%</td><td>0.1%</td><td>0.1%</td></tr><tr><td>Cladding diameter (BWR)</td><td>± 300 mm</td><td>100 µm</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Pin x and y position shifts</td><td>–</td><td>–</td><td>0.5%</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Fuel pellet density</td><td>≤ 2%</td><td>0.67%</td><td>0.5%</td><td></td><td></td><td></td><td>0.5%</td><td>1.65%</td></tr><tr><td>Fuel pellet diameter</td><td>± 20 µm</td><td>correlated with fuel pellet dens.</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Enr. (<sup>234</sup>U wt.%)</td><td></td><td></td><td></td><td>5%</td><td></td><td></td><td></td><td></td></tr><tr><td>Enr. (<sup>235</sup>U wt.%)</td><td>± 0.05%</td><td>0.0167%</td><td>0.2%</td><td>1%</td><td></td><td></td><td>0.02%</td><td>0.05%</td></tr><tr><td>Enr. (<sup>238</sup>U wt.%)</td><td></td><td></td><td></td><td>0.2%</td><td></td><td></td><td></td><td></td></tr><tr><td>Enr. (<sup>239</sup>Pu wt.%) (MOX)</td><td></td><td></td><td>0.2%</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>SNF power</td><td>–</td><td>1.67%</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Burnup depletion step</td><td></td><td></td><td>0.25%</td><td>1%</td><td></td><td></td><td></td><td></td></tr><tr><td>Water temp. (PWR only)</td><td>± 2 °C</td><td>2 °C</td><td>2%</td><td>2%</td><td>1.5<sup>U</sup>%</td><td></td><td></td><td></td></tr><tr><td>Water density (PWR only)</td><td>± 0.005 g/cm<sup>3</sup></td><td>0.005 g/cm<sup>3</sup></td><td></td><td>2%</td><td></td><td></td><td></td><td></td></tr><tr><td>Void fraction (BWR only)</td><td>± 6%</td><td>6%</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Fuel temperature</td><td>± 50 °C</td><td>50 °C</td><td>2%</td><td>2%</td><td></td><td></td><td></td><td></td></tr><tr><td>Boron content (PWR only)</td><td>± 10 ppm</td><td>10 ppm</td><td>2%</td><td>2%</td><td></td><td></td><td></td><td></td></tr><tr><td>Reactor pressure</td><td></td><td></td><td>1%</td><td></td><td>1<sup>N</sup>%</td><td></td><td></td><td></td></tr><tr><td>Sample power</td><td></td><td></td><td></td><td></td><td>1.5<sup>N</sup>%</td><td></td><td>2<sup>N</sup>%</td><td></td></tr><tr><td>Pin x and y position shifts</td><td>–</td><td>–</td><td>0.5%</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table><table><thead><tr><th rowspan="2">Probability</th><th colspan="6">Burnup-induced technological changes</th></tr><tr><th>1σ<br>Uniform</th><th>1σ [240]</th><th>1σ [241]</th><th>1σ [242]</th><th>1σ [243]</th><th>1σ [244]</th></tr></thead><tbody><tr><td>Fuel pin displacement</td><td>0.2 mm/burnup steps</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Moderator pin displacement</td><td>0.05 mm/burnup steps</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Pin radius</td><td>150 µm<br>@ 53 MWd/kg</td><td>15–50 µm<br>@ 28 MWd/kg</td><td>70 µm<br>@ 50 MWd/kg</td><td>80 µm<br>@ 30 MWd/kg</td><td>160 µm<br>@ 61 MWd/kg</td><td>150 µm<br>@ 60 MWd/kg</td></tr></tbody></table>

<table><caption>**Table 18.** Example of VVER-1000 nuclear fuel uncertainty ranges.</caption><thead><tr><th rowspan="2">No.</th><th rowspan="2">Parameter</th><th rowspan="2">Symbol</th><th rowspan="2">Units</th><th colspan="2">Range</th><th rowspan="2">Type</th></tr><tr><th>Negative</th><th>Positive</th></tr></thead><tbody><tr><td>1</td><td>Inner diameter ( only TVS-M and TVSA)</td><td>Δd<sub>fin</sub></td><td>mm</td><td>0</td><td>+0.2</td><td>Absolute</td></tr><tr><td>2</td><td>Outer diameter</td><td>Δd<sub>out</sub></td><td>mm</td><td>-0.03</td><td>0</td><td>Absolute</td></tr><tr><td>3</td><td>Density at 20°C</td><td>Δρ<sub>UO<sub>2</sub></sub></td><td>kg/m<sup>3</sup></td><td>-150</td><td>+150</td><td>Absolute</td></tr><tr><td>4</td><td>Length</td><td>L</td><td>m</td><td>n.a.</td><td>n.a.</td><td>n.a.</td></tr><tr><td>5</td><td>Effective full power days</td><td>δT<sub>eff</sub></td><td>–</td><td>-0.03</td><td>+0.03</td><td>Relative</td></tr><tr><td>6</td><td>Power of the most loaded FA</td><td>δN<sub>FA</sub></td><td>–</td><td>-0.05</td><td>+0.05</td><td>Relative</td></tr></tbody></table>

classified into different groups: nuclear data, operating conditions, manufacturing tolerances, and eventually burnup-induced technological changes. A number of values for uncertainties from the literature are presented in Table 17. Details are provided below.

### 6.6.1 Method of uncertainty propagation

In the case of depletion calculations followed by a cooling period, the method of uncertainty propagation mostly used for such non-linear and time-dependent cases is based on simple Monte Carlo: repeating a number of times the same calculation, but each time with a different set of input parameters. This approach allows to access an estimated

probability density function for the quantity of interest, from which the uncertainty can eventually be represented by a standard deviation. Such a method is applicable for decay heat, as well as for any other quantities similarly obtained, such as nuclide concentrations or assembly and reactor quantities. Many publications using this method can be found in the literature. One limitation of this approach is due to the required calculation time, which is virtually multiplied by the number of sampled cases. It is therefore not uncommon to estimate uncertainties based on a few hundred of sampled runs. A number of acceleration solutions exist, depending on the type of transport method used: Monte Carlo or deterministic.
Another method of uncertainty propagation is based on sensitivity calculations, under the linearity assumption, combined with the so-called “sandwich formula”. This method requires to calculate the sensitivity of the quantity of interest to a specific input parameter, and then to perform a matrix multiplication between sensitivity vectors and covariance matrices. If this method is in principle faster than the Monte Carlo uncertainty propagation, it is based on the linearity assumption and can become cumbersome in the case of a large number of input quantities (such as nuclear data).

### 6.6.2 Nuclear data

The impact of nuclear data on decay heat was presented in a number of publications, see for instance references [12,33,34,52,165,233–235,245,246]. The usual approach to assess the uncertainties due to nuclear data is to consider the covariance matrices provided in the evaluated nuclear data libraries, generate a set of sampled cross-sections or decay data, and repeat transport calculations with these different nuclear data. As long as covariance matrices are provided, the user is well-guided. In the case of missing information, such as for fission yields, users can apply different methods, effectively leading to different calculated uncertainties.

Nuclear data of importance for decay heat are generally cross sections (fission, capture, for both fission products and actinides) and fission yields (mostly from $^{235}$U and $^{239}$Pu). Decay data (released energy, decay mode, half-lives) induce uncertainties overall smaller than cross sections and fission yields, but their treatment is less rigorous due to the lack of covariance information in evaluated libraries.

Finally, it can be mentioned that nuclear data is one of the major contributors to the decay heat uncertainties, as generally indicated in the mentioned references.

### 6.6.3 Manufacturing tolerances

Manufacturing tolerances include uncertainties on the design dimensions of fuel assemblies (such as pin radius, pitch, and cladding thickness), on the fuel density, on the initial heavy metal contents (uranium, plutonium and eventually americium isotopes) and possibly neutron absorber contents (such as gadolinium isotopes). Their impact is not negligible, as shown in the previously mentioned references.

#### Initial geometry and density uncertainty

The proper evaluation of decay heat rates requires the selection of appropriate input uncertainties of stored FA in the Spent Fuel Pool and density uncertainties. For example, if the fuel is provided in the form of cylindrical pellets, the equation could be applied accounting geometrical and density uncertainties: $\frac{\pi}{4} \rho_{UO_2} (d_{\text{fout}}^2 - d_{\text{fin}}^2)$, where $d_{\text{fout}}$, $d_{\text{fin}}$, $\rho_{UO_2}$ and $l$ are the outer column diameter, inner column diameter (if exists), density, and fuel length, respectively. The fuel column length could be the total reactor core length or the individual fuel pellet. As

the column length uncertainty range, due to vendor procedures compensating the other geometrical uncertainties, then it is not accounted for in this section. An example of VVER-1000 fuel pellet uncertainty ranges is provided in Table 18 [247,248].

The analytical estimation of the geometry uncertainty generated the following equation (69):

$$ \pm \Delta m_{1,2} = \pm \frac{\pi}{4} l \rho_{UO_2} \left[ 2 \Delta d_{\text{fout}} \Delta d_{\text{fout}} + 2 d_{\text{fin}} \Delta d_{\text{fin}} + \frac{\Delta \rho_{UO_2}}{\rho_{UO_2}} (d_{\text{fout}}^2 + d_{\text{fin}}^2) \right]. \quad (69) $$

The description of the symbols, used in equation (69) is given in Table 18. In the case if no internal opening exists in fuel pellets by the design, then $\Delta d_{\text{fin}} = 0$. The uncertainty range is technologically bounded within the following limits: $m - \Delta m_1 \le m \le m + \Delta m_2$, where the particular absolute mass uncertainty could be evaluated based on data from Table 18. In this particular case, the lower and upper portions of the FA mass uncertainty range have different sizes, which indicates a non-symmetrical statistical distribution, due to the vendor's production of technological geometrical uncertainty ranges (provided in Tab. 18). However, for other types of fuel made by the same or other vendors, this analysis could be different.

#### Fuel mass uncertainty

The fuel mass uncertainty is an added value to the initial geometrical and density uncertainties, based on the operational procedure for the fuel production. The limit could be set so that the total mass should not exceed a threshold, which usually is the nominal value, considered by the fuel column design. In this case, the geometrical uncertainty could be upgraded to $m - \Delta m_1$, which transforms the original two-sided statistical uncertainty to one side. Therefore, the geometrical uncertainty is reevaluated, according to the fuel vendor procedures, in order to avoid additional non-realistic conservatism.

#### Enrichment uncertainty

Usually, the enrichment uncertainty is provided by the fuel vendors where the manufacturer's fuel enrichment uncertainty was considered to have a limit as presented in Table 17. However, during recent years some of the facilities, operating nuclear fuel adopted technology for enrichment evaluation by the use of a compact detector (as an example CdZnTe), which is inserted into the instrumentation channel of the examined assembly and registers the low-energy gamma-ray and characteristic radiation from the isotopes of uranium and their decay products (see Refs. [52,249]). As an example, the fuel enrichment measurement procedure at the Kozloduy NPP involves the measurement of 185.7 keV peak of $^{235}$U. In this particular case, the measurement time is 600 seconds, and the number of counts varies within the range of 5000 to 6000.
<figure><img src="image_37.png" /><figcaption>Fig. 34. Total decay heat generation rate in Kozloduy NPP unit 6, Spent Fuel Pool: (a) 2014-2015 uncertainty evaluation; (b) 2018-2019 uncertainty evaluation; (c) 2014-2015 bounding data; and (d) 2018-2019 bounding data.</figcaption></figure>

### 6.6.4 Operating conditions

Operating conditions of importance for decay heat uncertainties can be numerous, depending on the cooling time of interest: the fuel and moderator temperatures, reactor pressure, boron concentration, assembly power (or burnup), moderator density, and irradiation history.

#### Power history

The uncertainty in the power of the FAs can be obtained from the analytical core design software. Such software passes extended verification and validation procedures and depending on the results, a qualified uncertainty range is produced (generally in terms of “rms”, or root-mean-square). For example for the VVER-1000 core design software, the burnup relative uncertainty $\pm\delta B$ is evaluated as follows reference [250]:

$$ \pm\delta B = \pm\frac{\Delta B}{B} = \pm\left(\delta N_{FA} + \delta T_{eff} + \frac{\Delta m}{m}\right). \quad (70) $$

#### Moderator density (coolant pressure, temperature)

The decay heat rates can be influenced by the change of characteristics for the moderator, as it was indicated in a number of publications. As indicated in reference [246] for BWR calculations, the moderator density during irradiation has a non-negligible impact on local decay heat, being higher than 1% at 100 years of cooling time. For PWR cases, the impact of the moderator temperature reaches similar values [33].

### 6.6.5 Example: GRS best estimate method uncertainty analysis

This example of evaluation of the output uncertainties of FA decay heat generation rates is based on a method previously proposed by the Gesellschaft für

Anlagen- und Reaktorsicherheit (GRS) [251,252] and is currently applied for some safety analysis. This method allows realistic prediction of output uncertainties, calculated after the evaluation of selected input parameters uncertainty propagation. The selected input parameters are expected to be randomly varied within the defined ranges of their minimal and maximal values. The necessary number of calculations is determined as suggested previously by Wilks [253,254]: $1-u^n-n(1-u)u^{(n-1)} \le v$ for two-sided statistical tolerance intervals. Correspondingly for one side of statistical distribution, the proposed formula is: $1-u^{(n-1)} \le v$ [253,254].

The parameter $u$ is compared with the probabilistic acceptance criterion (percentile), while the parameter $v$ is the confidence level that accounts for the influence of the sampling error from a randomly chosen sample of limited size. The formula predicts $n = 93$ separated calculations with 95% probability for the 95% confidence level for two-sided statistical distribution and $n = 53$ separated calculations for one-sided statistical distribution. In general, there is no expectation that some values in the input parameters uncertainty ranges are more likely than others. Therefore, usually, uniform probability density functions (PDFs) are assumed, which is a widely accepted choice for this type of calculation. As a result by application of this method, a set of curves is expected to be obtained for the output decay heat generation rate, where the bounding maximal and minimal values form the propagated output uncertainty range.

Figure 34 presents an evaluation of the decay heat generation rate in Spent Fuel Storage onsite in Kozloduy NPP (Unit 6, VVER-1000). In the analysis were considered two storing periods: from 2014 to 2015 and from 2018 to 2019. The uncertainty calculations were performed based on the initial geometrical and density uncertainties, combined with power operation uncertainties evaluated using RG US NRC 3.54.

The uncertainty analysis was performed for two FA batches stored in the spent fuel pool. All calculations were based on the procedure provided by US NRC RG 3.54. The described input parameters in section three were varied randomly within their defined uncertainty ranges (see Tab. 18). The set of 98 calculations included $n = 93$ random combinations of the following parameters: fuel pellet inner diameter and outer diameters, $UO_2$ fuel density, full power effective days, burnup, and nominal operation power uncertainty (Tab. 18 [250]). The results of the calculation indicated that the decay heat rate uncertainty range is sensitive to the variation in the irradiation conditions, initial mass, and type of fuel assemblies, stored in the spent fuel pool. The results are plotted in Figure 34, where the independent variable of the plots is given in logarithmic scale. The conservative approach states a fixed upper bound of values, where a 10% conservatism was assumed, added to the estimated nominal values for the spent fuel pool [255].

These results were verified additionally by SCALE calculation. The evaluation of the uncertainties in the input values and their propagation to the output identifies the ranges of the variation in the overall spent fuel pool decay heat rates.
### 6.6.6 Alternative approach: safety factor definition

The uncertainties taken into account in the evaluation of decay heat generation rates can be produced by (1) uncertainty in nuclear data, (2) numerical computational models, and (3) inaccuracies from data numerical interpolation and handling of the irradiation history, see for instance reference [28]. These uncertainties are considered by the scaling coefficient and (safety factor) $F_s > 1$ in this reference. A reasonable amount of conservatism is then applied to the computed decay power values. The approved values of the safety factor in reference [28], derived from calorimetric measurements, are reported in Gauld and Murphy [52]. The uncertainty upper bound, which for instance is estimated numerically by reference [28], lay approximately +6% above the experimental data, without the use of the conservative assumption as a safety factor [52]. Furthermore, Gauld and Murphy [77] recommended a safety factor of 1.02 (maximum +8% uncertainty) for a PWR FA stored less than 25 years. Therefore, the safety factor ensures reasonable conservatism in the evaluation of the decay heat generation rate for a PWR. The decay heat rates calculated by the same authors using ORIGEN, compared with the measured values, produced an uncertainty of 4% [52].

## 7 Computer codes

In the following, a number of codes used in the estimation of assembly decay heat for LWRs are presented. For details, a number of references are provided in each case.

### 7.1 CEA computer codes: DARWIN/PEPIN2 and MENDEL

Two major codes devoted to nuclide inventory and decay heat calculation are available at CEA. They solve simultaneous ordinary differential equations describing the transmutation, the growth and the decay of the nuclide densities, and perform an accurate depletion calculation with a fine description of the irradiation history and the isotopic chain. From nuclide concentration results, a large range of physical quantities can be calculated like nuclide mass, radioactivity, decay heat, decays α-particle, β-electron, γ-ray emission, neutron source from spontaneous fission, delayed neutron, and (α,n) reaction. These physical quantities can be computed at any cooling time.

The first of the two codes is DARWIN/PEPIN2 [256] which is a modular code system, which means it is composed of different modules where each of them has specific features. It has been developed and maintained since the mid-90s. In the field of reactor physics studies, it can be linked to CEA neutron transport codes APOLLO2 which provide the necessary neutron data for the nuclide inventory calculation. It implements a numerical fourth-order Runge-Kutta method to solve the Bateman equation for irradiation period, and an analytical solution algorithm for the cooling period. In the irradiation periods, short half-life nuclides are computed apart (consid-

ering a saturated state) in order to obtain a faster solving of the numerical system. DARWIN/PEPIN2 code is a part of DARWIN package which is mainly used with APOLLO2 code in a validated calculation scheme for fuel cycle nuclide inventory or decay heat calculations, with consistent nuclear data libraries based on JEFF-3.1.1 evaluation. DARWIN/PEPIN2 is currently the reference code used by the CEA and its industrial partners EDF and FRAMATOME for fuel cycle studies through the DARWIN package.

In recent years, in the consistency of development of both the new generation 2D/3D transport code APOLLO3® [257] and the new nuclear data processing code GALILEE [258] at the CEA, a new depletion code-named MENDEL has been developed [259]. All of these codes are developed in a modern object-oriented architecture with C++ language, allowing a very strong interoperability between them. For example, MENDEL shares its depletion solver library with APOLLO3® and TRIPOLI-4® [260]. A specific neutron characteristics library named MPO (as Multiple Parameters Output) is provided by transport calculation with APOLLO3® and used by MENDEL for the nuclide inventory calculation. MENDEL implements CRAM (Chebyshev Rational Approximation Method [261]) in addition to the methods also available in DARWIN/PEPIN2 code (RK4 and analytical).

For both codes, once all nuclide concentrations $N_i(t)$ are known, other time-dependent physical quantities of interest (Q.O.I) can be computed: mass, activity, decay heat, radiotoxicity, α-particle emission, β-electron emission, γ-ray emission, neutron source from spontaneous fission, neutron source from (α,n) reaction, delayed neutron. At each cooling time $t$, the total decay heat is derived from the following expression:

$$DH_{\text{tot}}(t) = \sum_{i=1}^{n} \lambda_i \times N_i(t) \times (\overline{E}_i^{\alpha} + \overline{E}_i^{\beta} + \overline{E}_i^{\gamma}), \quad (71)$$

where:

* $\lambda_i$: is the decay constant of nuclide i

* $N_i(t)$: is the concentration of the nuclide i at cooling time $t$

* $\overline{E}_i^{\alpha}, \overline{E}_i^{\beta}, \overline{E}_i^{\gamma}$: are respectively the decay energy emitted by α, β and γ disintegration.

Q.O.I can be calculated in the three following operating modes:

* Nuclear reactor fuel depletion followed by cooling period.

* Neutron activation of material structure followed by cooling period.

* Cooling of initial radioactive material.

Uncertainty quantification propagating nuclear data uncertainties to Q.O.I (in particular for decay heat) is available in both codes:

* Deterministic direct first-order perturbation (one at a time approach) in both DARWIN/PEPIN2 (INCERD module) and MENDEL.
* Stochastic sampling and propagation in MENDEL (sampling can be done or by MENDEL either by an external sampler, scripts are given in order to use URANIE freeware [262]).

Uncertainty Quantification is the last step of the VVUQ process (Verification and Validation with Uncertainty Quantification) at the CEA. Its aim to give users a certain level of confidence in the numerical results, adding uncertainty bars and validating the fact that the observed discrepancies between numerical and experimental values are due to, respectively, nuclear data (and if possible which one), computer code accuracy and options, physical modeling (such as approximation on a geometry or a history of power), etc.

For the CEA computer codes, two approaches have been implemented. In DARWIN/PEPIN2, a first-order direct perturbation method (one at a time) approach has been introduced with the INCERD module. This approach consists of modifying one after the other all the parameters X (i.e. nuclear data) and to compute the Q.O.I. (Y) sensitivity to each parameter. In the second step, the variance of Y is computed through the classical sandwich formula:

$$ \operatorname{cov}(Y) = S_{Y/X}^{t} \operatorname{cov}(X) S_{Y/X}. \qquad (72) $$

Input uncertain parameters can be radioactive decay constants, radioactive decay energies, radioactive decay branching ratios, independent fission yields, multigroup microscopic cross sections, neutron-induced reaction branching ratios, α-particle or γ-ray intensities. Correlations can be taken into consideration, in particular for fission yields, branching ratios, or cross sections. Q.O.I. can be nuclide densities or activities, total activity, total decay heat, α-particle or γ-ray decay spectra.

In MENDEL, the first-order direct perturbation has been implemented for the same input parameters and Q.O.I. Furthermore, a stochastic approach based on random correlated sampling is also available. Sampling can be done either directly in MENDEL (from MENDEL v3.0, non-stratified sampling only for the moment), or using an external sampler like URANIE [263] (scripts provided with MENDEL are using LHS approach). Comparisons [94,264] between the two approaches and different codes have shown a global equivalence of both methods, and discrepancies can be explained.

The nuclear data provided for the two main CEA's evolution code systems, DARWIN/PEPIN2 and MENDEL, come from many international libraries as JEF-2, JEFF-3.1.1, JEFF-3.2, JENDL-4.0, JENDL-2015, ENDF/B-VII.1, ENDF/B-VIII.0, EAF-2003 and EAF-2010. The processing of these data is carried out by CEA's main verification and processing system GALILÉE (release 0.3). The libraries for DARWIN/PEPIN2 and MENDEL, issued from the data processing by GALILÉE, consist of two main nuclear data types: the multigroup cross section (GENDF format) and the radioactive decay data.

The uncertainty data of many physical quantities are also included in these libraries. These data are used to calculate the uncertainty of decay heat, nuclide concentration, and radiation spectra (α-particle, γ-ray, and neutron)

by DARWIN/PEPIN2 and MENDEL. The correlation between some physical quantities like independent fission yields, radioactive decay branching ratios, and neutron cross-sections, is taken into account in the uncertainty calculation. In the case of the neutron cross section, the covariance matrices come from COMAC (the nuclear database COvariance MAtrices from Cadarache).

## 7.2 SCALE

The SCALE modeling and simulation suite [182] for nuclear safety analysis and design, which is developed and maintained by Oak Ridge National Laboratory (ORNL), includes several depletion capabilities that enable calculation of decay heat in spent nuclear fuel for LWR or other reactor configurations. The ORIGEN point depletion and decay code, used worldwide for a variety of applications that require the calculation of nuclide inventories, decay heat, or neutron and gamma radiation emission rate and energy spectra, is the foundation of all depletion capabilities in SCALE. Since its first version [265] was deployed in 1973, ORIGEN [266] has evolved over the years, with continued modernization of the code, underlying methods, and associated nuclear data. ORIGEN includes two solvers for the depletion and decay equations: a hybrid matrix exponential/linear chains method (MATREX) and a Chebyshev Rational Approximation Method (CRAM). Under the TRITON [267] and Polaris [268] depletion sequences, ORIGEN serves as a depletion and decay solver, being coupled with a 1D, 2D, or 3D neutron transport solver to enable depletion and decay simulations for a variety of configurations. As a standalone code, ORIGEN has additional capabilities: (a) simulate continuous nuclide feed and chemical removal, which can be used to model reprocessing or liquid fuel systems and (b) generate alpha-particle, beta-electron, neutron, and gamma-ray emission spectra. ORIGEN is also the engine under the ORIGAMI (ORIGEN Assembly Isotopics) tool for performing rapid depletion and decay calculations to determine source terms for an LWR UO₂ fuel assembly, in a 0D, 1D, 2D, or 3D modeling approximation of the fuel assembly.

ORIGEN tracks 174 actinides, 1149 fission products, and 974 activation products. It explicitly models all available nuclides and transitions in the current nuclear data for decay and neutron-induced transmutation in ENDF/B-VII.1. Cross-section data not available in ENDF/B-VII.1 are obtained from the JEFF-3.0/A special purpose European activation library. Nuclear decay data, fission product yields, and gamma-ray emission data are based on ENDF/B-VII.1 evaluations. The SCALE cross-section libraries that are available for the neutron transport solvers under Polaris or TRITON, in multigroup (MG) or continuous energy (CE), are based on ENDF/B-VII.1 or -VIII.0 evaluated data. Pre-generated, burnup-dependent ORIGEN reactor libraries based on ENDF/B-VII.1 are available in SCALE for a variety of reactor types and assembly designs; for LWR assembly configurations, these libraries cover initial enrichments up to 8.5% and burnups up to 82.5 GWd/t.
<figure><img src="image_38.png" /><figcaption>**Fig. 35.** Decay heat simulation methodology with SCALE [269]. Left: generation of ORIGEN reactor library. Right: ORIGEN depletion and decay simulation.</figcaption></figure>

A rapid approach for determining an LWR assembly decay heat with SCALE is illustrated in Figure 35.

It includes two main computational steps: (1) generation of burnup-dependent ORIGEN reactor libraries for an LWR assembly of interest, using the TRITON or Polaris depletion sequence and (2) depletion and decay simulation with ORIGEN under the ORIGAMI tool. An alternate approach may consist of standalone ORIGEN decay-only simulations to a cooling time of interest, using as input the nuclide inventory binary file generated as a by-product by either TRITON or Polaris. An advantage of using ORIGAMI from an end-user perspective is the easiness of use and the ability to determine in seconds or minutes the assembly-average decay heat or the decay heat variation across the axial dimension of the fuel assembly. However, depending on the heterogeneity of the modeled configuration and the desired level of detail for the decay heat calculation, or the need for additional information, an end-user can take advantage of the flexibility of the depletion sequences or other tools in SCALE. For example, decay heat uncertainty analysis can be performed by using a model developed for the TRITON or Polaris depletion sequence within the Sampler uncertainty quantification in SCALE. That allows assessing the impact on the calculated decay heat of uncertainty in nuclear data (cross sections, fission yields, decay data) or input model parameters (e.g., initial enrichment, temperatures, geometry data).

## 7.3 SNF

### Overview

Studsvik's approach to spent nuclear fuel analysis combines nuclide concentrations, neutron fluxes, and cross sections, calculated by the neutron transport and depletion code CASMO5 [270], with irradiation history data from the reactor core simulator SIMULATE5 [271] and tabulated radionuclide decay data. These data sources are used and processed by the SNF code [272,273] to compute the basic spent nuclear fuel characteristics: nuclide concentrations; radioactivity; decay heat power; photon sources and

spectra; neutron source and spectra from spontaneous fission; neutron source and spectra from (α,n) reactions. This code system is fully integrated with Studsvik's Core Monitoring System (CMS5) and implicitly employs the three-dimensional (3D) discretization of the assemblies in the reactor core. Designed and implemented as a production tool, it provides comprehensive results per axial node of a fuel assembly or a fuel pin. In addition to the basic spent fuel characteristics, the SNF code is supplied with various optional modules intended to support and facilitate data preparation for analyses related to handling, transportation, and storage of spent fuel, including accident analyses, burnup credit in criticality safety evaluations, and dry cask loading scenarios. The decay heat power, computed by the summation method, is complemented by three standards for decay heat power in light water reactors: ANSI/ANS-5.1-2014 [136], DIN-25463-1/2-2014 [22,25], and US NRC RG 3.54 Rev.2 —citeguide354 implemented as integrated methodology or as stand-alone applications.

The latest set of decay data implemented in the SNF code is based on ENDF/B-VII.1 [274]. With the release of ENDF/B-VIII.0 [211] and following its implementation [275] in CASMO5, a decay data set based on ENDF/B-VIII.0 was generated and available with the SNF code.

### Methodology

The cross-section library, generated by the lattice physics code CASMO5, provides nuclide concentrations tabulated by exposure, moderator density history, control rod history, fuel temperature, and boron history. The nodewise exposure and accumulated history parameters, obtained from qualified operational reactor data or code simulation (SIMULATE5), are used as entry points in the interpolation routines and, with the power history model, are used to compute the end-of-life (EOL) nuclide concentration. The fuel assembly is divided into axial nodes (index: K) that may be further divided into sub-nodes (N). A sub-node is characterized by its associated lattice type (L). The concentration $C_{\text{base}}$ of a given nuclide (J) in sub-node
(N) of node (K) is calculated as:

$$C_{\text{base}}^{J,K,N} = C_{\text{tab}}^{J,L(N)}(E) + \sum_{h \in H} \Delta C_{\text{tab}}^{J,L(N),h}(E) \cdot \Delta H_h, \quad (73)$$

where

1. $C_{\text{tab}}^{J,L(N)}(E)$ is the concentration (g/t) of nuclide J in lattice type L(N) at exposure E;

2. $\Delta C_{\text{tab}}^{J,L(N),h}(E)$ is the history-differential (difference between base and alternate depletion histories) concentration of nuclide J in lattice type L(N) at exposure E;

3. $\Delta H_h$ is a nodal history effect (moderator density, control, boron, void, fuel temperature, etc.).

In this equation, all concentrations are obtained by interpolation in two-dimensional tables for lattice type L to the nodal exposure E. The nodal concentrations, $C_{\text{base}}^{J,K}$, are obtained by a weighted average over the sub-nodes in each node. Further, the individual nuclide concentrations are affected by the nodal power variations, the effect of which can be significant for many nuclides. The SNF code makes use of the nodal power histories (from the 3D core simulator) to apply “power history correction factors”, $C_{\text{PHIST}}^{J,K}$:

$$C_{\text{EOL}}^{J,K} = C_{\text{PHIST}}^{J,K} \cdot C_{\text{base}}^{J,K} \quad (74)$$

to the nodal concentrations. The power history correction factor for nuclide J in node K with burnup E is defined as the ratio of the final concentration, $C_J$, computed with the actual power history, to the reference concentration, $C_t^{\text{ref}}$, computed with the reference power level at the same burnup.

The EOL nuclide concentrations are then applied to the decay chains resolved by the Bateman system. The Bateman system is solved analytically by a method from reference [276] to obtain the final concentration, $C_t^{J,K}$, at decay time t. The decay chains are built exactly as defined by the parent-daughter information from ENDF/B, which allows all other parameters (such as recoverable energies, decay constants, and emission spectra) to be applied implicitly. With the final nuclide concentrations, it is straightforward to compute the disintegration rates at cooling time t, and the recoverable energy is computed by:

$$Q^{J,K}(t) = \sum_{m=\alpha,\beta,\gamma} Q_m^{J,K} \cdot R^{J,K}(t) \quad (75)$$

where

1. $R^{J,K}(t)$ is the disintegration rate of nuclide J in node K, after time t;

2. $Q_m^{J,K}$ is the average recoverable energy per disintegration for particle type m.

Decay data

The SNF library includes basic data, such as decay constants, atomic masses, and nuclide transmutation chains; radiation emission spectra for photons from radioactive

decay, ($\alpha$,n) reactions, bremsstrahlung, and spontaneous fission; electrons and alpha particles from radioactive decay; neutrons from radioactive decay, spontaneous fission, and ($\alpha$,n) reactions; decay heat production; and electro-atomic interaction data for bremsstrahlung production. These data are compiled from fundamental (ENDF/B, ENSDF [277], TENDL [278,279]) and processed (ESTAR and ASTAR [280]) sources. For consistency with CASMO5, ENDF/B-VII.1 data are complemented with data from TENDL-2012 [281], and the ENDF/B-VIII.0 data are complemented with data from TENDL-2017 [282]. The evaluation and validation procedure of the decay data implemented in the SNF library is reported in reference [281].

The basic decay data set per nuclide (half-life, atomic mass, decay modes, and branching ratios, Q-values, and daughter isotopes) is entirely derived from ENDF/B, as are the discrete emission spectra for $\alpha$, $\beta$, and x-rays. The uncertainty of line energies for the discrete spectra is provided by ENDF/B. When populating the photons in the user-defined energy groups, the intensity of a given line may be split equally between the adjacent groups if the group boundary falls within the corresponding line-uncertainty. For the continuum spectra for prompt photons emitted during spontaneous fission and $\alpha$-decay, empirical correlations from reference [282] are adopted.

The ($\alpha$,n) and ($\alpha$,γ) reaction data are processed from TENDL by NJOY-2012 [283], with the alpha particle stopping powers for the elements, including compounds and mixtures, derived from ASTAR. The TENDL data for the oxygen isotopes $^{17}$O and $^{18}$O is complemented by ($\alpha$,n) cross sections from Perry and Wilson [284]. The approach used for these isotopes was to rescale the neutron resonance data of Perry and Wilson to match the TENDL data at the upper energy of the resonance range (5.404 MeV for $^{17}$O and 5.18 MeV for $^{18}$O); the rescaling increased the cross sections below these energies by ≈25% and ≈4%, respectively. This allowed the cross-section representation to be extended to the upper limit of the SNF $\alpha$-decay data (≈12 MeV). The photon production data were not rescaled. Compared to TENDL-2012, the TENDL-2017 based ($\alpha$,n)-yields are on average about 2% smaller.

The continuum spectra for electrons and positrons were calculated by the method of Stamatelatos and England [285], coupled with the Fermi function approximation of Schenter and Vogel [286], with the spectrum endpoint energies provided in the ENDF/B decay sub-library. The external bremsstrahlung cross sections and production kernels were processed from the ENDF/B electron sub-library, while the electron-stopping power data were calculated with the ESTAR application.

7.4 FISPIN

FISPIN is a spent fuel inventory code developed in the UK. Initially developed by UKAEA [287], in the 1970s but subsequently developed by the National Nuclear Laboratory, it simulates the production and destruction of three groups of nuclides; heavy elements (associated with the production and decay chains from natural thorium and
uranium isotopes), fission products and products of activation of naturally occurring elements irradiated in reactors. The solution method is developed from the EXTRA code developed by Sidell [288]. The code uses burnup-dependent cross sections, fission product yields, energy per fission, atomic masses, isotopic abundances, and decay data to solve the production and destruction equations for each nuclide to determine the changes in the composition during irradiation. The code can use steps of irradiation or cooling defined by neutron flux in n/cm²/s and irradiation duration, the rating in MW/t and time to reach required irradiation in MWd/t (aka burnup), a set of fission rates and time, or for no neutron flux/cooling times a duration or cooling time. Additionally, at any step in the irradiation or cooling, the elements can be separated either by removing all but a set of specific elements of which a fraction passes into a product stream during chemical processing, or letting through all but a specific list of elements for which only a fraction is passed through to simulate a waste stream. Various output options exist to list at each step the masses, activities, decay heat, gamma-ray emissions, β activities and energy released, α activities and energy released, and neutron emission from both spontaneous fission and (α,n) reactions for an arbitrarily defined compound. The output can be split by nuclide or element to support chemical engineering applications. FISPIN10A is the latest version of the code and is released with a GUI under Windows or Linux that simplifies production of FISPIN runs. It includes a set of burnup-dependent neutron cross-section libraries for Magnox, AGR, PWR, and BWR fuels with natural enrichment for Magnox, up to 4% enrichment for AGR, and up to 5% for PWR and BWR. The validation for this code is available from reference [289]. FISPIN11 is a new version of the code currently being developed, but not yet released. It uses a purpose-developed Rosenbrock stiff problem integration method and is designed to be called an application via an API so that it can be incorporated into purpose-built wrappers around the FISPIN11 kernel. This simplifies QA of the code by only solving the production/destruction equations in FISPIN11, with data preparation and conversion of the output into user requirements being carried out with the developed wrappers.

## 7.5 EVOLCODE

The EVOLCODE 2.0 system is a burnup code developed and maintained at CIEMAT for more than 20 years. It combines neutron transport and depletion calculations to describe the burnup evolution of either critical reactors or subcritical systems operating in any neutron spectrum regime. The code has been focused on the estimation of a great variety of nuclear reactor parameters, with a particular interest in the nuclide composition evolution of the fuel in a nuclear reactor and also in the neutron transport characterization of the core [290].

The present version of the code allows using any published version of the general MC transport code MCNP/MCNPX (Pelowitz, 2014) for the neutron transport part of the calculation. This code is a Monte Carlo

radiation-transport code able to track, among others, neutrons, photons, and electrons. Complex statistical processes can be modeled without explicitly solving the Boltzmann equations, by means of simulating individual neutrons and then analyzing statistical averages of their interactions with the material medium. The depletion part of EVOLCODE can be performed by one of two codes: the activation code ACAB [291] and the ORIGEN code [292]. Both codes solve the Bateman equation using the matrix exponential method.

The EVOLCODE 2.0 system can to work with any basic database, with the only limitation that the cross sections are formatted in ACE format [293], usable by MCNP, to ensure coherence. Among the data that the code uses from the library, some can be highlighted, such as the cross sections, decay data (including the heat emission per decay), fission product yields, reaction branching ratios, and the probability tables to treat unresolved-resonance self-shielding as MCNP does. The nuclear data library most commonly used in EVOLCODE is the JEFF-3.3 nuclear data library [58].

## 7.6 VESTA

VESTA [294] is a Monte Carlo depletion code developed by IRSN, which couples a continuous energy Monte Carlo neutron transport code with a depletion module. VESTA validation calculations were performed using the JEFF-3.1, JEFF-3.2, ENDF/B-VII.0, and ENDF/B-VII.1 nuclear data libraries. VESTA 2.1.5 calculations were performed using MCNPX 2.6.0 and PHOENIX codes, with the JEFF-3.1 and ENDF/B-VII.0 nuclear data libraries. VESTA 2.2 calculations were performed using MCNP6.1 and PHOENIX codes, with the JEFF-3.2 and ENDF/B-VII.1 nuclear data libraries.

These libraries include the required cross-section data for the transport and reaction rate calculation but also the neutron induced fission yield data, isomeric branching ratio data, and decay data. The ORIGEN 2.2 library used as the base library is PWRU50.LIB. This base library is used only as a source for cross-section and isomeric branching ratio data when the information is missing in the nuclear data library.

No time-dependent temperature evolution has been modeled and an average boron concentration during irradiation is considered. The irradiation history is subdivided into depletion steps of 1 GWd/tHM, according to the VESTA validation procedure (see Ref. [295]). The number of particles per cycle and the number of active cycles for the MCNP calculation are determined to ensure statistical errors in the neutron flux of less than 0.1%.

## 7.7 Serpent

Serpent is a multi-purpose three-dimensional continuous-energy Monte Carlo particle transport code [296]. Serpent has been developed at VTT Finland since 2004. The original motivation for Serpent development was to provide spatial homogenization based on the high-accuracy Monte
Carlo method [297]. However, since the early years, Serpent's applications have been greatly expanded from reactor physics to fusion applications, photon transport, and multi-physics modeling. Currently, Serpent development is carried out in the context of the development of a new computational framework for coupled core physics calculations, Kraken [298].

The geometry routine in Serpent is based on a three-dimensional constructive solid geometry (CSG) model. The available selection of elementary quadratic and derived surface types enables the modeling of a large scale of different geometries and reactor types. CAD-based geometry type is also available for the modeling of complicated and irregular systems [299]. Particle tracking in Serpent is handled through a combination of conventional ray-tracing-based surface tracking and the rejection sampling-based delta-tracking method. Surface tracking is used in the presence of heavy absorbers when the efficiency of delta-tracking rejection sampling becomes low.

Serpent includes a built-in burnup calculation routine utilizing the CRAM method (Chebyshev Rational Approximation) to solve the Bateman depletion equations [300]. In a typical burnup calculation, a total number of 1200–1600 nuclide concentrations are tracked including 250–300 actinides and fission products with cross sections. The output includes depletion zone-wise, material-wise, and total nuclide concentrations, masses, activities, decay heat, and other parameters. The decay heat of a given material is based on nuclide concentrations, decay constants, and total energy emitted in a decay. The latter two components are read in Serpent from ACE format nuclear data libraries. Currently at VTT, nuclear data based on JEFF-3.2, JEFF-3.1.1, ENDF/B-VII.1, and JENDL-4.0 are available.

Serpent is available for academic users through the OECD/NEA databank and Radiation Safety Information Computational Center (RSICC). Commercial licenses are issued by VTT.

## 7.8 Codes and libraries developed by JAEA

Japan Atomic Energy Agency (JAEA) has developed some depletion calculation codes, among which two major codes: MVP-BURN [301] and SWAT4 [302] are briefly introduced in the following.

MVP-BURN is a depletion calculation code based on Bateman's method coupled with the continuous-energy Monte Carlo code MVP [303] for neutron flux and effective cross section calculation. MVP-BURN is now available as a part of the MVP3 package, which is the latest version of MVP. As latest nuclear data libraries, cross-section data libraries based on JENDL-4.0, ENDF/B-VII.1 and JEFF-3.2, and a burnup chain data library based on JENDL-4.0 are available. Recently, MVP-BURN (MVP3) was used in the blind benchmark exercise for spent nuclear fuel decay heat and its performance was compared [15].

SWAT4 is a depletion calculation code system driving ORIGEN2 [292] for depletion calculation and three neutron transportation calculation codes: the deterministic code SRAC [304], and the continuous energy Monte Carlo

codes MVP and MCNP. Users select a neutron transportation solver from these three codes depending on their purpose. SWAT4 works as an interface between ORIGEN2 and the selected neutron transportation solver, and various calculation capabilities of ORIGEN2 are available. Cross-section data of SWAT4 depends on a cross-section library of the selected neutron transportation code. Other nuclear data such as radioactive decay and fission yield data is applied from basic ORIGEN2 libraries or other libraries for ORIGEN2 such as ORLIBJ40 [305], which is a nuclear data library set based on JENDL-4.0 for ORIGEN2.

Regarding nuclear data libraries, JENDL-5, the latest version of the JENDL series, provides comprehensive evaluated data of nuclear reactions and decays, thermal neutron scattering kernel, fission product yields, etc. to facilitate nuclear science and technology covering a variety of applications not only for nuclear energy but also for radiation-related fields. The number of nuclides of neutron reaction data increased to 795 which is almost double of 406 of JENDL-4.0 to cover various applications. New evaluations were carried out for light nuclides, actinides, structure materials, fission products, minor-actinides using up-to-date experimental data and theoretical tools.

Thermal scattering kernel of light water, heavy water, etc. were also newly evaluated. The benchmark tests mainly for light-water and fast-neutron reactors and neutron shielding were carried out, showing globally a better performance to JENDL-4.0. New evaluations were also carried out for fission product yield and decay data, which are of particular importance for decay heat estimations.

The characteristic point of the evaluation of JENDL-5 fission product yields is the introduction of theoretically calculated nuclear shell corrections in the independent fission product yields and of the Hauser-Feshbach statistical model in calculating the fraction of isomer states [306]. The evaluated data of decay data is mainly based on ENSDF. In addition, newly evaluated data of delayed neutron emission and spectra were added, based on studies from the IAEA CRP [307], and theoretically calculated particle emission spectra.

The JENDL-5 fission product yield and decay data are benchmarked by using PIE data of Mihama-3 and Takahama-3, decay heat data of instant neutron irradiation carried out at the Oak Ridge National Laboratory, Yayoi fast reactor facility (University of Tokyo), delayed neutron yield data, and anti-neutrino spectra data [212,306].

Good improvement for delayed neutron yield evaluations was therefore obtained [212]. Verifications by comparing the calculated radioactivity and photon intensity from the debris of the Fukushima-Daiichi nuclear power plant were carried out, including calculations with the ENDF/B-VIII.0 and JEFF-3.3 libraries. One of the typical results is shown in Figure 36, where the fuel compositions are taken from reference [308]. It can be seen that both radioactivity and photon intensity of JENDL/DDF-2015 show smaller values than other libraries. This is due to the lack of several nuclear data, which are the isomer state of $^{235}$U, x-rays of $^{137m}$Ba, γ-rays of $^{241}$Am, etc.
<figure><img src="image_39.png" /><figcaption>**Fig. 36.** Radioactivity (left) and photon intensity (right) from debris of Fukushima-Daiichi nuclear power plant estimated from the decay data of JENDL-5, JENDL/DDF-2015, JEFF-3.3, and ENDF/B-VIII.0 (data are provided by courtesy of T. Matsumura (CLADS, JAEA))</figcaption></figure>

Those issues are fixed for JENDL-5 and results are much closer to ENDF/B-VIII.0 and JEFF-3.3.

## 7.9 COMPASS

COMPASS is a burnup analysis code system developed by CRIEPI for the prediction of neutron transport and nuclide composition in large-scale geometry such as multi-assembly or whole core geometry [309–312]. COMPASS uses a general-purpose continuous energy Monte Carlo code MVP3 [303] as neutron transport calculation solver, the CRAM method [300] and burnup chain-data based on JENDL-4.0 (ChainJ40) [313] in the depletion solver for nuclide composition calculation. All nuclear data libraries such as JENDL-4.0, ENDF/B-VII.1, and JEFF-3.2 bundled in MVP3 are available in COMPASS. A function of unified geometry treatment is implemented in COMPASS for coupled calculations with thermal-hydraulic and fuel performance calculation codes. One of the advantages of COMPASS is flexibility in setting calculation conditions and parameters. For example, COMPASS can change material temperature and geometry parameters in each calculation step and use all options of MVP3 code (e.g. exact resonance elastic scattering model [314]). The depletion solver in COMPASS has been validated enough by Post Irradiation Examination (PIE) analysis for PWR and BWR 55 assemblies.

## 7.10 ACAB

The ACAB code is a computer program designed to perform activation and transmutation calculations for nuclear applications [224]. ACAB considers decay transitions that proceed from the ground, first, and second isomeric states. All the nuclear reaction processes and all the corresponding reaction products are treated in ACAB code.

ACAB's primary result is nuclide concentrations as a function of time for each spatial interval or zone defined in the system. From the nuclide concentrations, ACAB

is able to generate radionuclide activities, decay heat (total and contributions from the different types of radiation), decay gamma spectra, contact dose rates, waste disposal ratings, offsite doses to the most exposed individual, as well as collective doses and associated consequences, radiotoxicity (committed effective dose equivalent by inhalation and ingestion) and neutron emission among other quantities that depend upon nuclide concentrations.

In addition to the calculational facilities, ACAB also performs sensitivity analysis of the activation results. Critical radionuclides are identified and pathways contributing to their production are evaluated. ACAB is able to compute uncertainties both on inventory calculations and inventory-related quantities. The method is based on the application of the Monte Carlo technique, it allows dealing efficiently with the synergic/global effect of the uncertainties of the total set of cross sections to obtain the overall uncertainty on the radiological calculations.

The potential of ACAB for fission reactor applications has been proven to be reliable in computing accurate nuclide inventory. Examples of results using ACAB coupled to a Monte Carlo neutron transport code are: a simple UO₂ pin-cell benchmark [315], an ATRIUM-10XP design BWR fuel assembly [316], a high temperature gas-cooled reactor plutonium burnup benchmark [317]. In addition, ACAB is able to predict decay heat and delayed neutron emission rate, average neutron energy and neutron delayed spectra after a neutron fission pulse [318].

## 8 Validation

Validation is an essential step to gain confidence in simulation results. It is the link between experiment and calculation. The latter uses the former to derive prediction power, biases, and uncertainties. The former is adapted based on the needs of the latter. Decay heat validation is a field that cannot be extensively be described in this report, given the large amount of related publications; only an
overview, with some relevant references is provided in the following.

Not all codes and nuclear data are included, and the focus is here on decay heat. It is rather common that the same codes are also validated for other quantities, such as $k_{\text{eff}}$, nuclide concentrations, or radiation emissions (gamma, beta, alpha, neutron) in SNF. Such validation will be partially referred to.

## 8.1 Definitions

A number of terms are often used together, such as verification, quality assurance, validation, and uncertainty quantification. The definition of these processes can vary from one organization to the other, but global notions are generally common. Detailed examples can be found in references [319–321].

Verification consists in the actions that demonstrate that a specific code (or its parts) is performing as expected, given the initial requirements from the design phase. It can for instance be a well-defined formalized process having the goal of determining if specific equations are correctly solved from the numerical and software aspects.

Quality Assurance (QA) refers to the process of determining and documenting whether or not items, processes, or documents which support the computer code conform to specified requirements [322].

Validation refers to the process of evaluating a computer code after the completion of its development and prior to its application, to ensure compliance with the requirements of the intended application. Validation provides the evidence that the computer code is fit for the intended purpose, by comparison of the computed results with data from experiments or other trusted sources [322].

Quite often, verification and validation of a computational tool go hand in hand and are referred to as V and V. Herein validation, the focus of this section, refers to the comparison of computed results with data from experiments. Uncertainty quantification, which is an important part of the validation process, is discussed in Sections 9 and 6.6.

An important aspect of the decay heat validation is assessing the effect of nuclear data used with a computational tool on the comparison calculation experiment. For the same computational tool, results obtained with different sets of nuclear data (cross sections, fission yield, decay data), could be significantly different. When comparing the performance of various computational tools in predicting decay heat, it is recommended for completeness to cite both the computational tool and associated nuclear data used for the simulation, rather than the name of the computational tool only.

## 8.2 Validation based on direct decay heat measurements

The measurement data serving as the basis of the validation studies summarized here were obtained from decay heat integral experiments (full assembly or fuel rods

within an assembly), or from pulse fission experiments for fissile materials (foils). Full-assembly decay heat experiments are generally feasible and performed for cooling times greater than 2–3 years; these cooling times are of particular importance to spent nuclear fuel storage, transportation, and disposal applications. Full-assembly decay heat experiments are not feasible for cooling times up to 105 s that are of interest for evaluation of postulated loss-of-coolant accident scenarios; for this time range, pulse fission experiments serve for determining the energy released following reactor shutdown to derive decay heat. These experiments generally involve gamma and beta spectroscopy measurements to derive decay heat from the individual gamma and beta energy components.

### 8.2.1 Decay heat validation of ALEPH2

ALEPH2 has been extensively validated based on various sets of experimental data relevant to spent nuclear fuel applications, including for nuclide inventories [323–325] and decay heat. For LWR spent nuclear fuel, ALEPH2 was validated using based on full-assembly decay heat measurements performed at the Clab facility for 25 PWR and 34 BWR spent fuel assemblies. Simulation results obtained by running ALEPH2 with various neutron transport and fission yields showed the calculated decay heat values underestimate the measurement data by less than 4% for PWR assemblies and less than 6% for BWR assemblies [326].

It is also worth mentioning that additional validation of ALEPH2 for decay heat calculations for fast reactor spent nuclear fuel was performed, using limited information available from the JOYO MK-II core experimental program. The obtained results showed an underestimation of less than 7% [327].

### 8.2.2 Decay heat validation of VESTA 2.2

VESTA 2.2.0 has been extensively validated for decay heat calculations for intermediate and long cooling times applications, based on full-assembly decay heat measurements for 52 PWR UO₂ assemblies. These 52 assembly measurements include: 34 measurements performed by the Clab program, 14 measurements from the General Electric program, and 4 measurements from the HEDL program.

Using the JEFF-3.1, JEFF-3.2, ENDF/B-VII.0 and ENDF/B-VII.1 nuclear data libraries, VESTA calculates the assembly decay heat for almost all cases within 4% of the measured decay heat.

On average, the calculated decay heat values that used JEFF-3.1 appear to result in a systematic underestimation of about 2%, while those that used ENDF/B-VII.0 led to a systematic underestimation of about 0.5% with almost all points within the two sigma experimental uncertainty. While the results for the comparison calculation-measurement for the GE experimental data are similar to those obtained for the Clab data, it appears that some specific measurements from the HEDL data are problematic; specifically, this concern an unrealistic large uncertainty on the final discharge burnup for two assemblies, which will require further investigation. When omitting
these two assemblies from the HEDL data set, the average calculated-to-experimental decay heat for the aggregate set of the Clab, GE and HEDL measurements is 0.979 for JEFF-3.1 and 0.996 for ENDF/B-VII.0 calculations.

The validation calculations performed with the JEFF-3.2 and ENDF/B-VII.1 libraries considered only the Clab measurements. Using the JEFF-3.2 nuclear data library, except for a few outliers, VESTA calculates the assembly decay heat for almost all cases within 3% of the measured decay heat. On average, the JEFF-3.2 calculated decay heat values underestimate the measured data by about 2%. Using the ENDF/B-VII.1 nuclear data library, all VESTA results overestimate the measured data by about 0.5%, with almost all points being within the two sigma of the reported experimental uncertainty.

### 8.2.3 Validation of DARWIN/PEPIN2 and MENDEL

Validation tests should be the responsibility of another group of persons (i.e., not the developers). Generally, such tests are realized by the developers and integrated in the non-regression base. Nevertheless, the full validation of the code system needs a group of users able to use the code differently from what the developers were initially thinking of. They should have access to a large data bank of experimental data, and use it for the code validation. At the CEA, this work is mainly done by the group elaborating on the DARWIN2 package. Validation tests consist of PIE, activation experiments, decay heat measurements after fission pulse or irradiation, and the scope tends to be as complete as possible, to fully validate the model (computer code and nuclear data).

### 8.2.4 Decay heat validation of CASMO5 and SNF

The validation effort with the deterministic transport code CASMO5, regarding integral measurements of assembly decay heat, can be found in a number of publications, for instance, references [13,14,34]. In parallel, indirect validation was also performed using Post Irradiation Examination and measured nuclide concentrations, such as in references [33,34,234,245,246]. Different nuclear data libraries were used, indicating spread of data similar to other codes. In parallel, the impact of the nuclear data uncertainties was also calculated using the irradiation histories from either the Clab assemblies, or ARIANE PIE data (see mentioned references).

Regarding the SNF code, it calculates nuclide concentrations and decay heat (among other quantities) either based on irradiation histories specified by the user or based on core follow-up simulations from the SIMULATE-3 or SIMULATE5 codes. As for CASMO5, it was validated with PIE nuclide concentrations [245,328,329], and with Clab decay heat measurements (including different nuclear data libraries) [330].

### 8.2.5 Decay heat validation of SCALE

Validation of computational capabilities in the SCALE code system for decay heat prediction has been documented for various versions of SCALE, starting with 5.1 [266] and continuing with 6.1 [79], 6.2 [337], [269], and

6.3 [11]. This continuing validation effort has involved the use over the years of different ENDF nuclear data libraries, from ENDF/B-V through ENDF/B-VIII.0.

The most recent validation for decay heat at cooling times of relevance to transportation, storage, and disposal of spent nuclear fuel was performed based on 91 PWR and 145 BWR full-assembly decay heat measurements. These measurements cover a cooling time ranging from 2 to 27 years for PWRs and BWRs, a burnup range of 18-51 GWd/MTU for PWRs and a burnup range of 5-45 GWd/MTU for BWRs. The decay heat predicted with SCALE 6.2 and ENDF/B-VII.1 libraries is on average within 1% of the experiment results ($\sigma = 1.6\%$) for PWRs, and within 2% of the experiment results $\sigma = 7.7\%$) for BWRs, see reference [269]. The validation study also looked at the effect of using a different nuclear data library. The C/E results obtained with SCALE 6.2.4/ENDF/B-VII.1 and those obtained with SCALE 6.1/ENDF/B-VII.0 when using a same basis of 121 LWR assembly measurements differed on average by 0.6% ($\sigma = 0.5\%$) for PWR assemblies and by 1.2% ($\sigma = 0.8\%$) for BWR assemblies. This difference was found to be primarily due to the change in the $^{238}$Pu capture cross-section between ENDF/B-VII.0 and ENDF/B-VII.1 [269].

Validation of SCALE 6.2 with ENDF/B-VII.1 libraries for very short cooling times was based on data from pulse fission experiments designed to determine the energy release from fission of $^{233}$U, $^{235}$U, $^{238}$U, $^{239}$Pu, $^{241}$Pu, and $^{232}$Th. Overall, this validation basis covered cooling times from 0.2 s to 11 h following fission. The results obtained showed generally good agreement between calculated and measured values, within the level of the reported measurement uncertainties [269].

## 8.3 Indirect validation

As mentioned, a number of transport and depletion codes are also validated for nuclide concentrations, based on Post Irradiation Experiments as for instance described in the SFCOMPO database [32], or with criticality-safety benchmarks, as from the ICSBEP database. In the latter case, the fuel under consideration is fresh and does not directly overlap with the SNF characteristics. Such validation is nevertheless useful for decay heat as it demonstrates the code (and nuclear data) performances for neutron transport. In the case of PIE validation, the prediction of nuclide concentrations after irradiation and cooling periods defines the code predictions for calculating the fission product and actinide concentrations, which play a key role in the calculation of decay heat. In both cases, such indirect validation for decay heat is also of prime importance and demonstrates the complexity of a full validation scheme. A number of publications can be found for PIE validations, with codes also used in decay heat calculations, for instance in references [245,328,331].

Another source of indirect validation for SNF decay heat is the calculation of pulse irradiations and comparisons with measurements, as presented in Section 6.4.4. The gamma-ray and beta-electron emissions on different actinide targets were measured at short cooling time
(down to 1 second), with well-characterized neutron beam, which represent the only available experimental data for a short time range. An example of such decay heat validation can be found in reference [186].

# 9 Recommendations

Based on the current understanding of the physical phenomena at stakes in the decay heat for current spent nuclear fuel assemblies, the following recommendations are proposed. This list is not comprehensive, yet represents the main priorities recognized at the time of this report.

## 9.1 Extending the current experimental database

As mentioned in Section 5, only one full assembly calorimeter is nowadays in activity worldwide. The Clab calorimeter is a unique facility, dedicated to measuring assembly decay heat coming from Swedish nuclear power plants. The high-quality published measurements, with the necessary information on the assemblies, allow for a number of code validations. Most of the measured values come from reference [76], complemented with information from ORNL reports. They correspond to assemblies irradiated up to the beginning of the year 2000, with initial enrichments and end-of-life burnup lower than today's typical values. Additionally, these assemblies contained UO₂ fuel solely.

A new experimental program is currently being analyzed with the support of SKB and EPRI, expanding the range of burnup and initial enrichment to higher values (possibly up to 55 MWd/kgHM burnup and 4.1% $^{235}$U enrichment). Such results, when available will help in covering part of the present gap with respect to the SNF characteristics as employed nowadays. Assemblies currently used in reactor cycles generally have higher initial enrichment (close or slightly higher than 5.0%) and are reaching average burnup values higher than 60 MWd/kgHM. The cooling times of interest, in the case of deep geological repository are also different than the ones measured (from 2 to 27 years), mainly due to the delays in the start of such facilities (as well as encapsulation facilities). Finally, the magnitude of the measured decay heat rate is relatively restricted: the Clab and GE-Morris measurements correspond to values between 20 and 940 W, with the addition of the recent measurements on five assemblies, from 660 to 1660 W [15]. Most of the measured values are in the range of 150 to 500 W, which is rather limited for a complete validation work.

There is therefore a clear need for additional integral measurements on SNF with higher burnup, higher initial enrichment, higher decay heat rate values, and different cooling times, ideally at different facilities. Another pressing need concerns MOX assemblies. There are currently no integral measurements in the public literature on MOX SNF. Such assemblies have been used in more than 40 reactors worldwide and a number of such SNFs are currently in storage. The lack of decay heat measurements

for this type of assemblie is limiting the current validation capabilities. Finally, other types of SNF do not have measured decay heat: VVER and CANDU assemblies and more recently ATF assemblies.

In parallel to integral decay heat measurements, fission pulse decay measurements also bring valuable information, as presented in Section 6.4.4, helping to disentangle different components of the total decay heat, with direct information for the underlying nuclear data (fission yields and decay data). The compilations used for validations, mainly references [44,107,115,125,126] for thermal neutron-induced fission, are based on measurement methods and reported uncertainties which can also be nowadays improved. Better estimation of SNF decay heat will therefore benefit from new fission pulse decay measurements based on modern techniques.

## 9.2 Improving theoretical understanding and nuclear data

As presented in Section 6, various possibilities exist for the prediction of the SNF decay heat (e.g. simple formula, summation methods, standard methods). Some of them rely on a physical description of decay phenomena and therefore depend on the correct estimation (or simply the availability) of specific nuclear data (half-lives, decay energies, decay modes). In the case of heavy actinides produced from many subsequent neutron captures, such as curium isotopes, a number of fission and capture cross sections are also of importance. For fission products, both independent fission yields and capture cross-sections can affect nuclide concentrations and consequently the energy released during the sum of their decay. Improving the knowledge of such nuclear data will impact the prediction capabilities for nuclide concentrations and consequently for decay heat. It is therefore recommended to improve the quality of nuclear data for nuclides of importance for decay heat.

In parallel, the handling of such decay data can be mathematically challenging, for instance, due to the large difference in the order of magnitude of half-lives. Different numerical solvers are used such as the matrix exponential method [220], the Ordinary Differential Equations (ODE) [97], the CRAM matrix exponential method [296], or the transmutation trajectory analysis (TTA) [215]. The implementation of such methods (even for the same method) can lead to differences, affecting calculated quantities such as the decay heat. Benchmarking these methods and their implementations can help in improving the precision of calculations, effectively removing biases due to numerical methods.

Finally, physical phenomena might need to be better taken into account, such as delayed fission. After a fission event, a number of delayed neutrons are emitted, up to minutes following the prompt release of energy (for instance from the decay of $^{87}$Br). Such neutrons can also generate fission events. It was shown in Section 6.2 that these delayed fission events can be of significant fraction of the total decay heat up to one or two minutes after a reactor shutdown. The question of how well is this process taken into account by simulation codes is of relevance for
short-term cooling after a reactor transient, as neglecting it would seriously undermine the estimation of the core decay heat.

## 9.3 Best estimate calculations and conservative estimations

As mentioned in Section 6, the prediction of decay heat can be achieved with detailed calculations based on summation methods, but also with standard methods such as ISO, ANSI, or DIN [22,26,27]. Such standard methods have the advantage to provide fast answers, with a certain amount of conservatism. Still, some of them require fission fractions as a function of irradiation history (for the main four actinides contributing to fission events), and these quantities are often derived from dedicated transport calculations. It is recommended to use more precise methods for final assessments, performing dedicated transport and depletion calculations based on specified irradiation history and assembly characteristics. If an amount of conservatism is still required, it is preferable to explicitly add it as a separate quantity.

## 9.4 Biases and uncertainties

Biases and uncertainties are as important to estimate as the decay heat values themselves. They depend on the input data and assumptions, as well as on model simplifications and calculation methods. Needless to say biases also depend on the reference data, their uncertainties, and correlations. As many possibilities can be selected, it is of prime importance to indicate the maximum amount of information regarding the performed calculations: code versions, nuclear data libraries, calculation methods and experimental database.

Regarding calculated uncertainties, users are generally facing many choices for the assumptions regarding input values, such as the recommendations provided in reference [232], various uncertainties found in the open literature or in-house values. There are for instance currently no recommendations on fission yield correlations, or decay data correlations, which directly impact calculated uncertainties on decay heat. It is not uncommon to find differences larger than factor 2 or 3 for the uncertainties on the same calculated quantities, such as nuclide concentrations. Similarly, operational and technological parameters are not systematically recorded with recommended uncertainties (and probability density functions).

It is therefore of important to provide the user community with systematic recommendations for uncertainties, correlations and probability density functions for all parameters and physical quantities used in simulations. Additional recommendations shall cover sampling and uncertainty propagation methods.

## 9.5 Other recommendations

The list of recommendations can be extended to many other features, most not only applicable to decay heat

calculations but solely linked to simulation methods. Such aspects include but are not limited to:

- **Redundancy.** During the characterization of the SNF for various practical applications (safeguard, storage, handling, final disposal), redundancy of calculations, measurements, or verification of any nature is a best practice in order to identify potential errors or mistakes, for instance from human nature. In the example of selecting SNFs for canister loading, a redundancy can consist of comparing the declared SNF burnup (from plant data) with gamma-ray or neutron emission measurements at the time of loading. Unexpected differences may indicate potential mishandling at a given stage of identification. In the case of decay heat calculation, predictions by independent parties can help identify various types of mistakes.

- **Comparison of code implementation.** It is known that different coding implementations of the same mathematical model can lead to differences, sometimes beyond numerical precision. Comparing modules of codes, which are expected to perform similar tasks, can be useful to improve speed and efficiency, and detect potential errors.

- **Blind benchmarks.** The term of "blind benchmark" is understood here as the exercise of performing similar simulations between different (independent) partners, without accessing any information concerning the calculations of other participants, or measurements to be later compared to, for the full duration of the project. The main advantage of this approach is to avoid having interference effects between participants themselves and between simulations and measurements. The consequence is a true independence of calculations, performed with a common set of information. If such strict conditions are accepted, different iterations are possible, while keeping track of changes: a second step can be performed after a limited comparison between calculated quantities, or with measurements. The evolution of calculated values as a function of iteration can be an indication of improved knowledge, which can be clearly quantified. The consequences can be interpreted as simulating the realistic situation where no experimental will be available, and other calculations to be compared to. It is therefore recommended to perform new benchmarks as blind, helping to better assess results obtained in realistic assignments.

- **User effect.** The user effect can lead to non-negligible variations of calculated or measured quantities. It can happen due to different interpretations of the same information, ambiguity of definitions, selections of different methodological options (not described in a task definition), or simple mistakes. It is often difficult to quantify, in part due to the absence of blind benchmarks, or to the absence of redundancy. It was nevertheless observed in different situations and it is necessary to include a measure of this effect in dedicated best estimate calculations, with uncertainties. Missing to recognize it would lead to an overestimation of the total uncertainty.
Some of these four items are often partially addressed within institutes, using verification procedures, so-called logbooks, archiving of data, often regulated by ISO or national standards. The validation of the adequacy of these steps can nevertheless be tested within the international community.

## 10 Conclusion

The activities performed in this subgroup 12 (SG12), started in 2022 under the auspices of the NEA's WPNCS, aimed at describing the current knowledge on decay heat (or thermal decay power) for existing SNF assemblies. A number of different areas were reviewed: theoretical description, simulation, standard methods, and measurements. One of the motivations for this work is the current needs of the nuclear industry for the estimation of decay heat during short-term cooling (such as in core transient) and in the back-end of the fuel cycle (SNF wet and dry storage, up to and including final repository).

As highlighted in this report, the evolution of fuel characteristics for current LWR over the past decades is following a trend of higher initial enrichment, higher average assembly burnup, and extended storage. It is foreseen that in the future, the initial fuel enrichment will commonly go beyond 5% in $^{235}$U, also leading to higher average burnup values. The need for code and method validation (as well as for nuclear data libraries) is as acute as in the past, and the current experimental databases, used for validation, are not covering all SNF characteristics (neither fuel type, e.g. MOX, CANDU or VVER fuel, or ATF fuel). It was therefore recognized in this report that validation efforts need to continue, with a number of consequences in the areas mentioned above.

Experimentally, it was acknowledged that the existing set of decay heat measurements is highly valuable yet limited, as only one calorimeter is today in use (at SKB, Sweden). Such limited availability has two main consequences: (1) all recent (and eventually new) measurements come from a unique facility, rendering redundancy and independent checking impossible, and (2) all recent (and eventually new) measurements are performed on a limited set of SNF, corresponding to the assemblies in use in Sweden.

From the simulation aspect, new needs related to existing SNF assemblies mainly arise from stringent economic requirements and ever stricter safety rules related to core transient, extended storage, and deep geological repositories. Calculations of SNF decay heat need to be accurately described in terms of biases, outliers, confidence intervals, or uncertainties. When possible, direct validation (mostly for cooling time less than a few decades) is still necessary, but for time periods of interest for deep geological repositories, other questions may come to light. Generally, it was acknowledged that inputs used in simulations (e.g. fuel and moderator temperatures, void, irradiation history) can also be uncertain, leading to uncertainties in calculated decay heat values. New needs for the understanding of important parameters were also recognized, such as assembly average or local burnup, and the impact

of the knowledge of detailed irradiation history. Consequently, the improvements in calculated decay heat values pushed for a better knowledge of experimental methods and uncertainties, a better assessment of key parameters for simulations, and also for making use of various sources of information not directly linked to decay heat, such as Post Irradiation Experiments and full core simulations.

As mentioned, the SG12 brought to light a number of areas where additional efforts will be beneficiary and recommendations are provided in Section 9. As decay heat is one of many characteristics of SNF assemblies, the SG12 activities were also closely connected with the needs in criticality-safety studies, as well as for safeguard requirements. Other international activities were progressing in parallel with common subjects of interest, such as the European Joint Programme on Radioactive Waste Management (EURAD), and the International Atomic Energy Agency Coordinated Research Project on Spent Fuel Characterization (T13018). In the United States, it was also recognized that the existing experimental database needed to grow, and at the time of this subgroup, a report from EPRI regarding new experimental campaigns using the SKB calorimeter was being finalized. The latter will provide new measured decay heat values for SNF assemblies generally having different cooling times, higher burnup, and initial enrichment compared to the currently available set. Finally, during the course of the SG12, two projects for new calorimeters were presented, one in Switzerland at the Gösgen power plant, and one in France, developed by the CEA and EDF. If realized, such projects will address many of the existing issues, regarding the need for measurements on MOX fuel, on a variety of different assembly designs, and on redundancy.

The largest success of this work was to bring together actors from the nuclear industry, research, and safety organizations, having possibly contrasting horizons, requirements and needs. Experiences were shared in constructive ways and a consensus was found to deliver a significant message: better estimations of decay heat are needed.

## Acknowledgements

The authors and contributors of this work are very grateful for the support of the NEA WPNCS.

## Funding

This research did not receive any specific funding.

## Conflict of interests

The authors declare that they have no competing interests to report.

## Data availability statement

This article has no associated data generated.

## Author contribution statement

All authors contributed equally to the paper.

## Appendix A Details on the PUIREX experiment

PHENIX is an experimental reactor as well as an electricity production plan. It started in 1973 with a 565 MWth
<figure><img src="image_40.png" /><figcaption>Fig. A.1. Schematic view of the PHENIX reactor.</figcaption></figure>

power which was reduced to 350 MWth in the 1990s. PHENIX is a pool-type sodium-cooled fast reactor. The primary circuit is composed of three pumps and four inter-mEDIATE sodium-sodium heat exchangers connected to two secondary circuits (see Fig. A.1).

It has a purification circuit with cold traps and a decay heat emergency system which consists of a water loop around the primary pool. Each secondary circuit includes a pump, a steam generator, and a purification circuit. Under nominal conditions, the sodium temperature in the primary circuit is 520°C, the rotational frequency of the primary pumps is 540 rpm and the speed of the secondary pumps is 700 rpm. The purification circuit is operated with its cold traps. The tertiary circuit is operated as well. Just before the test, the reactor was operated under nom-inal conditions but cold traps were by-passed. Then, a quick shutdown was carried out, and the sodium temper-ature plummeted from 520°C to 368°C. An hour after the shutdown, the tertiary circuits were drained. The speed of the primary pumps was reduced from 540 to 300 rpm and down to 110 rpm after two days of cooling. Both secondary circuits were operated at the beginning of the experiment but one was emptied after five days of cooling. Those two events (after 2 and 5 days of cooling) led to discontinu-ities in the temperature evolution. One and a half hours after the shutdown, the reactor was considered to be in isothermal conditions and behaved as a calorimeter.

In this test, the measured quantity is the sodium temper-ature in the primary circuit. The decay heat of the entire core is derived from the measurement through the energy balance of equation (A.1) where DH stands for decay heat, P$_{losses}$ is the inertia depending on the temperature of the sodium θ, P$_{sources}$ are the sources of power and P$_{losses}$ the losses of power.

$$DH = I(\theta) \frac{d\theta}{dt} - P_{\text{sources}} + P_{\text{losses}}. \quad (A.1)$$

The sources and losses terms are derived from a model which depends on the temperature. Knowing the masses of the studied system, the inertia can be computed following equation (A.2) where M$_{i}$ is the mass of the i$^{th}$ reactor component and C$_{p}$$_{i}$ is its heat capacity.

$$I(\theta) = \sum_{i} M_{i} \cdot C_{p_{i}}. \quad (A.2)$$

<figure><img src="image_41.png" /><figcaption>Fig. A.2. Comparison of calculated and measured decay heat [332].</figcaption></figure>

The masses contributing to the total inertia I(θ) are the sodium in the primary circuit (800 tons) and in the sec-ondary circuit (200 tons), the fuel (16 tons), the steel (980 tons) in the sub-assemblies and pins, vessel, structures, pipes, pumps, exchangers, etc. and the graphite in the core neutron shielding (75 tons).

The sources of power are the heat produced by the pri-mary and secondary pumps. The losses of power in the form of radiation and conduction phenomena in the pipes are due to the decay heat emergency system and reactor slab in the primary circuit, and to residual losses in the steam generator in the secondary circuit.

The entire core was equipped with thermocouples in order to measure the sodium temperature and its homogeneity and to check the thermal behavior of the system.

The uncertainty associated with the decay heat derived from equation (A.1) was determined by deterministic propagation of the input parameters uncertainties: ±5% in the calculation of the inertia, 5% due to the pumps, 10% set for the losses in the primary and secondary cir-cuits. The uncertainty associated with the thermocouples is equal to 1.4°C and the uncertainty of the rate of change of the temperature is calculated from the dispersion values of the 13 thermocouples and is equal to 0.5%. The final uncertainty to be taken on the experimental decay heat is between 5 and 11% at one standard deviation depending on the cooling time. The major contributors are the losses and the sodium inertia.

The previous decay heat measurements on PHENIX and SUPERPHENIX concerned cooling times that do not exceed 3 days, and they were all interpreted with older generations of codes and libraries between 1979 and 2002. The PUIREX experiment performed in 2008 was inter-preted with the ERANOS-DARWIN2.3 package and the JEFF-3.1.1 library in 2012 [332]. The comparison of the decay heat values for the PUIREX experiment derived from the measurements to the calculated decay heat with the ERANOS-DARWIN2.3 package and the JEFF-3.1.1 library is shown in Figure A.2.

The decay heat of the PHENIX reactor is predicted within a ±8% discrepancy between 1.5 hours and 12 days. These discrepancies are below the experimental uncer-tainty except for the first day. Figure A.2 also suggests a systematic bias as the discrepancy decreases steadily.
The fuel is the main contributor to the total decay heat, compared to the sodium or steel structures in the vessel. Fission products represent more than 60% of the decay heat, heavy nuclides about 10% around 12 days of cooling (curium and $^{239}$Np contributions), and 30% around 1 day ($^{239}$U and $^{239}$Np). The activated steel structures contribute to less than 10% of the total decay heat. The contribution of the internal storage accounting for assemblies placed at the core edges in order to be cooled before unloaded does not exceed 5%. The sodium contribution is negligible.

# Appendix B Details on decay heat formula

## B.1. Borst-Wheeler function, circa 1946

The Borst-Wheeler function is an approximate and useful expression providing the total power emitted by fission products, given a number of assumptions. It is also called the Wigner-Way formula, based on rough approximations [38]. Such expression can for instance be found explicitly in reference [2]. It provides the total power $P(t, T)$ emitted by fission products, in the form of $\beta$ and $\gamma$-ray emissions, for a reactor operated for $T$ seconds at a constant power $P_0$, and for a cooling time of $t$ seconds after shutdown.

If $\beta(t)$ and $\gamma(t)$ represent the average energy emitted by $\beta$ and $\gamma$-ray emissions, for $t$ seconds after one fission event, one can use the following expressions:

$$ \beta(t) = 1.26 \times t^{-1.2} \text{ MeV/sec} \quad (\text{B.1}) $$

$$ \gamma(t) = 1.40 \times t^{-1.2} \text{ MeV/sec.} \quad (\text{B.2}) $$

According to reference [2], such expressions are valid from about $t = 1$ second to $t = 10^6$ seconds (or almost 12 days), and their integrations lead to the total power $P(t, T)$. Indeed, the total average energy emitted is the sum of equations (B.1) and (B.2), being approximate to $2.7 \times t^{-1.2}$; with the additional approximation that the recoverable energy per fission is 200 MeV, the expression of $P(t, T)$ for a reactor of $P_0$ Watts becomes

$$ P(t,T) = \frac{2.7 P_0}{200 \times 1.6 \times 10^{-13}} \int_t^{t+T} \tau^{-1.2} d\tau \quad \text{MeV/sec} \quad (\text{B.3}) $$

$$ = 4.1 \times 10^{11} P_0 [t^{-0.2} - (t+T)^{-0.2}] \quad \text{MeV/sec} \quad (\text{B.4}) $$

$$ = 0.065 P_0 [t^{-0.2} - (t+T)^{-0.2}] \quad \text{Watts.} \quad (\text{B.5}) $$

Equation (B.5) corresponds to the one mentioned in reference [2], and is widely used in reactor physics.

The origin of the semi-empirical equations (B.1) and (B.2) can be traced back to the 40s, as presented in reference [333]. The original experiments were performed prior to 1946 and presented to the American Physical Society in Chicago in June 1946 (recorded on page 115 of volume 70 of the Physical Review, July 1946), but the available documentation can be found in reference [38] from

1948. It is interesting to note that reference [333] presents a number of disagreements between the results from the Borst-Wheeler formula and other values from the Naval Radiological Defense Laboratory.

Although the type of fuel is not mentioned for the applicability of equation (8), it is based on a recoverable energy per fission of 200 MeV, close to the one for $^{235}$U (about 202 MeV).

## B.2. Untermeyer-Weills function, 1952

The report from Untermeyer and Weills [67] presents an empirical function for the decay heat of irradiated uranium, from 0.1 second up to $10^8$ seconds (about 3.2 years). It is based on the analysis of measured data, provided in a series of other reports, using a calorimeter as well as radiochemical measurements, combined with calculations for the decay of $^{239}$U and $^{239}$Np. The justification for propose a different function than equation (8) is that such an equation is based on a statistical analysis of measurements, which might not hold for short and long cooling times. Untermeyer and Weills propose the following “experimental data fit”, over the entire available range:

$$ P(t, T) = 0.1 P_0 \left[ (t+10)^{-0.2} - 0.87(t+2 \times 10^7)^{-0.2} - (t+T+10)^{-0.2} - 0.87(t+T+2 \times 10^7)^{-0.2} \right] \quad \text{Watts} \quad (\text{B.6}) $$

with, as previously all times in seconds, $t$ being the time since the operation (cooling time), $T$ the operation time, $P_0$ the reactor constant power. Reference. [67] also provides an estimate of accuracy (using the reference's vocabulary):

* $t < 1$ second: large error

* $1 < t < 10^2$ seconds: $\pm 50\%$

* $10^2 < t < 10^4$ seconds: $\pm 30\%$

* $10^4 < t < 10^6$ seconds: $\pm 10\%$

* $10^6 < t < 10^8$ seconds: $\pm 50\%$.

## B.3. Patterson-Schlitz function, 1955

The function for the decay heat proposed by Patterson and Schlitz in reference [334] is similar to the one of Untermeyer and Weills, presented in equation (B.6). No further information is available.

## B.4. Glasstone function, 1963

Two different expressions can be found in reference [129], with specific ranges of applicability. The solutions proposed in this reference also follow a pragmatic approach, given that the emission rates of beta particles and gamma-ray photons can be expressed by means of empirical expressions. It is stated in reference [129] that the applicability of such expressions is within a factor of two or less (Ref. [38] is used to justify such a statement).
For periods between 10 seconds up to several weeks after fission, the beta and gamma emission rates for one fission are presented as

$$ \beta(t) = 3.8 \times 10^{-6} t^{-1.2} \quad \text{particles/sec} \quad \text{(B.7)} $$

$$ = 1.5 \times 10^{-6} t^{-1.2} \quad \text{MeV/sec} \quad \text{(B.8)} $$

(considering that about 0.4 MeV are emitted by beta particles), and

$$ \gamma(t) = 1.9 \times 10^{-6} t^{-1.2} \quad \text{photons/sec} \quad \text{(B.9)} $$

$$ = 1.3 \times 10^{-6} t^{-1.2} \quad \text{MeV/sec} \quad \text{(B.10)} $$

(considering that the photon energy is about 0.7 MeV). The total energy released by beta particles and photons is given by the sum of (B.8) and (B.10). The total energy observed at a time *t* after shutdown, due to the emission during an interval *dT* is then given by $2.8 \times 10^{-6}t^{-1.2}$, given in MeV/sec These expressions are very close to equations (B.1) and (B.2) (although with the values for $\beta$ and $\gamma$ being exchanged), and leads to the same value of $P(t, T)$ from equation (8). It is mentioned that this expression can be applied for cooling time above 10 seconds, and it is assumed that it is applicable to reactor with a constant power.

Another expression is also given in reference [129], for shorter cooling time (from 1 sec upwards), and is obtained from the original work of Untermeyer and Weills [67].

## B.5. El-Wakil function, 1971

A different equation was proposed by M.M. El-Wakil for the decay heat, which is to be relevant for $t > 200$ seconds [335]:

$$ P(t, T) = P_0 \times 0.095t^{-0.26} \left[ 1 - \left( 1 + \frac{T}{t} \right)^{-0.2} \right] . \quad \text{(B.11)} $$

No additional information is available.

## B.6. Glasstone function, 1981

In the new edition of his book [129,132], Glasstone proposed a new approximation, based on the ANS standard (from 1971, or 1973, or 1979 ?):

$$ P(t, T) = 0.005 P_0 \times a \left[ t^{-b} - (t+T)^{-b} \right] \quad \text{Watts} \quad \text{(B.12)} $$

where *a* and *b* are two constants which depend on the time after shutdown *t* (see Table B.1).

It is indicated that such fit follows that ANS recommendation within ±6%.

## B.7. Todreas-Kazimi function, 1990

The expressions provided in reference [130] by Todreas and Kazimi basically refer to Glasstone (1963) for the

<table><caption>**Table B.1.** Parameter values for *a* and *b* as a function of the time after shutdown *t* for equation (B.12).</caption><thead><tr><td>Time after shutdown <i>t</i> (seconds)</td><td><i>a</i></td><td><i>b</i></td></tr></thead><tbody><tr><td>0.1 &lt; <i>t</i> &lt; 10</td><td>12.05</td><td>0.0639</td></tr><tr><td>10 &lt; <i>t</i> &lt; 150</td><td>15.31</td><td>0.1807</td></tr><tr><td>150 &lt; <i>t</i> &lt; 8 &times; 10<sup>8</sup></td><td>27.43</td><td>0.2962</td></tr></tbody></table>

decay heat from fission products and actinides decays. But they also present a practical expression for the decay heat due to delayed neutrons. In this case, solving neutron kinetic equations with a large negative insertion of reactivity leads to such decay heat value, expressed as $P_{\text{delayed neutrons}}(t, T)$:

$$ \begin{aligned} P_{\text{delayed neutrons}}(t, T) &= P_{\text{delayed neutrons}}(t) \\ &= P_0 \times (0.0625 \exp^{-0.0124t} \\ &\quad + 0.9375 \exp^{-960t}) . \end{aligned} \quad \text{(B.13)} $$

This expression is valid using typical quantities for $^{235}$U fuel, in the case of a water-moderated reactor. Due to the exponential values, the second term is negligible after 0.01 seconds, and the decay heat due to delayed neutrons exponentially decays over a period of about 80 seconds. In the case of the decay heat from fission products and actinides, reference [130] proposes the same expressions for $\beta(t)$ and $\gamma(t)$ as Glasstone [129], leading to equation (8).

## Appendix C A numerical simulation on Bremsstrahlung Ray Measurement for β-decay in Spent Nuclear Fuel: a Japanese example

This work on the simulation of Bremsstrahlung measurements was presented in detail in reference [336]. In Japan, dry storage of a spent nuclear fuel assembly (FA) is not allowed during a pre-determined time length after it is discharged from the core. Mainly due to the limitation of surface temperature of the fuel cladding, the time length is 15 years for an FA of which the maximum discharge burn-up is 48 GWd/t, for example. The time length is too conservative for a spent FA with lower burn-up. Then we proposed to allow storing FA in a dry cask provided the residual heat decays below a limitation value. The residual heat is estimated based on the reactor record and a burn-up calculation. Then we need validation of the residual heat estimation. However, it costs so much to equip a calorimeter to SFPs of every NPP. Then we consider applying gamma-ray spectroscopy to assess the decay heat. Currently, we can quantify the radioactivities of $^{137}$Cs→$^{137m}$Ba, $^{134}$Cs, $^{154}$Eu, $^{106}$Ru→$^{106}$Rh, and $^{144}$Ce→$^{144}$Pr. Their relative contribution to the total decay heat is around 40%.

With the ORIGEN-ARP module in SCALE-6.1, we have calculated time-dependent decay heat for FA discharged after from one to four cycles of irradiation in
<figure><img src="image_42.png" /><figcaption>Fig. C.1. Relative contribution of $^{90}$Sr→$^{90}$Y to total decay heat.</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>Fig. C.2. Calculated energy deposition in HPGe detector.</figcaption></figure>

PWR. The burn-up increment per cycle was assigned 13 GWd/t. The relative contribution of $^{90}$Sr→$^{90}$Y to the total decay heat is around 40%, see Figure C.1.

Then we picked up 27 nuclides of which the relative contribution to the total radio activities exceeds 0.01%. Then we prepare a calculation model where a single fuel rod with 2 cm in length and a cylindrical HPGe detector with 8 cm in length and 8 cm in diameter. The distance between the center of the cylindrical rod and the HPGe detector was 20 cm. Then, we assigned the radioactivities of the 27 nuclides homogeneously in the fuel pellets and calculated photon and electron transports with the PHITS

<figure><img src="image_44.png" /><figcaption>Fig. C.3. Relative contribution of Bremsstrahlung ray from beta decay of $^{90}$Y to the total count rate in the pulse height region from 1.6 to 1.87 MeV.</figcaption></figure>

code. The source spectra of gamma and beta rays listed in ICRP-107 were used for the calculation. The calculated energy deposition spectra in the HP-Ge detector 13 years after discharge from the core are shown in Figure C.2.

In addition to peak gamma rays from $^{134,137}$Cs, $^{154}$Eu, $^{144}$Pr, and $^{106}$Rh, the count rate of the Bremsstrahlung ray is significant, especially in the pulse height range from 1.5 to 1.8 MeV. The relative contribution of Bremsstrahlung ray from $^{90}$Y to the total count rate in pulse height region from 1.6 to 1.7 MeV is plotted in Figure C.3.

For the cooling times of 13 years, the contribution is greater than 94% and that is almost independent of the burn-up. Also, for FAs irradiated more than or equal to two cycles, the summation of the contributions of $^{90}$Sr→$^{90}$Y, $^{106}$Ru→$^{106}$Rh, and $^{144}$Ce→$^{144}$Pr reaches 98%. The radioactivities of $^{90}$Sr→$^{90}$Y can be quantified by the count rate of the pulse height region from 1.6 to 1.7 MeV and the peak count rates of $^{106}$Rh and $^{144}$Pr gamma rays. Summation of the relative contribution of $^{90}$Sr→$^{90}$Y, $^{137}$Cs→$^{137m}$Ba, $^{134}$Cs, $^{154}$Eu, $^{106}$Ru→$^{106}$Rh, and $^{144}$Ce→$^{144}$Pr to the total decay heat is 80%. Therefore, focusing on the Bremsstrahlung rays from the beta decay of $^{90}$Y, the decay heat estimation by photon spectrum measurement might become more credible.
## Appendix D Tables for TAGS data

<table><caption>**Table D.1.** List of nuclei updated in JEFF-3.1.1 with the Greenwood's TAGS mean decay energy for light particle component (ELP). The values from the different versions of JEFF and ENDF/B libraries are also included for comparison.</caption><thead><tr><th rowspan="2">Isotope</th><th colspan="6">ELP in keV</th></tr><tr><th>JEFF-3.1</th><th>JEFF-3.1.1</th><th>JEFF-3.3</th><th>ENDF/B-VII.1</th><th>ENDF/B-VIII.0</th><th>TAGS (Greenwood)</th></tr></thead><tbody><tr><td><sup>90</sup>Rb</td><td>2049</td><td>1916</td><td>1904</td><td>1905</td><td>1905</td><td>1916</td></tr><tr><td><sup>90</sup><sup>m</sup>Rb</td><td>1403</td><td>1118</td><td>1118</td><td>1115</td><td>1115</td><td>1118</td></tr><tr><td><sup>91</sup>Rb</td><td>1612</td><td>1368</td><td>1368</td><td>1370</td><td>1370</td><td>1368</td></tr><tr><td><sup>95</sup>Sr</td><td>2208</td><td>1901</td><td>1901</td><td>1892</td><td>1892</td><td>1901</td></tr><tr><td><sup>139</sup>Cs</td><td>1640</td><td>1671</td><td>1671</td><td>1664</td><td>1664</td><td>1671</td></tr><tr><td><sup>141</sup>Cs</td><td>1935</td><td>1506</td><td>1506</td><td>1490</td><td>1513</td><td>1506</td></tr><tr><td><sup>143</sup>Ba</td><td>1417</td><td>1195</td><td>1195</td><td>1182</td><td>1199</td><td>1195</td></tr><tr><td><sup>144</sup>Ba</td><td>1040</td><td>930</td><td>930</td><td>902</td><td>902</td><td>930</td></tr><tr><td><sup>145</sup>Ba</td><td>1860</td><td>1285</td><td>1285</td><td>1458</td><td>1458</td><td>1285</td></tr><tr><td><sup>142</sup>La</td><td>868</td><td>962</td><td>962</td><td>954</td><td>954</td><td>962</td></tr><tr><td><sup>143</sup>La</td><td>1237</td><td>1235</td><td>1235</td><td>1225</td><td>1232</td><td>1235</td></tr><tr><td><sup>144</sup>La</td><td>1382</td><td>986</td><td>986</td><td>1021</td><td>1031</td><td>986</td></tr><tr><td><sup>145</sup>La</td><td>1499</td><td>762</td><td>762</td><td>827</td><td>854</td><td>762</td></tr><tr><td><sup>148</sup>Ce</td><td>713</td><td>586</td><td>586</td><td>618</td><td>618</td><td>586</td></tr><tr><td><sup>147</sup>Pr</td><td>899</td><td>669</td><td>669</td><td>673</td><td>673</td><td>669</td></tr><tr><td><sup>148</sup>Pr</td><td>1679</td><td>1348</td><td>1348</td><td>1280</td><td>1280</td><td>1348</td></tr><tr><td><sup>148</sup><sup>m</sup>Pr</td><td>1701</td><td>1239</td><td>1239</td><td>1059</td><td>1059</td><td>1239</td></tr><tr><td><sup>149</sup>Pr</td><td>1286</td><td>811</td><td>811</td><td>780</td><td>780</td><td>811</td></tr><tr><td><sup>151</sup>Pr</td><td>1394</td><td>1114</td><td>1114</td><td>1137</td><td>1137</td><td>1114</td></tr><tr><td><sup>153</sup>Nd</td><td>1112</td><td>1153</td><td>1153</td><td>1176</td><td>1176</td><td>1153</td></tr><tr><td><sup>154</sup>Nd</td><td>937</td><td>856</td><td>856</td><td>837</td><td>837</td><td>856</td></tr><tr><td><sup>155</sup>Nd</td><td>1500</td><td>1085</td><td>1085</td><td>1215</td><td>1215</td><td>1085</td></tr><tr><td><sup>152</sup>Pm</td><td>1326</td><td>1304</td><td>1304</td><td>1305</td><td>1305</td><td>1304</td></tr><tr><td><sup>153</sup>Pm</td><td>686</td><td>663</td><td>663</td><td>661</td><td>661</td><td>663</td></tr><tr><td><sup>154</sup><sup>m</sup>Pm</td><td>1321</td><td>849</td><td>849</td><td>810</td><td>810</td><td>849</td></tr><tr><td><sup>156</sup>Pm</td><td>1717</td><td>1206</td><td>1206</td><td>1194</td><td>1194</td><td>1206</td></tr><tr><td><sup>157</sup>Pm</td><td>1453</td><td>1540</td><td>1540</td><td>1456</td><td>1456</td><td>1540</td></tr><tr><td><sup>157</sup>Sm</td><td>912</td><td>839</td><td>839</td><td>839</td><td>839</td><td>839</td></tr><tr><td><sup>158</sup>Sm</td><td>667</td><td>512</td><td>512</td><td>514</td><td>514</td><td>512</td></tr></tbody></table>
<table><caption>**Table D.2.** List of nuclei updated in JEFF 3.1.1 with the Greenwood's TAGS mean decay energy for the electromagnetic component (EEM). The values from the different versions of JEFF and ENDF/B libraries are also included for comparison.</caption><thead><tr><th rowspan="2">Isotope</th><th colspan="6">EEM in keV</th></tr><tr><th>JEFF-3.1</th><th>JEFF-3.1.1</th><th>JEFF-3.3</th><th>ENDF/B-VII.1</th><th>ENDF/B-VIII.0</th><th>TAGS (Greenwood)</th></tr></thead><tbody><tr><td><sup>90</sup>Rb</td><td>1982</td><td>2271</td><td>2280</td><td>2272</td><td>2272</td><td>2270</td></tr><tr><td><sup>90m</sup>Rb</td><td>3241</td><td>3869</td><td>3869</td><td>3866</td><td>3866</td><td>3869</td></tr><tr><td><sup>91</sup>Rb</td><td>2269</td><td>2706</td><td>2706</td><td>2708</td><td>2708</td><td>2706</td></tr><tr><td><sup>95</sup>Sr</td><td>1118</td><td>1790</td><td>1790</td><td>1790</td><td>1790</td><td>1790</td></tr><tr><td><sup>139</sup>Cs</td><td>345</td><td>305</td><td>305</td><td>303.</td><td>303</td><td>305</td></tr><tr><td><sup>141</sup>Cs</td><td>770</td><td>1701</td><td>1701</td><td>1708</td><td>1708</td><td>1701</td></tr><tr><td><sup>143</sup>Ba</td><td>1417</td><td>1341</td><td>1341</td><td>1343</td><td>1343</td><td>1341</td></tr><tr><td><sup>144</sup>Ba</td><td>1040</td><td>785</td><td>785</td><td>785</td><td>785</td><td>785</td></tr><tr><td><sup>145</sup>Ba</td><td>1860</td><td>1833</td><td>1833</td><td>1831</td><td>1831</td><td>1833</td></tr><tr><td><sup>142</sup>La</td><td>2325</td><td>2121</td><td>2121</td><td>2117</td><td>2117</td><td>2121</td></tr><tr><td><sup>143</sup>La</td><td>252</td><td>439</td><td>439</td><td>424</td><td>438</td><td>439</td></tr><tr><td><sup>144</sup>La</td><td>2330</td><td>3085</td><td>3085</td><td>3158</td><td>3158</td><td>3085</td></tr><tr><td><sup>145</sup>La</td><td>624</td><td>2144</td><td>2144</td><td>2144</td><td>2144</td><td>2144</td></tr><tr><td><sup>148</sup>Ce</td><td>713</td><td>486</td><td>486</td><td>486</td><td>486</td><td>486</td></tr><tr><td><sup>147</sup>Pr</td><td>899</td><td>929</td><td>929</td><td>929</td><td>929</td><td>929</td></tr><tr><td><sup>148</sup>Pr</td><td>938</td><td>1776</td><td>1776</td><td>1777</td><td>1777</td><td>1776</td></tr><tr><td><sup>148m</sup>Pr</td><td>937</td><td>2050</td><td>2050</td><td>2332</td><td>2332</td><td>2050</td></tr><tr><td><sup>149</sup>Pr</td><td>305</td><td>1332</td><td>1332</td><td>1332</td><td>1332</td><td>1332</td></tr><tr><td><sup>151</sup>Pr</td><td>1394</td><td>1363</td><td>1363</td><td>1363</td><td>1363</td><td>1363</td></tr><tr><td><sup>153</sup>Nd</td><td>1112</td><td>509</td><td>509</td><td>443</td><td>443</td><td>509</td></tr><tr><td><sup>154</sup>Nd</td><td>937</td><td>551</td><td>551</td><td>546</td><td>546</td><td>551</td></tr><tr><td><sup>152</sup>Pm</td><td>285</td><td>345</td><td>345</td><td>318</td><td>318</td><td>345</td></tr><tr><td><sup>153</sup>Pm</td><td>73</td><td>123</td><td>123</td><td>121</td><td>121</td><td>123</td></tr><tr><td><sup>154m</sup>Pm</td><td>1321</td><td>1876</td><td>1876</td><td>1877</td><td>1877</td><td>1876</td></tr><tr><td><sup>156</sup>Pm</td><td>1717</td><td>2206</td><td>2206</td><td>2205</td><td>2205</td><td>2206</td></tr><tr><td><sup>157</sup>Pm</td><td>1453</td><td>849</td><td>849</td><td>849</td><td>849</td><td>849</td></tr><tr><td><sup>157</sup>Sm</td><td>912</td><td>585</td><td>585</td><td>585</td><td>585</td><td>585</td></tr><tr><td><sup>158</sup>Sm</td><td>667</td><td>590</td><td>590</td><td>594</td><td>594</td><td>590</td></tr></tbody></table>
<table><caption>**Table D.3.** List of nuclei measured per the Total Absorption Gamma Spectroscopy method per the TAS and MTAS collaborations with the associated mean decay energy for light particle (ELP). The values from the different versions of JEFF, JENDL and ENDF/B libraries are also included for comparison. The TAGS values already included in the different versions of libraries are in grey.</caption><thead><tr><th rowspan="2">Isotope</th><th colspan="7">ELP in keV</th><th rowspan="2">Reference</th></tr><tr><th>JEFF-3.1.1</th><th>JEFF-3.3</th><th>JENDL-5.0</th><th>ENDF/B-VII.1</th><th>ENDF/B-VIII.0</th><th>TAGS</th></tr></thead><tbody><tr><td><sup>86</sup>Br</td><td>1943</td><td>1943</td><td>1687</td><td>1944</td><td>1944</td><td>1687-1730</td><td>[202,208]</td></tr><tr><td><sup>87</sup>Br</td><td>1577</td><td>1170</td><td>1170</td><td>1657</td><td>1657</td><td>1170</td><td>[203]</td></tr><tr><td><sup>88</sup>Br</td><td>2394</td><td>1706</td><td>1706</td><td>2393</td><td>1702</td><td>1706</td><td>[203]</td></tr><tr><td><sup>89</sup>Kr</td><td>1377</td><td>1468</td><td>1250</td><td>1468</td><td>1468</td><td>1250</td><td>[208]</td></tr><tr><td><sup>90</sup>Kr</td><td>1459</td><td>1332</td><td>1139</td><td>1368</td><td>1368</td><td>1130</td><td>[208]</td></tr><tr><td><sup>89</sup>Rb</td><td>929</td><td>929</td><td>935</td><td>969</td><td>969</td><td>930</td><td>[208]</td></tr><tr><td><sup>90</sup>Rb</td><td>1916</td><td>1904</td><td>1916</td><td>1905</td><td>1905</td><td>1920</td><td>[208]</td></tr><tr><td><sup>90m</sup>Rb</td><td>1118</td><td>1118</td><td>1153</td><td>1115</td><td>1115</td><td>1100</td><td>[208]</td></tr><tr><td><sup>91</sup>Rb</td><td>1368</td><td>1368</td><td>1821</td><td>1370</td><td>1370</td><td>1389</td><td>[202]</td></tr><tr><td><sup>92</sup>Rb</td><td>2875</td><td>3496</td><td>3498</td><td>2867</td><td>2867</td><td>3496-3570</td><td>[201,208]</td></tr><tr><td><sup>94</sup>Rb</td><td>3102</td><td>2450</td><td>2452</td><td>2020</td><td>2020</td><td>2450</td><td>[203]</td></tr><tr><td><sup>95</sup>Rb</td><td>2829</td><td>2829</td><td>2580</td><td>2300</td><td>1380</td><td>2573</td><td>[205]</td></tr><tr><td><sup>96</sup>Y</td><td>3205</td><td>3203</td><td>2656</td><td>3206</td><td>3206</td><td>3200-3209</td><td>[206,207]</td></tr><tr><td><sup>96m</sup>Y</td><td>1851</td><td>1851</td><td>1610</td><td>1851</td><td>1630</td><td>1739</td><td>[206]</td></tr><tr><td><sup>100</sup>Nb</td><td>2493</td><td>2554</td><td>2422</td><td>2552</td><td>2252</td><td>2414</td><td>[204]</td></tr><tr><td><sup>100m</sup>Nb</td><td>2047</td><td>2047</td><td>1715</td><td>2010</td><td>2010</td><td>1706</td><td>[204]</td></tr><tr><td><sup>101</sup>Nb</td><td>1863</td><td>1962</td><td>2045</td><td>1960</td><td>1960</td><td>1797</td><td>[199]</td></tr><tr><td><sup>102</sup>Nb</td><td>2402</td><td>2402</td><td>1948</td><td>2306</td><td>2306</td><td>1948</td><td>[204]</td></tr><tr><td><sup>102m</sup>Nb</td><td>2276</td><td>2276</td><td>2829</td><td>2420</td><td>2420</td><td>2829</td><td>[204]</td></tr><tr><td><sup>105</sup>Mo</td><td>1922</td><td>1049</td><td>1087</td><td>1049</td><td>1049</td><td>1049</td><td>[199]</td></tr><tr><td><sup>102</sup>Tc</td><td>1945</td><td>1945</td><td>1935</td><td>1945</td><td>1945</td><td>1935</td><td>[199,200]</td></tr><tr><td><sup>104</sup>Tc</td><td>1595</td><td>931</td><td>941</td><td>931</td><td>931</td><td>931</td><td>[199,200]</td></tr><tr><td><sup>105</sup>Tc</td><td>1310</td><td>764</td><td>764</td><td>764</td><td>764</td><td>764</td><td>[199,200]</td></tr><tr><td><sup>106</sup>Tc</td><td>1943</td><td>1457</td><td>1467</td><td>1457</td><td>1457</td><td>1457</td><td>[199]</td></tr><tr><td><sup>107</sup>Tc</td><td>2056</td><td>1263</td><td>1311</td><td>1263</td><td>1263</td><td>1263</td><td>[199]</td></tr><tr><td><sup>137</sup>I</td><td>1862</td><td>1862</td><td>1934</td><td>1920</td><td>1920</td><td>1934-1870</td><td>[205,209]</td></tr><tr><td><sup>137</sup>Xe</td><td>1695</td><td>1705</td><td>1706</td><td>1700</td><td>1700</td><td>510</td><td>[209]</td></tr><tr><td><sup>139</sup>Xe</td><td>1686</td><td>1720</td><td>1616</td><td>1808</td><td>1808</td><td>1580</td><td>[208]</td></tr><tr><td><sup>142</sup>Cs</td><td>2899</td><td>2918</td><td>2482</td><td>2918</td><td>2456</td><td>2489</td><td>[208]</td></tr></tbody></table>
<table><caption>**Table D.4.** List of nuclei measured per the Total Absorption Gamma Spectroscopy method per the TAS and MTAS collaborations with the associated mean electromagnetic decay energy (EEM). The values from the different versions of JEFF, JENDL and ENDF/B libraries are also included for comparison. The TAGS values already included in the different versions of libraries are in grey.</caption><thead><tr><th rowspan="2">Isotope</th><th colspan="7">EEM in keV</th><th rowspan="2">Reference</th></tr><tr><th>JEFF-3.1.1</th><th>JEFF-3.3</th><th>JENDL-5.0</th><th>ENDF/B-VII.1</th><th>ENDF/B-VIII.0</th><th>TAGS</th></tr></thead><tbody><tr><td><sup>86</sup>Br</td><td>3297</td><td>3297</td><td>3782</td><td>3297</td><td>3297</td><td>3782-3720</td><td>[202,208]</td></tr><tr><td><sup>87</sup>Br</td><td>3089</td><td>3938</td><td>3938</td><td>3345</td><td>3345</td><td>3938</td><td>[203]</td></tr><tr><td><sup>88</sup>Br</td><td>3112</td><td>4609</td><td>4609</td><td>3134</td><td>3134</td><td>4609</td><td>[203]</td></tr><tr><td><sup>89</sup>Kr</td><td>1836</td><td>1836</td><td>2240 2240</td><td>1931</td><td>1931</td><td>2240</td><td>[208]</td></tr><tr><td><sup>90</sup>Kr</td><td>1459</td><td>1282</td><td>1690</td><td>1320</td><td>1320</td><td>1690</td><td>[208]</td></tr><tr><td><sup>89</sup>Rb</td><td>2234</td><td>2234</td><td>2228</td><td>2236</td><td>2236</td><td>2260</td><td>[208]</td></tr><tr><td><sup>90</sup>Rb</td><td>2271</td><td>2280</td><td>2227</td><td>2272</td><td>2272</td><td>2300-3240</td><td>[208]</td></tr><tr><td><sup>90m</sup>Rb</td><td>3869</td><td>3869</td><td>3240</td><td>3866</td><td>3866</td><td>4000</td><td>[208]</td></tr><tr><td><sup>91</sup>Rb</td><td>2706</td><td>2706</td><td>2672</td><td>2708</td><td>2708</td><td>2669</td><td>[202]</td></tr><tr><td><sup>92</sup>Rb</td><td>1750</td><td>464</td><td>461</td><td>2147</td><td>2147</td><td>464-385</td><td>[201,208]</td></tr><tr><td><sup>94</sup>Rb</td><td>2747</td><td>4063</td><td>4063</td><td>1895</td><td>1895</td><td>4063</td><td>[203]</td></tr><tr><td><sup>95</sup>Rb</td><td>2642</td><td>2642</td><td>3110</td><td>2162</td><td>2162</td><td>3110</td><td>[205]</td></tr><tr><td><sup>96</sup>Y</td><td>80</td><td>80</td><td>1205</td><td>80</td><td>80</td><td>100 67</td><td>[207] [206]</td></tr><tr><td><sup>96m</sup>Y</td><td>4486</td><td>4486</td><td>4307</td><td>4308</td><td>4308</td><td>4467</td><td>[206]</td></tr><tr><td><sup>100</sup>Nb</td><td>708</td><td>708</td><td>959</td><td>708</td><td>708</td><td>959</td><td>[204]</td></tr><tr><td><sup>100m</sup>Nb</td><td>2064</td><td>2064</td><td>2209</td><td>2213</td><td>2213</td><td>2763</td><td>[204]</td></tr><tr><td><sup>101</sup>Nb</td><td>245</td><td>270</td><td>471</td><td>270</td><td>270</td><td>445</td><td>[199]</td></tr><tr><td><sup>102</sup>Nb</td><td>2402</td><td>2402</td><td>2764</td><td>2094</td><td>2094</td><td>2764</td><td>[204]</td></tr><tr><td><sup>102m</sup>Nb</td><td>2094</td><td>2094</td><td>1023</td><td>2420</td><td>2420</td><td>1023</td><td>[204]</td></tr><tr><td><sup>105</sup>Mo</td><td>551</td><td>2407</td><td>2410</td><td>2407</td><td>2407</td><td>2407</td><td>[199]</td></tr><tr><td><sup>102</sup>Tc</td><td>81</td><td>81</td><td>106</td><td>81</td><td>81</td><td>106</td><td>[199] [200]</td></tr><tr><td><sup>104</sup>Tc</td><td>1890</td><td>3229</td><td>3229</td><td>3229</td><td>3229</td><td>3229</td><td>[199] [200]</td></tr><tr><td><sup>105</sup>Tc</td><td>668</td><td>1825</td><td>1825</td><td>1825</td><td>1825</td><td>1825</td><td>[199] [200]</td></tr><tr><td><sup>106</sup>Tc</td><td>2191</td><td>3132</td><td>3132</td><td>3132</td><td>3132</td><td>3132</td><td>[199]</td></tr><tr><td><sup>107</sup>Tc</td><td>515</td><td>1822</td><td>1822</td><td>1822</td><td>1822</td><td>1822</td><td>[199]</td></tr><tr><td><sup>137</sup>I</td><td>1220</td><td>1220</td><td>1220</td><td>1135</td><td>1135</td><td>1220-1250</td><td>[205,209]</td></tr><tr><td><sup>137</sup>Xe</td><td>190</td><td>189</td><td>189</td><td>190</td><td>190</td><td>190</td><td>[209]</td></tr><tr><td><sup>139</sup>Xe</td><td>1686</td><td>1017</td><td>1341</td><td>1027</td><td>1027</td><td>1340</td><td>[208]</td></tr><tr><td><sup>142</sup>Cs</td><td>675</td><td>676</td><td>1720</td><td>952</td><td>1704</td><td>1725</td><td>[208]</td></tr></tbody></table>
## References

1. A. Bevilacqua, J.P. Carreton, A. Guskov, T. Hornebrant, C. Jones, Y. Kumano, et al. Storage of Spent Nuclear Fuel, specific safety guide. (International Atomic Energy Agency, IAEA, Vienna, Austria, 2020), SSG-15, Rev. 1. https://www-pub.iaea.org/MTCD/publications/PDF/P1882_web.pdf

2. J.R. Lamarch, Introduction to Nuclear Reactor Theory (Addison-Wesley Publishing Company, 1966)

3. C. Poinssot, P. Toulhoat, J.P. Grouiller, J. Pavageau, J.P. Piron, M. Pelletier, et al. Synthesis on the Long Term Behavior of Spent Nuclear Fuel – Volume 1 and 2 (Commissariat à l'Energie Atomique, CEA, 2001). https://inis.iaea.org/collection/NCLCollectionStore/_Public/33/009/33009840.pdf

4. B. Duchemin, C. Nordborg, Decay heat calculation, An international Nuclear Code Comparison. Nuclear Energy Agency, NEACRP-319 and NEANDC-275; 1990. https://inis.iaea.org/collection/NCLCollectionStore/_Public/22/009/22009691.pdf

5. NEA Committees on reactor physics and nuclear data. Proceedings of a specialists' meeting on Data for Decay Heat Predictions (Nuclear Energy Agency, NEACRP-302 and NEANDC-245, 1987). https://www.oecd-nea.org/upload/docs/application/pdf/2020-01/neacrp-1-1987-302.pdf

6. C. Nordborg, NEACRP/NEANDC Task Force Meeting on Decay Heat Predictions (Nuclear Energy Agency, EACRP-A-1039 and NEANDC-A-26, 1989). https://www.oecd-nea.org/upload/docs/application/pdf/2020-07/neacrp-a-1989-1039.pdf.

7. M.A. Kellett, A.L. Nichols, O. Bersillon, H. Henriksson, R. Jacqmin, B. Roque, et al. Assessment of Fission Product Decay Data for Decay Heat Calculations, A report by the Working Party on International Evaluation Co-operation of the NEA Nuclear Science Committee, 2007. ISBN 978-92-64-99034-0.

8. Decay heat of spent nuclear fuels, Accessed: 2022-09-01, https://www.oecd-nea.org/download/wpncs/sg12/

9. K. Asfaw, K. Kamimura, T. Nakajima, S. Schultz, P. Shaw, K. Sim, et al. Design of the Reactor Core for Nuclear Power Plants (International Atomic Energy Agency, IAEA, Vienna, Austria, 2019). https://www-pub.iaea.org/MTCD/publications/PDF/PUB1859web.pdf

10. S.J. Baik, J. Beard, E. Courtin, N. Fil, M. Gasparini, C. Jackson, et al. Design of the Reactor Coolant System and Associated Systems for Nuclear Power Plants (International Atomic Energy Agency, IAEA, Vienna, Austria, 2020). https://www-pub.iaea.org/MTCD/publications/PDF/PUB1878_web.pdf

11. G. Ilas, J. Burns, LWR Decay Heat Analysis with SCALE 6.3 and ENDF/B-VIII.0 Nuclear Data Libraries (International Conference TopFuel 2021, Santander, Spain, October 2021), pp. 24-28

12. H. Akkurt, H. Liljenfeldt, G. Ilas, S. Baker, Phenomena Identification and Ranking Table (PIRT) for Decay Heat (Electric Power Research Institute, EPRI, USA, 2020), 3002018440

13. T. Yamamoto, D. Iwahashi, Validation of decay heat calculation results of ORIGEN2.2 and CASMO5 for light water reactor fuel, J. Nucl. Sci. Technol. 53, 2108 (2016). https://doi.org/10.1080/00223131.2016.1183528

14. A. Shama, D. Rochman, S. Caruso, A. Pautz, Validation of spent nuclear fuel decay heat calculations using Polaris, ORIGEN and CASMO5, Ann. Nucl. Energy. 165, 108758 (2022). https://www.sciencedirect.com/science/article/pii/S0306454921006344

15. P. Jansson, M. Bengtsson, U. Bäckström, F. Àlvarez-Velarde, D. Calic, S. Caruso, et al. Blind Benchmark Exercise for Spent Nuclear Fuel Decay Heat, Nucl. Sci. Eng. 196, 1125 (2022)

16. F. Laugier, C. Diop, S. Ebalard, C. Garzenne, A. Sargeni, A new decay heat standard proposition based on a technical specifications guide for computation codes. In: Proc.

Int. Conf. PHYSOR 2008, Interlaken, Switzerland, 14-19 September (2008)

17. J.C. Wagner, Spent Nuclear Fuel Disposition (Idaho National Laboratory, USA, 2016)

18. IAEA, Management of Spent Fuel from Nuclear Power Reactors Learning from the Past, Enabling the Future (International Atomic Energy Agency, Vienna, Austria, 2019). http://www.iaea.org/bulletin

19. NWMO, What is used Nuclear Fuel? (Nuclear Waste Management Organization, Ontario, Canada, 2022). https://www.nwmo.ca

20. ASTM american standards, Standard Guide For Characterization Of Spent Nuclear Fuel In Support Of Interim Storage, Transportation And Geologic Repository Disposal (American Nuclear Society; 2021), ASTM C1682-21. https://webstore.ansi.org/Standards/ASTM/astmc168221

21. IAEA. Joint Convention on the Safety of Spent Fuel Management and on the Safety of Radioactive Waste Management (IAEA International Law Series, 2006)

22. DIN Standards Committee Materials Testing, Calculation of the decay power in nuclear fuels of light water reactors - Part 1: Uranium oxide nuclear fuel for pressurized water reactors, English translation of DIN 25463-1:2014-02 (DIN Standards Committee Materials Testing, Germany, 2014)

23. FSC-6910, Decay Heat. In: DOE Fundamentals Handbook, No. DOE-HDBK-1012/2-92 in 1 (U.S., Department of Energy, 1992)

24. Technical Committee ISO/TC 85, Nuclear energy, Subcommittee SC 3, Power Reactor Technology, Nuclear Energy – Light Water Reactors – Calculation of the Decay Heat Power in Nuclear Fuels (ISO International Standard, 1992)

25. DIN Standards Committee Materials Testing, Calculation of the Decay Power in Nuclear Fuels of Light Water Reactors – Part 2: Mixed-uranium-plutonium oxide (MOX) nuclear fuel for pressurized water reactors, English translation of DIN 25463-2:2014-02 (DIN Standards Committee Materials Testing, Germany, 2014)

26. American Nuclear Society Standards Committee Working Group ANS-5.1 Decay Heat Power in Light Water Reactors (American Nuclear Society, 2014), ANSI/ANS-5.1-2014

27. Technical Committee ISO/TC 85, Nuclear energy, Subcommittee SC 6, Power Reactor Technology, Nuclear Energy – Light Water reactors – Calculation of the Decay Heat Power in Nuclear Fuels (ISO International Standard, 2022), ISO 10645:2022

28. A. Sotomayor-Rivera, Spent Fuel Heat Generation in an Independent Spent Fuel Storage Installation (U.S. Nuclear Regulatory Commission, USA, 2018), US NRC Regulatory Guide 3.54, Rev.2

29. Fuel Burnup, Accessed: 2022-06-23, https://www.nuclear-power.com/nuclear-power/reactor-physics/reactor-operation/fuel-burnup/

30. M. Aissa, Fuel Burnup Plant Records: Generation and Accuracy, In: *Proceedings of International Workshop on Advances in Applications of Burnup Credit for Spent Fuel Storage, Transport, Reprocessing, and Disposition, 27-30 October. Cordoba, Spain* (2009). https://www.oecd-nea.org/science/wpncs/buc/workshop-2009/

31. B.B. Bevard, J.C. Wagner, C.V. Parks, M. Aissa, Review of Information for Spent Nuclear Fuel Burnup Confirmation (Oak Ridge National Laboratory, USA, 2009), NUREG/CR-6998

32. F. Michel-Sendis, I. Gauld, J.S. Martinez, C. Alejano, M. Bossant, D. Boulanger, et al. SFCOMPO-2.0: An OECD NEA database of spent nuclear fuel isotopic assays, reactor design specifications, and operating data, Ann. Nucl. Energy 110, 779 (2017). https://www.sciencedirect.com/science/article/pii/S0306454917302104

33. D. Rochman, A. Vasiliev, H. Ferroukhi, M. Hursin, R. Ichou, J. Taforeau, et al. Analysis for the ARIANE GU3 sample: nuclide
inventory and decay heat, EPJ Nucl. Sci. Technol. **7**, 14 (2021).
https://doi.org/10.1051/epjn/2021013

34. D. Rochman, F. Álvarez-Velarde, R. Dagan, L. Fiorito, S. Häkkinen, M. Kromar, et al., On the estimation of nuclide inventory and decay heat: a review from the EURAD European project, EPJ Nucl. Sci. Technol. **9**, 14 (2023). https://doi.org/10.1051/epjn/2022055

35. K.S. Smith, S. Tarves, T. Bahadir, R. Ferrer, Benchmarks for Quantifying Fuel Reactivity Depletion Uncertainty (Electric Power Research Institute, EPRI, USA, 2011), p. 1022909

36. S. Börresen, A. Becker, Calculation of neutron emission from spent LWR fuel assemblies: SNF method and validation, In: *International Conference on Reactor Physics, Nuclear Power: A Sustainable Resource, September 14-19*. Interlaken, Switzerland (2008)

37. P. Schillebeeckx, M. Verwerft, P. Romojaro, G. Žerovnik, N. Messaoudi, G. Alaerts, et al. An absolute measurement of the neutron production rate of a spent nuclear fuel sample used for depletion code validation, Front. Energy Res. **11**, 1162367 (2023). https://www.frontiersin.org/articles/10.3389/fenrg.2023.

38. K. Way, E.P. Wigner, The Rate of Decay of Fission Products, Phys. Rev. **73**, 1318 (1948). https://link.aps.org/doi/10.1103/PhysRev.73.1318

39. Y. Liu, B. Kochunas, W. Martin, T. Downar, Delayed fission energy effect on LWR normal operation and transients, Ann. Nucl. Energy **128**, 84 (2019). https://www.sciencedirect.com/science/article/pii/S0306454918307126

40. O.W. Hermann, C.V. Parks, J.P. Renier, Technical Support for a Proposed Decay Heat Guide Using Sas2H/Origen-S Data (Oak Ridge National Laboratory, 1994), NREG/CR-5625, ORNL-6698. https://www.osti.gov/servlets/purl/43752

41. U.S. NRC, Spent Fuel Heat Generation in an Independent Spent Fuel Storage Installation, U.S. Nuclear Regulatory Commission, USA, US NRC Regulatory Guide 3.54, Rev.1., 1999. https://www.nrc.gov/docs/ML0037/ML003761667.pdf

42. I.C. Gauld, J.C. Ryman, Nuclide Importance to Criticality Safety, Decay Heating, and Source Terms Related to Transport and Interim Storage of High-Burnup LWR Fuel (Oak Ridge National Laboratory, USA, 2001)

43. G. Žerovnik, P. Schillebeeckx, K. Govers, A. Borella, D. Ćalić, L. Fiorito, et al. Observables of interest for the characterisation of Spent Nuclear Fuel (Publications Office of the European Union, Luxembourg, 2018), EUR 29301 EN, ISBN 978-92-79-90347-2

44. A. Tobias, Decay heat, Prog. Nucl. Energy **5**, 1 (1980). https://www.sciencedirect.com/science/article/pii/0149197080900025

45. F. Storrer, Review of Decay Heat Predictions and Standards (Meeting of the JEF Working Group on Radioactive Decay and Fission Yield Data, BNFL Office, London, UK., 1994), JEF/DOC-473

46. K. Shure, Decay Heat and Decay Rate of Actinides in Highly Neutron-Irradiated Uranium Initially of High <sup>235</sup>U Content, Nucl. Sci. Eng. **85**, 51 (1983). https://doi.org/10.13182/NSE83-A17151

47. T.R. England, W.B. Wilson, TMI-2 Decay Power: LASL Fission-Product and Actinide; Decay Power Calculations for the President's Commission on the Accident at Three Mile Island TMI-2 (Los Alamos Scientific Laboratory, NM, USA, 1979) LA-8041-MS. https://inis.iaea.org/collection/NCLCollectionStore/_Public/11/528/11528011.pdf

48. H. Smith, The Relative Significance of Actinide Decay Heat for Nuclear Fuel Characterization, In: *Transactions of the American Nuclear Society, Vol. 85, 2001 Winter Meeting, Reno, Nevada* (2001), p. 318

49. K. Tasaka, J. Katakura, T. Yoshida, T. Kato, R. Nakasima, Recommended values of decay heat power and method to utilize the data (Japan Atomic Energy Research Institute, 1991), JAERI-M 91-034

50. J.H. Mairs, S. Nair, The Inventories of Actinides and Fission Product Arising in Spent Nuclear Fuel: Results Using the Rice Code (Berkeley Nuclear Labs, USA, 1979), CEGB-RD-B-N-4579

51. D.H. Stoddard, E.L. Albenesius, Radiation Properties of <sup>238</sup>Pu Produced for Isotopic Power Generators (Savannah River Laboratory, USA, 1965), TID-4500. https://www.osti.gov/servlets/purl/4616364

52. I.C. Gauld, B.D. Murphy, Proposed Expansion of Regulatory Guide 3.54-Decay Heat Generation in an Independent Spent Fuel Storage Installation (Oak Ridge National Laboratory, 2010) ORNL/TM-2007/231

53. M.W. Francis, C.F. Weber, M.T. Pigni, I.C. Gauld, Reactor Fuel Isotopics and Code Validation for Nuclear Applications (Oak Ridge National Laboratory, USA, 2014), ORNL/TM-2014/464

54. K. Debertin, H. Ramthun, FPND needed for Nondestruc- tive Burnup Determination, In: *Proceedings of a panel on fission product nuclear data organized by the International Atomic Energy Agency, November 26-30*. Bologna, Italy (1974). https://www-nds.iaea.org/publications/tecdocs/iaea-0169-volume3.pdf

55. D. Rochman, J. Taforeau, T. Simeonov, A. Shama, Com- parison of calculated and measured spent nuclear fuel decay heat with CASMO5, SNF and standard methods, Nucl. Eng. Design **410**, 112392 (2023). https://www.sciencedirect.com/science/article/pii/S0029549323002418

56. B.L. Broadhead, M.D. Dehart, J.C. Ryman, J.S. Tang, C.V. Parks, Investigation of Nuclide Importance to Functional Requirements Related to Transport and Long-Term Storage of LWR Spent Fuel (Oak Ridge National Laboratory, 1995), ORNL/TM-12742

57. K. Tasaka, J. Katakura, T. Yoshida, Recommended values of decay heat power and method to utilize the data (Japan Atomic Energy Research Institute (JAERI), Japan, 1991), JAERI-M 91-034, NEANDC(J)-161/U and INDC(JPN)-149/L

58. A.J.M. Plompen, O. Cabellos, C. De Saint Jean, M. Fleming, A. Algora, M. Angelone, et al., The joint evaluated fission and fusion nuclear data library, JEFF-3.3, Eur. Phys. J. A **56**, 181 (2020). https://doi.org/10.1140/epja/s10050-020-00141-9

59. F.P. Roberts, The Radiation Characteristics of Reactor Pro- duced Rhodium, Palladium, Ruthenium and Technetium (Pacific Northwest Laboratory, Richland, Washington, USA, 1972), BNWL-1693

60. Subcommittee on Radiochemistry; Subcommittee on the Use of Radioactivity Standards; Committee on Nuclear Science; National Research Council, Users' Guide for Radioactivity Standards (The National Academies Press, Washington, DC, 1974). https://nap.nationalacademies.org/catalog/20179/users-guide-for-radioactivity-standards

61. S. Nakamura, H. Harada, T. Katoh, Measurement of Thermal Neutron Capture Cross Section and Resonance Integral of the Reaction <sup>133</sup>Cs(n,γ)<sup>134m,134g</sup>Cs, J. Nucl. Sci. Technol. **36**, 847 (1999). https://doi.org/10.1080/18811248.1999.9726275

62. B.F. Rider, J.L. Russell, Jr., D.W. Harris, J.P. Peterson Jr., The Determination of Uranium Burnup in MWd/ton (Vallecitos Atomic Laboratory General Electric Company Pleasanton, California, 1960), GEAP-3373. https://www.osti.gov/servlets/purl/4156626

63. D. Harris, J. Epstein, Properties of Selected Radioisotopes (NASA Goddard Space Flight Center Greenbelt, MD, US, 1968), NASA SP-7031. https://ntrs.nasa.gov/citations/19680020487

64. Practices and Developments in Spent Fuel Burnup Credit Applications. No. 1378 in TECDOC Series (CD-ROM). Vienna, Austria: International Atomic Energy Agency, 2003

65. S.E. Skutnik, Proposed re-evaluation of the 154Eu thermal (n,γ) capture cross-section based on spent fuel benchmarking studies, Ann. Nucl. Energy **99**, 80 (2017).