Nuclear Data Uncertainty Quantification in Criticality Safety Evaluations for Spent Nuclear Fuel Geological Disposal

Presently, a criticality safety evaluation methodology for the final geological disposal of Swiss spent nuclear fuel is under development at the Paul Scherrer Institute in collaboration with the Swiss National Technical Competence Centre in the field of deep geological disposal of radioactive waste. This method in essence pursues a best estimate plus uncertainty approach and includes burnup credit. Burnup credit is applied by means of a computational scheme called BUCSS-R (Burnup Credit System for the Swiss Reactors–Repository case) which is complemented by the quantification of uncertainties from various sources. BUCSS-R consists in depletion, decay and criticality calculations with CASMO5, SERPENT2 and MCNP6, respectively, determining the keff eigenvalues of the disposal canister loaded with the Swiss spent nuclear fuel assemblies. However, the depletion calculation in the first and the criticality calculation in the third step, in particular, are subject to uncertainties in the nuclear data input. In previous studies, the effects of these nuclear data-related uncertainties on obtained keff values, stemming from each of the two steps, have been quantified independently. Both contributions to the overall uncertainty in the calculated keff values have, therefore, been considered as fully correlated leading to an overly conservative estimation of total uncertainties. This study presents a consistent approach eliminating the need to assume and take into account unrealistically strong correlations in the keff results. The nuclear data uncertainty quantification for both depletion and criticality calculation is now performed at once using one and the same set of perturbation factors for uncertainty propagation through the corresponding calculation steps of the evaluation method. The present results reveal the overestimation of nuclear data-related uncertainties by the previous approach, in particular for spent nuclear fuel with a high burn-up, and underline the importance of consistent nuclear data uncertainty quantification methods. However, only canister loadings with UO2 fuel assemblies are considered, not offering insights into potentially different trends in nuclear data-related uncertainties for mixed oxide fuel assemblies.


Introduction
By 2024, a general license application for a deep geological repository for Swiss spent nuclear fuel (SNF) is planned for submission by the Swiss National Technical Competence Centre in the field of deep geological disposal of radioactive waste (NAGRA) [1]. In collaboration with NAGRA, the Paul Scherrer Institute (PSI) currently develops a criticality safety evaluation (CSE) methodology for an assessment of NAGRA's disposal concept addressing the criticality safety requirements for such an installation. The methodology is based on a best-estimate-plus-uncertainty (BEPU) approach with certain bounding and conservative assumptions and conditions for fuel assembly (FA) burnup profiles.
Previous scoping studies performed by NAGRA concluded that the nuclear criticality safety cannot be guaranteed without taking credit for the SNF reactivity change with burnup [2]. For this reason, the CSE method is taking burnup credit (BUC) into account.
The criticality evaluations at PSI are performed by means of a computational scheme called BUCSS-R (Burnup Credit System for the Swiss Reactors-Repository case) [3] and complemented by a comprehensive quantification of uncertainties originating from a variety of sources. The information of which spent (or used) nuclear fuel assemblies can be loaded into the final disposal canister, fulfilling the regulatory safety limit from the criticality point of view, is finally provided in the form of loading curves.
The objective of the present study is to enhance the existing uncertainty quantification (UQ) method with respect to nuclear data (ND) uncertainties. ND uncertainties affect the computation of keff values through the calculation chain of BUCSS-R in three ways: They cause uncertainties in the SNF composition from the transport and depletion calculations, in the decay calculations as well as in the criticality calculations themselves. In a previous study [4], the uncertainties in the final keff values emerging from each of these steps have been quantified independent of each other. This caused an excess in conservatism since a full correlation of the uncertainty contributions had to be assumed between the depletion and criticality calculations. The present work, on the other side, pursues a consistent ND-UQ approach so that ND uncertainties in the criticality results could be determined in a more direct and realistic manner.
Other UQ methods in the area of criticality safety that are worthwhile mentioning here are, for instance, the Total Monte-Carlo method (TMC) based on reaction model parameters [5], non-intrusive and adaptive polynomial chaos techniques [6] and the linear perturbation method based on sensitivity calculations [7]. The code sequence called "Sampler" for uncertainty analysis with SCALE and XSUSA is an example of a statistical sampling method [8]. The direct application of measurements in the criticality models is the approach used in [9]. However, a consistent quantification of ND-related uncertainties in depletion and criticality calculations has not been so far subject of a dedicated analysis to the best of our knowledge and this is presented for the first time in this work.

Engineered Barrier System
The geological disposal scenario under consideration in this study is based on the repository concept for the disposal of SNF developed by NAGRA [10]. This concept envisages a horizontal emplacement of disposal canisters in excavated tunnels with Opalinus Clay as host rock and a bentonite backfilling of the voids around and between the canisters. The bentonite is then expected to gradually take up water from the host rock, swell and form a low permeability barrier around the canisters.
The preliminary disposal canister design essentially consists in a carbon steel cylinder, almost 5 m long, with an outer diameter of 105 cm and ready to fit 4 PWR FAs in separate carbon steel boxes inserted and welded, see Figure 1. No strong neutron absorbers are included by design.

Canister Loading
A bounding fuel case assessment in the context of criticality safety constraints for Swiss SNF loading in disposal canisters showed that loadings exclusively with UO2 FAs and a fuel enrichment of 4.94 w / 0 235U are the most critical ones [11]. The comparison covered loadings with all relevant PWR FA-designs: UO2 fuel, MOX as well as enriched reprocessed uranium (ERU), but also mixed loadings with FAs of different design. With UO2 FAs being the most critical type of loading, only this design was chosen for a more detailed ND-UQ in the present work.
The canister is considered to be loaded with 4 PWR-FAs having the same nuclear and mechanical design. The FAs correspond to Gösgen nuclear power plant (KKG) fuel design. Thus, a specific KKG FA representing the batch of 4.94 w / 0 235U enriched fuel and operated up to five cycles and, therefore, reaching the highest burnups at the end of life, was selected for the analysis. The irradiation history was reconstructed based on real plant operational data, see Section 3.1. In consequence, this study is based on a best estimate modelling approach for the fuel operating history.
The FA has a 15 × 15 array of fuel pins containing UO2 enriched to 4.94 w / 0 235U and 20 guide tubes. The ND-UQ was performed for the fresh FA as well as for assembly averaged discharge burn-ups of 17.61, 33.82, 50.47, 61.92 and 72.75 GWd/tHM, corresponding to five operation cycles.

Model Development and Assumptions
For the criticality calculations of this study, MCNP6 has been employed as one of the most validated Monte-Carlo (MC) codes for this purpose [12]. The fully 3-dimensional MCNP6 model consists in the disposal canister loaded with 4 FAs and surrounded by a 35 cm layer of water saturated bentonite clay, albeit the bentonite's effect on the criticality of the scenario is negligible. The remaining volume of the bentonite backfilling and the host rock formation are not anymore part of the MCNP model, i.e., modeled as void.
The Swiss regulations for a deep geological repository of high-level waste (HLW) require a proof of safety of the repository over a time span ("Nachweiszeitraum") of 1,000,000 years, explicitly including the criticality safety of the installation [13]. On the other hand, safe containment of the HLW by the disposal canisters has to be shown only for at least 1000 years. Although it is intended to exceed this minimum lifetime of the canisters by a suitable margin setting the lifetime design target to 10,000 years [10], and although it is not excluded that the canister sealing remains intact for many years longer, the containment of the HLW cannot be guaranteed by the canisters over the entire period of 1,000,000 years. For this reason, a breach of the canister is postulated and a full immersion of the FAs with water from the surrounding host rock is assumed by the MC model. Furthermore, the water properties such as composition, density and temperature were modelled as optimal for neutron moderation:

•
The water is modelled as pure H 2 O, ignoring the presence of diluted minerals. This leads to a better simulated neutron moderation and lower neutron absorption eventually resulting in higher keff values; • In previous preliminary assessments for canisters loaded with UO2 FAs, higher water densities have been found to correspond with higher keff values [2]. In view of the prevailing pressure and temperature regimes in the envisioned deep geological repository, a conservative water density of 1.005 g/cm 3 has, therefore, been set in the MCNP model; • The material temperatures of all model components are an item of required information in the neutron transport calculations, i.e., for the selection of ND libraries employed by the MCNP model. In the present model, these temperatures were assumed to be all equal to 293 K. This is a rather low temperature leading to a better neutron moderation in the MC-simulations.
A comparison of these conservative assumptions for the properties of water with more realistic conditions regarding their effects on the system criticality is currently in work at PSI. Criticality calculations are carried out for Opalinus Clay and bentonite porewaters instead of pure H 2 O, investigating their influence on keff simulation output.
Despite the postulated canister breach, the geometry of the canister is considered as intact in all simulations of this work. Also degradation issues and damages of the disposal canister and the FAs, which are inevitable over longer post-closure time scales, and corresponding changes in the geometry or the materials are not taken into account by the MC model. The BUCSS-R scheme applies BUC starting with the depletion calculation of a FA in the in-core burnup phase considering also shutdown cooling between cycles, followed by cooling during interim storage as well as cooling in the long term of final disposal. Based on detailed information about the irradiation conditions, this computational sequence allows for the determination of a realistic isotopic inventory in the FAs and its evolution over time. BUCSS-R is described in detail in [5]. Nevertheless, the elements of the methodology shall be briefly outlined here in order to provide the background for this study, see also Figure 2. Each step of the computational sequence is listed and outlined below, headed by the main code package(s) used for the current task: 1.
CASMO-SIMULATE The first task is a realistic estimation of operation histories, including in-core depletion and shutdown cooling between cycles. The calculations are based on reference CASMO-SIMULATE core models, developed and validated for all Swiss reactors and almost all operated cycles within the PSI core management system (CMSYS) platform [14][15][16][17].

BOHR
The detailed irradiation history can be extracted from SIMULATE3 results by means of an in-house code called BOHR [18]. State parameter values for every fuel assembly and axial elevation at numerous cycle instants are retrieved and corresponding CASMO5 [19] input files are produced.

CASMO5
New CASMO5 calculations are carried out based on BOHR output. Pin-wise and axially resolved isotopic compositions of the burnt fuel are determined. This step allows obtaining detailed spent nuclear fuel compositions, rather closely corresponding to the actual fuel assembly operating history.

SERPENT2
The SNF isotopic composition change due to decay is calculated by means of SER-PENT2 [20]. In this way, the long-term nuclide content evolution inside the loaded disposal canister can be determined. The choice of the SERPENT decay module was stipulated by its good performance and also due to the possibility to easily use different sources of the decay data, facilitating the different data libraries assessment [21].

COMPLINK
The SERPENT2 output is processed with the in-house code COMPLINK [22]. The isotopic composition data are automatically transferred into a pre-existing MCNP6-model template of the disposal canister loaded with FAs, producing MCNP6-input files ready to be executed [3].

MCNP6
The criticality calculation is finally performed with MCNP6 and the neutron multiplication factor keff is determined.
Nuclear data uncertainty is propagated through the execution of the calculation chain for n samples of perturbed nuclear data. The effect of ND uncertainties on criticality is then quantified by the standard deviation in the n-tuple of calculated keff values.
The BUCSS-R methodology has been used already for a preliminary criticality safety assessment of UO2 spent fuel assemblies from a Swiss pressurized water reactor (PWR) loaded into a disposal canister of NAGRA's design [11,23]. The results of this assessment were presented in the form of so called loading curves, indicating what minimum fuel assembly burnup is required for a fuel assembly with a given initial fuel enrichment to ensure that the effective neutron multiplication factor of the flooded disposal canister would comply with the imposed nuclear criticality safety (NCS) criterion.
A distinctive feature of PSI's CSE method in general has been outlined in [23]: "Normally, the fuel depletion calculations for BUC evaluations are done with a set of bounding parameters in terms of power density, fuel and coolant temperatures, densities, etc., so that the reactivity of the fuel at discharge for such conservative assumptions will be higher than the reactivity obtained with any possible real irradiation history. This path, however, was not fully adopted for the joint PSI/NAGRA BUCSS-R [ . . . ] project. The approach being utilised is different because real operational data could be employed (using available CASMO/SIMULATE3 validated core follow models at PSI) for all fuel assemblies operated in Switzerland. This allows estimating the loading curves on the basis of best-estimate assessments integrated with a conservative but rational treatment of uncertainties." A thorough comparison of the present method with others in the field of criticality safety is indicated at this point, however, corresponding efforts would certainly go beyond the scope of this work.

Uncertainty Propagation and Quantification
Following the BEPU concept, the fuel depletion and canister criticality calculations with BUCSS-R are furthermore complemented by UQs. CSEs in preparation of long-term geological disposal are in general affected by the following uncertainty sources: Operating conditions of the FA in the reactor core; • Design and technological parameters of the FAs and the disposal canister; • Radiation-induced changes in the FA geometry.
A comprehensive, although still preliminary, quantification of the uncertainties from all of these sources was summarised in [11]. ND-UQ results therein have been adopted from another study, also conducted at PSI [4]. Since ND are, in one format or another, essential input required for steps 3, 4 and 6 of the BUCSS-R scheme, several tools were developed and eventually used to quantify ND-related uncertainties stemming from each of those steps, see red frame in Figure 2. Using these tools, ND uncertainties can be propagated through the calculation chain by repeating the calculations for a large number of times, each time with a different set of ND, randomly sampled following information on the ND uncertainties provided in the selected ND library. The methods of the tools are based on stochastic sampling (SS), which has been explored extensively in the context of ND-uncertainty assessments for reactor physics applications and validation studies [24]. The UQ finally consists in the statistical analysis of the outcome of all the simulations. The following tools can be employed independently from each other: 1.
SHARK-X ND (i.e., cross sections, fission neutron yields ν and fission spectra) can be perturbed with the SHARK-X [16,[25][26][27] tool, which also takes care of the correct formatting for usage in CASMO5. First, ENDF covariance files are processed with NJOY into a number of energy groups from thermal energy up to 20 MeV. These matrices are then used to produce random realizations of ND by means of SS from joint probability distributions defined with normal distributions. SHARK-X has been used for a previous study [28] assessing uncertainties for Swiss LWR spent nuclear fuels due to ND, however, with a focus on other quantities of interest: assembly burnup, decay heat, neutron and gamma sources, as well as isotope inventory for realistic fuel histories.

2.
Fission yields sampling The UQ due to fission yields is similar to the case of the above ND. A stochastic sampling of evaluated fission yields is performed, based on fission yield uncertainties (as provided in the evaluated libraries), and on an in-house correlation matrix (as no fission yield correlations are provided in such evaluated libraries). Details on the correlation matrix can be found in [29].

ENDF2C
The impact of decay data uncertainties on keff results can be investigated with a modified version of the ENDF2C tool [30]. This routine uses stochastic sampling for a perturbation of the decay constants and branching ratios, whose uncertainties are contained and read from original ENDF decay data files. Furthermore, the physical relationship between the branching ratios is respected. The perturbed decay data are formatted for use in SERPENT2.
A previous study with ENDF2C [4] showed that the effect of decay data uncertainties on system criticality in the present context is negligible when compared with the effects of uncertainties in cross section or fission yield (FY) data. This finding is in line with other publications addressing decay data uncertainties in this respect [31][32][33]. An assessment of the effects of decay data uncertainties is, hence, not part of this work.

NUSS
The Nuclear data Uncertainty Stochastic Sampling method (NUSS) [5,25,34,35] applies perturbations in groups to pointwise nuclear data given in the specific ACE format which is used for criticality calculations with MCNP6. NUSS hereby follows the descriptions of ND sampling in [25]. The pointwise ND are divided into n segments with n being the specified number of energy-groups. Perturbation factors, i.e., the ratio between sampled and nominal ND, are then applied to each of these ND segments varying the data uniformly within each energy group. In order to ensure the consistency of the ACE file, total cross sections and the reaction cross sections are adjusted accordingly. With NUSS, perturbations of (n,n), (n,n ), (n,2n), (n,f ), (n,γ) cross sections and of ν, χ data can be applied.
An example of the functionality of NUSS and a comparison of SS with other UQ methods, e.g., the Total Monte Carlo (TMC) method and sensitivity based calculations, can be found in [5].

Consistent Sampling Approach Realisation
For each of steps 3, 4 and 6 of the BUCSS-R scheme, it is thus possible to introduce ND uncertainties into the calculation chain and to propagate them through it to evaluate their effect on the criticality of the system. In this way, ND-UQs have been carried out in an independent manner for each step, with BUCSS-R executions using perturbed ND sets for only one step at a time [4]. Different sets of sampled perturbation factors were applied in separate batches of BUCSS-R computations addressing either the depletion, decay or the criticality aspect. The present work, however, allowed for an all-encompassing ND-UQ featuring a simultaneous perturbation of ND input for step 3 (depletion) and step 6 (criticality), ignoring the negligible uncertainties in step 4 (decay). The simultaneous perturbation could be achieved by linking SHARK-X and NUSS to produce a consistent set of perturbed ND, formatted for usage in both, CASMO5 and MCNP6, respectively. In the following, this new approach is referred to as 'consistent ND-UQ'. Both, the results obtained with consistent ND-UQ and those of an independent ND-UQ are subject of a more detailed comparison and discussion in Section 4.

Consistent Sampling Approach Verification
Prior to this work, the capability of producing and using consistent sets of perturbed ND for both depletion (with Shark-X/CASMO5) and criticality (with NUSS/Serpent) calculations has been demonstrated. The Peach Bottom 2 pin-cell model of UAM in Hot Zero Power (HZP) conditions served as a test case [25]. The same Shark-X model was used while a consistent Serpent model was established based on the MCNP model used in [25]. The ENDF/B-VII.0 nuclear data library is used in both codes to calculate kinf of the pincell.
Various sources of covariance information were considered for testing the consistent ND sampling, namely ENDF/B-VII.1, ENDF/B-VIII.0 and JEFF-3.3. 300 perturbed ACE formatted files were generated with the NUSS scripts for each of the covariance libraries. Covariance matrices with a 19 energy group structure were applied in the case of ENDF/B-VIII.0 and JEFF-3.3 while for ENDF/B-VII.1 they had a 44 energy group structure. Only 235U, 238U and 1H have been considered for this test albeit more nuclides are available. The following reaction channels were perturbed, when available: (n,n), (n,n ), (n,2n), (n,f ), (n,γ), total ν and χ.
In order to provide the consistency of the ND sets employed in the CASMO5 and Serpent calculations, the perturbation factors produced by NUSS were used to generate the corresponding ND perturbations in Shark-X.
All Serpent calculations were performed with 100,000 histories per cycle, 200 active cycles and 20 inactive cycles. This leads to a kinf uncertainty of about 15 pcm, e.g., much less than the ND uncertainty. Furthermore, the seed number of each perturbed calculation is set to that of the nominal calculation.
Uncertainties in the kinf values, calculated with CASMO5 and Serpent on the basis of 300 sampled ND sets for each library, ENDF/B-VII.1, ENDF/B-VIII.0 and JEFF-3.3, are reported in Table 1. The kinf relative uncertainties are consistent, even though small discrepancies are observed, most likely due to the different mean values for kinf. The consistency is further demonstrated by ordering the Serpent ratio of the perturbed kinf with respect to the nominal value in increasing order and plotting it against the same quantity produced by Shark-X. This is shown in Figure 3a for ENDF/B-VII.1.  As expected, a slope of 1 is observed. Fluctuations around the unity slope are visible in Figure 3a. Further analysis on these fluctuations is carried by computing the "distance" between the Shark-X and Serpent perturbed kinf for a consistent perturbation; this distance is expressed as the ratio defined by Equation (1): The standard deviation of the distances is in the 50 pcm range and appears normally distributed, see Figure 3b. However, such fluctuations are somewhat larger than the inherent statistical spread expected from the ratio of two independent Serpent calculations, each calculation having a 15 pcm standard deviation (the resulting standard deviation, assuming that the correlation between both Serpent k-eff values in Equation (1) is zero, should then be √ 2 σ Serpent ≈ 20 pcm). The larger fluctuations are most likely caused by the different neutron transport approaches of Serpent and CASMO5, MC and deterministic method of characteristics (MoC).

Nuclear Data
This subsection provides essential information about the nuclear data considered in this work, details of the performed uncertainty propagation and all relevant options applied in the calculations. Tables 2 and 3, respectively, give an overview about the used ND libraries, code version numbers and execution options for each of the BUCSS-R steps and UQ modules. (1) Usage of perturbed nuclear data (ND) provided by SHARK-X. Actinides-plus-fission-product (AC + FP) isotopes with perturbed ND listed in Table 3.
( (1) Usage of perturbed ND provided by NUSS. AC + FP isotopes with perturbed ND listed in Table 3.  Table 3. Overview about materials and data used in nuclear data uncertainty quantification (ND-UQ) modules of this study. (2) Also group structure (44 energy groups) adopted from NUSS

. Nuclear Data Libraries
Both ND-library releases, ENDF/B-VII.0 and ENDF/B-VII.1 [36,37], were applied in the BUCSS-R calculations, see Table 2. However, covariance matrices only from ENDF/B-VII.1 were used for the ND perturbations, where they have been performed. In step 3 of BUCSS-R (depletion calculations with CASMO5), ND from ENDF/B-VII.0 were thus perturbed based on the covariance information from ENDF/B-VII.1. Nevertheless, this inconsistency is not of prime importance for the purpose of this study, the quantification of ND uncertainties.

Isotopic Composition and ND Perturbation
Two approaches, the actinides only (AC) and the actinides plus fission products (AC + FP) approach, are commonly used for the application BUC, taking into account either only major actinide nuclides or, in addition, also fission product nuclides in the criticality safety analysis. A detailed specification of these different sets of nuclides can be found in [38]. In order to cover the maximum number of nuclides as potential sources of ND uncertainties, the AC + FP approach has been chosen for the BUC criticality calculations in the present study. The final MCNP models, created within the BUCSS-R process, contain a total of 90 isotopes. For 56 isotopes thereof, covariance matrices were available from the ENDF/B VII.1 library so that ND-perturbations could be performed with NUSS/SHARK-X. All isotopes are listed in Table 4, including an overview of the perturbed reaction channels for each isotope. ND for the remaining isotopes were not perturbed and only nominal data were used. The reaction channel MT = 13 in Table 4 is not used by the official ENDF MT system, but is used here to represent a "combination" of elastic (MT = 2) and inelastic (MT = 4) scattering.
Uncertainties in thermal scattering data S(α,β) were not considered in this work. However, their contribution in a ND-UQ for keff is not negligible, although certainly not dominant. Uncertainties in thermal scattering data for H in H 2 O, for instance, account for a keff uncertainty in the range of about 150 pcm in the ICSBEP LEU-COMP-THERM benchmark LCT6 [39,40]. The possible application of NUSS using covariance information about thermal scattering data, once available, is discussed in [35].
FY data (MT = 454) were perturbed for 235U, 238U, 239Pu and 241Pu using an advanced normalization technique, see discussion in 4.3. Table 4. List of 90 isotopes present in the MNCP6 input files, 56 isotopes thereof with perturbed nuclear data and reaction channels. AC + FP nuclides highlighted in red and yellow, respectively, as specified in [38].

Results
200 BUCSS-R executions were performed for each burn-up step. Furthermore, decay times of 0 and 50,000 years were considered since the evolutions of keff values over a time span of 1,000,000 years exhibit a global and local maximum at these time points, regardless of the BU [4,11,41]. Note that this is not only true for the AC + FP case presented here, but also for a AC only case (without FP) in which the keff values' global maximum would be reached not at 0 but at 50,000 years of decay.
The resulting neutron multiplication factors are presented in Figure 4 (Note that the PWR canister filled with fresh fuel assemblies (BU = 0 GWd/tHM) and completely flooded with pure water can be critical. This is one of the main reasons why the disposal of fresh fuel makes no sense and is by no means intended. Nevertheless, the case of fresh fuel in a pure water-filled canister sets the parameters for the most reactive (bounding) scenario and, in this context, is merely of academic interest without any practical relevance for the final disposal. Several technical and procedural solutions to prevent the unlikely event of a criticality excursion are currently under consideration. For instance, the minimum average fuel assembly BU required to ensure that the keff value of the water-filled canister complies with the imposed criticality safety criterion can be given in the form of loading curves, which will be determined within the framework of the joint BUCSS-R project between PSI and NAGRA. A preliminary loading curve for FAs of KKG is reported in [11]). Shown is the cask eigenvalue keff, averaged over the 200 values obtained with an isotopic content corresponding to each burn-up and decay time, as well as error bars indicating the lowest and the highest calculated keff values, respectively. At zero burn-up, no depletion calculations had to be carried out. Since no noticeable change in the fresh UO2 fuel composition can take place over a decay time of 50,000 years, also no decay calculations were performed for this scenario, so that the values for 0 years of decay were simply adopted. For non-zero burnup, on the other side, keff values are lower after 50,000 years of decay than directly at discharge mainly due to 239Pu decay (T 1/2 ≈ 24,100 years).  Figure 5b. For better comparison with the previous method (Figure 5a, Section 4.1) and a detailed investigation of the prevailing correlations (Section 3.2), two more UQs have been performed taking into account only ND-related uncertainties in the depletion or criticality calculations, respectively. For each of these UQs, again 200 BUCSS-R executions were carried out at each burnup level and for 0 years or 50,000 years of decay. The nomenclature used is denoted and explained as follows: • σ CASMO ND (BU): Standard deviation in keff values due to nuclear data-related uncertainties associated with the spent fuel composition. These uncertainties are propagated through SNF compositions from depletion calculations with CASMO5, employing perturbed ND sets produced by SHARK-X. All these quantities are burn-up dependent, but for the sake of brevity the argument notation '(BU)' is omitted hereafter.
A detailed breakdown of the computational expense for an UQ based on 200 BUCSS-R executions is shown in Table 5. The depletion step of the UQ, for instance, requires a 30 min lasting CASMO5 calculation for each of the axial meshes times the number of samples (30 min × 41 × 200 = 170 d). For the criticality step, on the other hand, on average 17 h are needed for every considered BU value and each decay time, i.e., 0 or 50,000 years in our study (17 h × 5 × 2 × 200 = 1420 d). Compared to the other steps, depletion and criticality calculations with CASMO5 and MCNP6, respectively, consist in the computationally most expensive parts of the UQ. Nevertheless, by means of massive parallelization of the entire process, the effective time needed for the UQ can be reduced to a few weeks.

Comparison of New and Previous Approaches
The new approach of simultaneously and consistently perturbing ND input for depletion and criticality calculations made it possible to calculate σ TOTAL ND directly. Formerly, the total ND-related uncertainty in keff had to be determined indirectly using the independently computed quantities for σ CASMO ND and σ MCNP ND , applying Equation (2) for uncertainty propagation.
Since the correlation between σ CASMO ND and σ MCNP ND was unknown, a conservative correlation coefficient of r = 1 had to be assumed, as in a previous study [4]. In this case, the equation above simplifies to: Total ND-related uncertainties obtained via Equation (3) are shown in Figure 5a. σ TOTAL ND values are starting with about 800 pcm, steadily increasing with burnup to finally approximately 1030 pcm and 1140 pcm at 0 and 50,000 years of decay, respectively (It shall be noted here that in the course of the present work, a deficiency of the previous ND-UQ study reported in [4] has been discovered. Namely, the neutron multiplicity and spectrum, ν and χ data, have been unintentionally excluded from the random sampling of the ENDF/B-VII.1 ACE files with the NUSS tool in [4]. It is known from other studies, in particular those also performed with NUSS [25,28,34,42], that these parameters significantly contribute to the total keff uncertainties. In fact, the keff uncertainties (at one standard deviation level) for typical LWR UO2 lattices in water have been assessed in works [17,24,43] to be reaching values up to~600-700 pcm, while in [4] the σ MCNP ND value was erroneously estimated as only~380 pcm. Therefore, the ND uncertainties results reported in [4] shall be considered as noticeably underestimating the σ MCNP ND component). This clear upward trend for total uncertainties with increasing burnup level is driven by the σ CASMO ND component while σ MCNP ND is slightly decreasing over fuel depletion. In the present study, on the other hand, σ TOTAL ND results obtained with the new consistent ND-UQ reveal a slight downward trend with increasing burn-up, see Figure 5b. Total standard deviation amounts to about 800 pcm for zero burn-up, decreasing to a minimum of about 520 pcm for 0 years of decay and 590 pcm for 50,000 years of decay, at high burn-ups between 62 and 72 GWd/tHM. The new approach, intrinsically respecting correlations of ND-related uncertainties in depletion and criticality calculations, yields hence significantly lower σ TOTAL ND values for high burnups as compared with the previous method. The next subsection provides more details about this correlation aspect of the new ND-UQ.
Furthermore, it can be noticed that the results reported in Figure 5b are consistent with the previous results obtained at PSI for depletion of similar UO2 lattices with the CASMO code and applying stochastic sampling-based nuclear data uncertainties propagation with the PSI tool SHARK-X, when using the ENDF/B-VII.1 nuclear data and covariance matrix library [44,45].
The standard errors (SE) of the calculated standard deviations (σ TOTAL ND , σ CASMO ND and σ MCNP ND ) from both studies can be estimated using the following approximation, with N being the sample size [27].
This formula is valid assuming a normal distribution of the underlying population and a sufficiently large sample size N. For a sample size of 200, the equation becomes: Note that σ CASMO ND is equal to 0 pcm for a burn-up of 0 GWd/tHM, i.e., fresh fuel, as no depletion calculation has to be performed with CASMO5. For the same reason σ TOTAL ND coincides with σ MCNP ND in this case.

Correlation between Depletion and Criticality Perturbation
The Pearson correlation coefficient between ND-related uncertainties in depletion and criticality calculations can now be determined, analyzing the sets of keff values used for the calculation of σ CASMO ND and σ MCNP ND , respectively. This is possible due to the separate quantification of ND-related uncertainties stemming from either the depletion or the criticality step in the calculation chain and exploiting the consistent set of perturbation factors for both instances. The results, presented in Figure 6, show a significant anticorrelation which is even more pronounced at higher burnup. No Pearson correlation coefficient is calculated for zero burnup as there were no CASMO5 simulations. The standard error of the correlation values can be estimated using Equation (6) from [46]: The SE estimations for r range from about 5.6% (BU: 72.75 GWd/tHM, decay: 0 years) to 7.0% (BU: 17.61 GWd/tHM, decay: 50,000 years).

Discussion
This study considered a rather limited number of FA-averaged discharge burnups. This circumstance can be of relevance, especially regarding the quantification of ND-related uncertainties in depletion calculations. For low burnups, the few data points of σ CASMO ND at 0 GWd/tHM BU and the next discharge BU of 17.61 GWd/tHM in Figure 5 do not offer a complete picture of the ND-related uncertainties. In fact, assembly burnup calculations performed for UO2 PWR-FAs show the highest ND-related uncertainties in calculated assembly BUs at low exposures of about 10 GWd/tHM, as reported in [28]. Of course, those uncertainties cannot be translated directly into uncertainties in keff. Furthermore, with a maximum of about 2% they are rather moderate. Nevertheless, those findings indicate a possibly poor resolution of the present ND-UQ in this BU-range.
The UQ and further analyses in this study concentrated on a disposal canister loading with fuel enriched to 4.94 w / 0 235U. The conservatism of this approach in terms of ND uncertainties in the criticality calculations for the underlying SNF disposal scenario has been examined by performing ND-UQs also for lower enriched fuel and could eventually be verified. Compared to fuel with an enrichment of 3.2 w / 0 235U, for example, σ TOTAL ND is about 30 pcm (fresh fuel) to 90 pcm (at~32 GWd/tHM BU) higher in case of 4.94 w / 0 . ND play an essential role already in the first step of the BUCSS-R calculation scheme, the reconstruction of the operation history and corresponding depletion calculations based on CASMO-SIMULATE core models. ND uncertainties affecting the calculation parameters, however, were not considered by the current ND-UQ method, leaving room for its further development.
Knowledge of the precise reasons for the observed anti-correlations may be of general value for a better understanding of the underlying nuclear data effects. A follow up of this study with an investigation targeting specifically the anti-correlations from depletion and criticality calculations would thus be worthwhile.
Furthermore, it shall be noted that the question of how exactly to integrate the NDrelated uncertainties corresponding to the criticality calculations in the overall CSE methodology goes beyond the given paper. It can be just mentioned here in passing that normally the nuclear criticality safety criteria involve the so-called upper subcritical keff limits (USL), which are derived on the basis of associated validation studies using criticality safety benchmarks.
Thus, by their nature the USL values shall be correlated with the application case keff values. The stronger such correlation, the stronger can be the cancellation effect for the net ND-related uncertainties affecting the CSE and, for instance, for the loading curves in the context of the given study. More detailed explanation of this phenomenon and some quantitative examples can be found in [23,43].
Finally, as pointed out in [11], besides the BUCSS-R calculation scheme "an alternative way to derive the loading curves recently became available based on the exploitation of the advanced features of the additional CMS code "SNF", which can be seamlessly integrated into the PSI CMSYS system. A description of the new methodology developed at PSI outside the BUCSS-R project, which is based on the CASMO/SIMULATE/SNF/COMPLINK/MCNP system of codes and which was named the CS2M method has been presented in Reference [47]". The CS2M calculation scheme consists of two specific couplings: • CMSYS with SNF to reconstruct isotopics and calculate their evolution during decay; • SNF to MCNP using the in-house COMPLINK methodology to initialize the material compositions of all fuel rods in the MCNP models.
A discussion on potential advantages and verification of the new approach can be found in Section 7.2 "Alternative Advanced Way to Derive the Loading Curves under Development and Verification" of Reference [11].

Conclusions
PSI's quantification methodology for nuclear data uncertainties in the context of criticality safety evaluations for geological disposal of Swiss SNF has been enhanced by a new approach of perturbing ND input for both depletion as well as criticality calculations. The application of consistent sets of perturbation factors for the ND input in both depletion and criticality calculations lead to an improved understanding of the correlation of the two contributions to overall uncertainties in the keff results, eventually revealing an anticorrelation between them and a decreasing correlation-coefficient with higher burn-up. The suspected excessive conservatism of the previous method, applied e.g., in [4], in which a full correlation of ND-related uncertainties in keff stemming from depletion and criticality calculations has been assumed, could thus be confirmed.
The ND-UQ in this study was performed only for discharge burnups, leading to a rather coarse resolution of the uncertainty curves obtained over the BU domain. Consulting the literature on this topic, it seems that this is not a negligible issue, at least for low degrees of depletion. Nevertheless, no major deviations of the observed trends have to be expected. Furthermore, with the study's exclusive focus on canister loadings with PWR UO2 FAs, no statements can be made as to whether the findings presented on ND-related uncertainties are also valid for MOX assemblies. Moreover, the presented analysis was carried out in this detail only for one UO2 FA-design with an enrichment of 4.94 w / 0 . Nevertheless, a comparison with fuel of lower enrichment showed that this is a conservative approach and ND uncertainties for higher enriched fuel are for all burnup values well above those of lower enriched fuel.
As an integral part of the CSE methodology at PSI, the updated ND-UQ results will be employed by coming criticality assessments and in the derivation of loading curves for Swiss SNF assemblies.

Acknowledgments:
The authors would like to thank Madalina Wittel from NAGRA for her constructive review of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.