Evaluation of BWR Burnup Calculations Using Deterministic Lattice Codes SCALE-6.2, WIMS-10A and CASMO5

: The UK nuclear innovation programme supported by the government includes preparation for future ABWR construction. The UK has signiﬁcant expertise in building and operating gas-cooled nuclear reactors and some experience with PWRs, while there is limited knowledge in BWR technologies. Hence, an important aim of this work is to understand the discrepancies between codes to assess uncertainties in BWR lattice and depletion calculations, while identifying speciﬁc development demands to progress existing tools into extended applications. The objective of the study is to quantify the discrepancy between SCALE-6.2, CASMO5 and the UK WIMS-10A deterministic lattice code for BWR lattice and burnup modelling. Two models of BWR systems were considered for this new systematic comparison. They are a single BWR pin-cell with UO 2 fuel only, and a 3 by 3 array of BWR UO 2 fuel rods with gadolinia rod in the centre. Criticality over burnup was estimated for both models using these codes. Spectral indexes, number densities and neutron spectrum were compared for several burnup stages using SCALE-6.2 and WIMS-10A. The study showed that k inf obtained with CASMO5 was in a good agreement with the SCALE-6.2. A clear discrepancy in behaviour was observed between WIMS-10A and SCALE-6.2 as well as CASMO5.


Introduction
The increasing levels of Light Water Reactor (LWR) spent fuel (SF) and the requirement of the dry storage of SF for long term interim storage have resulted in industry exploring the potential of application of the burnup credit method (BUC) to Pressurized Water Reactor (PWR) and Boiling Water Reactor (BWR) SF. BWRs are of particular interest since they are widely deployed internationally, but their fuel characteristics are studied less than PWRs due to complexity of the modelling associated with their highly heterogeneous fuel assemblies, as well as core arrangements known for BWR [1,2]. This has been recently confirmed for an advanced BWR in [3] with a published core design which had input (and was reviewed) by experts in reactor physics at HITACHI. The fuel assembly and core complexities not only reflect reactor operation but also into the following handling of the spent fuel. Even if BWR fuel is less reactive and, thus, does not bring the same criticality safety concerns in relation to the BUC as PWR, the topic has gained significant traction in the last years [1,2].
Historically, the criticality safety analysis of SF dry storage facilities was performed with methods which treated SF in the system as fresh fuel without burnable poison-the so-called fresh fuel approach. This approach is considered to be highly conservative, and significantly overestimates the calculated reactivity of the cask filled with SF [2] as long as no burnable poison is apparent in the fuel. BUC tends to utilize realistic spent fuel characteristics and considers the depletion of fissile materials, as well as burnable poisons in the fuel assemblies during the reactor operation. For the fuel discharged at the targeted burnup, the application of BUC results in a reduction in the estimated reactivity of the system with SF when compared with an approach based on the assumption that the fuel is fresh. As was pointed out in [4], experimental measuring of the spent fuel reactivity is not convenient for BUC application; thus, it has to be demonstrated that the employed computational tools deliver robust results. To summarise, with the development of BUC, the capacity of the SF storage facilities can be expanded without affecting the safety margins by just relying on a more sophisticated modelling approach removing unnecessary conservative processes.
In the past few years, a number of studies were performed regarding the burnup credit method application to BWR spent fuel storage and transport. In [2], the computational benchmark was developed to explore the influence of fission products and minor actinides on reactivity in burnup credit calculations in comparison with those which use major actinides only. The study showed that reactivity in burnup calculations reduces by 70% due to major actinides, and by 30% due to additional nuclides. The application of the peak reactivity method to the transport and storage of BWR spent fuel was extensively investigated in [5]. Various parameters and modelling approaches were considered to determine their impact on peak reactivity. The study showed that void fraction and control blade insertion/position has the strongest influence, while other parameters such as fuel temperature or operation history have only a minor effect. The report [6] covers the effect of various axial moderator density distributions, control blade usage, and axial burnup distributions on BWR BUC beyond peak reactivity (for burnups 30-50 GWd/tU). The study showed that changes in axial moderator density and burnup profiles significantly affect the final storage cask reactivity, while the effect of control blade usage is not as large as was initially assumed.
As described above, the important part of developing BUC is a deep understanding of the fuel assembly behaviour in the reactor core over the burnup, as well as their robust modelling. The modern industrial modelling and simulation approach of steady state and transient modelling of the reactor core operations has two stages, which include the use of lattice and nodal codes. Nodal codes represent a reactor core as a set of homogenous coarse grid elements called nodes. Each node has homogenized, constant nuclear properties such as diffusion coefficient and material cross-sections. The nodal approach also does not reflect any geometric details of the fuel assembly such as fuel rods or control blades. This information has to be brought into the input constants based on higher-quality detailed modelling. In the applied sequence, lattice codes are used to prepare cross-section sets for each node providing multi-group transport informed constants as an input library. Fuel assemblies consist of a lattice of fuel rods in the radial direction, which is represented in a two-dimensional unstructured mesh which is the currently used industrial standard, while the axial layers of an assembly with unique fuel characteristics, such as enrichment or design, are represented in the nodal core simulator. In the cross section matrix for the core simulator, the sets of condensed and homogenized constants are usually reduced to two groups for LWRs simulations from the multi-group flux distribution delivering the fast and thermal energy group [7]. After cross-section sets are generated in 2D lattice code, they are transferred to 3D nodal code where the reactor core is modelled as coupled neutronics and thermal-hydraulics system. This approach is widely used across the industry and is considered to be an optimum in terms of minimizing computational burden whilst maintaining sufficient accuracy. The quality of cross-section sets prepared in lattice codes directly affects the accuracy of the reactor core simulation results. Thus, for the future application of the UK code system WIMS to a future reactor system ABWR, it is essential to evaluate the quality of the lattice calculation as a first step for the robust BWR modelling and simulation. The current study is used to investigate the quality and level of discrepancy between various deterministic lattice codes for burnup calculations, to evaluate the quality of the foreseen UK code system for the application to BWRs. The evaluation of different codes also supports the confidence into its quality for the cross-section preparation and, consequently, the following nodal analysis. The well-known lattice and burnup code system SCALE-6.2 [8] and the independent industrial standard code CASMO5 [9] are used to evaluate the quality of the WIMS-10A code [10] to answer the question 'Can WIMS be directly used for the modelling and simulation of future BWR reactors or is there some demand for upgrading of models to provide robust results?'. The answer to this question will provide guidance for further investment requests within the nuclear innovation programme of the Department for Business, Energy and Industrial Strategy [11].
Lattice code analysis is performed for two infinite BWR reactor systems, containing either a single UO 2 fuel rod (pin-cell) or a 3 by 3 array of 8 UO 2 fuel rods and a gadolinia-poisoned rod in the centre. The Gd-poisoned fuel rod is surrounded by uranium fuel rods to appropriately estimate its characteristics over the burnup.
Infinite multiplication factors calculated with SCALE-6.2 and WIMS-10A for burnups up to 57.5 GWd/tU were compared with each other and with the industrial standard code CASMO5. The limit of the burnup was chosen as the licensed maximum assembly average burnup of 55 GWd/tU [12] plus a safety gap. Deeper analysis, such as spectral indexes and number densities estimation, was performed using SCALE-6.2 to compare with WIMS-10A results, to deeper investigate the effect of the discrepancies on follow up nodal calculations and to evaluate the quality of the results and possible sources of discrepancies in detail.

Code Descriptions
Burnup calculations were performed using the deterministic lattice codes CASMO5, SCALE-6.2 and WIMS-10A. SCALE-6.2 is a modelling and simulation tool developed by the Oak Ridge Nuclear Laboratory (ORNL, USA) for criticality analysis, reactor and lattice physics calculations, and sensitivity and uncertainty analysis [8]. Multigroup (MG) neutron libraries based on ENDF/B-VII.0 and ENDF/B-VII.1 nuclear data libraries were used in the calculations. They have 238 and 252 energy group structure, respectively, and will be referred to as v7_238 and v7_252. Burnup calculations were performed with the TRITON module of SCALE-6.2. TRITON uses 2D NEWT module for the transport calculations and ORIGEN, for depletion analysis. The NEWT transport solver is based on the discrete ordinates method (S_N) on unstructured mesh. The S_N order was set to 6 in the current calculations.
WIMS-10A is a long term existing code for cell and lattice calculations for a broad range of reactor types, which has been developed and is currently upgraded by Wood plc, and their predecessors in the UK [10]. The MG neutron library is based on ENDF/B-VII.0 nuclear data. The improved decay data from ENDF/B-VII.1 library version were used in the calculations. WIMS-10A utilizes 172 energy group scheme. WIMS-10A performs transport calculations in 2D using the CACTUS module which implements Method of Characteristics (MOC), while burnup analysis is made in a BURNUP module.
CASMO5 is a lattice physics code for PWR and BWR fuel modelling [9] developed by Studsvik Scandpower. The code is based on the proprietary E7R0LIB neutron library with 586 energy groups which is used for the calculations of the BWR pin-cell. While for the 3 by 3 array with gadolinia model, the same library was used with the number of groups condensed to 19 [13] in the second step of the calculation. The CASMO5 transport solver is based on Method of Characteristics. The commercial version of CASMO5 was used for this research.

Model Descriptions for BWR Burnup Calculations
As already described, two models were created for the burnup calculations in application to BWR systems. The first one is a modern BWR pin-cell with a 4% enriched uranium dioxide fuel rod. The second is a hypothetical 3 by 3 array containing 8 BWR UO 2 fuel rods and one central rod with Energies 2020, 13, 2573 4 of 14 UO 2 mixed with gadolinia. In both cases, the depletion calculations were performed under a constant power assumption.

BWR Pin-Cell Model
The pin-cell of a BWR GE14 assembly containing UO 2 fuel rod was chosen for the first level of code comparison. A comprehensive description of its characteristics is provided in [2]. Table 1 shows geometrical and operational parameters used in the simulations. Figure 1 displays the SCALE-6.2 model of the BWR pin-cell which consists of fuel, clad and moderator. Reflective boundary conditions were applied to the system.

BWR Pin-Cell Model
The pin-cell of a BWR GE14 assembly containing UO2 fuel rod was chosen for the first level of code comparison. A comprehensive description of its characteristics is provided in [2]. Table 1 shows geometrical and operational parameters used in the simulations. Figure 1 displays the SCALE-6.2 model of the BWR pin-cell which consists of fuel, clad and moderator. Reflective boundary conditions were applied to the system.  3 10.5216 Moderator density, g/cm 3 0.6 Average power density, W/g 40 Figure 1. SCALE-6.2 model of BWR pin-cell.

BWR 3 by 3 Array Model with Gadolinia
The burnup calculations for a system containing burnable poison were evaluated based on an artificial BWR 3 by 3 array model with 8 UO2 fuel rods and one gadolinia containing fuel rod in the centre as shown on Figure 2. The UO2 fuel rods have the same specification as for BWR pin-cell case described in 3.1. The central poisoned rod consists of 4% UO2 fuel with 5% of Gd2O3 burnable absorber. It has the same geometrical and temperature parameters as UO2 fuel rod (Table 1). (U,Gd)O2 fuel density is 10.24 g/cm 3 . Figure 2 displays the BWR 3 by 3 array model created in SCALE-6.2 with the additional discretization introduced into the (U,Gd)O2 fuel rod. Discretization was introduced to capture the so called "onion-skin" effect [7] in a Gd-poisoned fuel rod. This effect occurs since gadolinia isotopes have a large absorption cross-section which leads to the instant absorption of thermal neutrons entering the surface of the rod. It causes the Gd-poisoned rod to deplete by layers [7]. Hence, to model this effect, the (U,Gd)O2 fuel rod was divided into 10 annular regions of equal area. Again, reflective boundary conditions were applied to the model.

BWR 3 by 3 Array Model with Gadolinia
The burnup calculations for a system containing burnable poison were evaluated based on an artificial BWR 3 by 3 array model with 8 UO 2 fuel rods and one gadolinia containing fuel rod in the centre as shown on Figure 2. The UO 2 fuel rods have the same specification as for BWR pin-cell case described in 3.1. The central poisoned rod consists of 4% UO 2 fuel with 5% of Gd 2 O 3 burnable absorber. It has the same geometrical and temperature parameters as UO 2 fuel rod (Table 1). (U,Gd)O 2 fuel density is 10.24 g/cm 3 . Figure 2 displays the BWR 3 by 3 array model created in SCALE-6.2 with the additional discretization introduced into the (U,Gd)O 2 fuel rod. Discretization was introduced to capture the so called "onion-skin" effect [7] in a Gd-poisoned fuel rod. This effect occurs since gadolinia isotopes have a large absorption cross-section which leads to the instant absorption of thermal neutrons entering the surface of the rod. It causes the Gd-poisoned rod to deplete by layers [7]. Hence, to model this effect, the (U,Gd)O 2 fuel rod was divided into 10 annular regions of equal area. Again, reflective boundary conditions were applied to the model.

Meshing
The spatial discretization of the cell in burnup calculations in deterministic codes is necessary to produce reliable results. The research in [14] outlines that the difference between kinf in LWR-burnup calculations can vary from −400 to +400 pcm for the basic (3 regions) and reference (328 regions) discretization cases.
In the current research, the standard CASMO5 settings have been used for depletion calculations for both models. CASMO5 is well validated for BWR and PWR cases [9] and has been used, in conjunction with the nodal code SIMULATE, for modelling ABWR nuclear reactor cores [3]. Sensitivity studies were performed in WIMS-10A to optimize the number of tracks used in the method of characteristics solution scheme and meshing. This resulted in 9 polar and 17 azimuthal angles used in the WIMS calculation with the XFINE and YFINE options also enabled in the WIMS-10A to subdivide each unit cell of the problem into 100 equi-volumes to increase fidelity. Sensitivity analysis in SCALE 6.2 showed that 16 by 16 grid (Figure 3) for the pin-cell is enough to produce accurate results. Furthermore, like in the SCALE-6.2 calculation, the UO2 rods were modelled as a single material region in WIMS-10A, but the UO2-Gd2O3 rods were subdivided into 10 equi-volume materials.

Infinite Multiplication Factor
The infinite multiplication factor (kinf) for both models was calculated from 0 to 57.5 GWd/tU range of burnup. CASMO5 and WIMS-10A results were compared against SCALE-6.2 with ENDF-VII.1 library. The discrepancy in kinf between codes was defined as:

Meshing
The spatial discretization of the cell in burnup calculations in deterministic codes is necessary to produce reliable results. The research in [14] outlines that the difference between k inf in LWR-burnup calculations can vary from −400 to +400 pcm for the basic (3 regions) and reference (328 regions) discretization cases.
In the current research, the standard CASMO5 settings have been used for depletion calculations for both models. CASMO5 is well validated for BWR and PWR cases [9] and has been used, in conjunction with the nodal code SIMULATE, for modelling ABWR nuclear reactor cores [3]. Sensitivity studies were performed in WIMS-10A to optimize the number of tracks used in the method of characteristics solution scheme and meshing. This resulted in 9 polar and 17 azimuthal angles used in the WIMS calculation with the XFINE and YFINE options also enabled in the WIMS-10A to subdivide each unit cell of the problem into 100 equi-volumes to increase fidelity. Sensitivity analysis in SCALE 6.2 showed that 16 by 16 grid (Figure 3) for the pin-cell is enough to produce accurate results. Furthermore, like in the SCALE-6.2 calculation, the UO 2 rods were modelled as a single material region in WIMS-10A, but the UO 2 -Gd 2 O 3 rods were subdivided into 10 equi-volume materials.

Meshing
The spatial discretization of the cell in burnup calculations in deterministic codes is necessary to produce reliable results. The research in [14] outlines that the difference between kinf in LWR-burnup calculations can vary from −400 to +400 pcm for the basic (3 regions) and reference (328 regions) discretization cases.
In the current research, the standard CASMO5 settings have been used for depletion calculations for both models. CASMO5 is well validated for BWR and PWR cases [9] and has been used, in conjunction with the nodal code SIMULATE, for modelling ABWR nuclear reactor cores [3]. Sensitivity studies were performed in WIMS-10A to optimize the number of tracks used in the method of characteristics solution scheme and meshing. This resulted in 9 polar and 17 azimuthal angles used in the WIMS calculation with the XFINE and YFINE options also enabled in the WIMS-10A to subdivide each unit cell of the problem into 100 equi-volumes to increase fidelity. Sensitivity analysis in SCALE 6.2 showed that 16 by 16 grid (Figure 3) for the pin-cell is enough to produce accurate results. Furthermore, like in the SCALE-6.2 calculation, the UO2 rods were modelled as a single material region in WIMS-10A, but the UO2-Gd2O3 rods were subdivided into 10 equi-volume materials.

Infinite Multiplication Factor
The infinite multiplication factor (kinf) for both models was calculated from 0 to 57.5 GWd/tU range of burnup. CASMO5 and WIMS-10A results were compared against SCALE-6.2 with ENDF-VII.1 library. The discrepancy in kinf between codes was defined as:

Infinite Multiplication Factor
The infinite multiplication factor (k inf ) for both models was calculated from 0 to 57.5 GWd/tU range of burnup. CASMO5 and WIMS-10A results were compared against SCALE-6.2 with ENDF-VII.1 library. The discrepancy in k inf between codes was defined as: where d is the discrepancy between code i and SCALE-6.2 code in pcm; k SCALEv7_252 is the k inf calculated with SCALE-6.2 (v7_252); k i is the k inf calculated with code i; i is one of the considered codes CASMO5, WIMS-10A or SCALE-6.2 (v7-238). Figure 4 (left) shows the discrepancy in k inf for BWR pin-cell model obtained using the SCALE-6.2 (v7-238), WIMS-10A and CASMO5 codes in relation to SCALE-6.2 (v7-252). It can be seen that results obtained using CASMO5 are in good agreement with SCALE-6.2 and v7_252 library, which gives an independent evaluation of the quality of the used reference code. The discrepancy between SCALE-6.2 (v7_252 library) and CASMO5 codes lies in the interval from −30 to +140 pcm. The infinite multiplication factor calculated by SCALE-6.2 with v7_238 library has a negative bias of around 300 pcm from SCALE-6.2 (v7_252 library), with discrepancy lying in the interval from −170 to −330 pcm. This bias is obviously caused by the usage of the somewhat older nuclear data libraries which have also slightly different energy structure, since the used solver is identical for these two cases. WIMS-10A and SCALE-6.2 (v7_252 library) comparison showed that the discrepancy between codes systematically increases with burnup and lies in the interval from −850 to 270 pcm, with a clear negative gradient over the whole burnup period.    gives information about the discrepancy in k inf between the considered codes for BWR 3 by 3 array with one gadolinia rod. In this much more challenging case, larger differences are observed. According to the results, there is an initial bias between SCALE-6.2 (v7_238), CASMO5 and SCALE-6.2 (v7_252) k inf , which is equal to −110 and 230 pcm, respectively. The discrepancy between given codes over burnup lies in the intervals from −325 to −110 pcm (v7_238 library) and from 160 to 280 pcm (CASMO5). Thus, it can be considered that SCALE-6.2 and CASMO5 codes are also in good agreement with each other for the much more challenging systems with burnable absorbers. As for the BWR pin-cell model, a large discrepancy has been observed between WIMS-10A and SCALE-6.2 Energies 2020, 13, 2573 7 of 14 results. It varies from −1600 to 200 pcm for v7_252 library of SCALE-6.2 and increases with burnup and requires some further, in depth investigation. In this second case, the discrepancy for the WIMS-10A code once more grows with a strong negative gradient. Figure 4 also shows that the discrepancy in k inf between WIMS-10A and SCALE-6.2 significantly increases for the case with the burnable absorber which is more challenging for the solver. Consequently, these large discrepancies will propagate through the cross section set into the core simulator and can lead to shifts in the moderator density and burnup profiles in the full core analysis, which will affect the spread of the burnup distribution over the reactor core.
The university available, SCALE-6.2 and CASMO5 standard codes indicated a very good agreement in k inf , so an investigation between WIMS-10A and SCALE-6.2 was performed for the following sections.
Moderator density in the BWR reactor core varies significantly in terms of distance from the bottom of the core as it can be seen in Figure 5. Thus, it is necessary to identify the impact of the moderator density on the level of discrepancy between codes. Burnup calculations for the BWR pin-cell model with several typical moderator densities such as 0.2, 0.3, 0.4 and 0.6 g/cm 3 were performed in the follow up analysis. Figure 6 shows that the level of discrepancy between SCALE-6.2 and WIMS-10A increases slightly when the moderator density decreases. This means that for higher density configurations, closer to PWR operation conditions, WIMS-10A performs better than for lower density BWRs. However, all considered cases show the same tendency of a clear negative gradient which is much stronger than the effect of the water density changes. Thus, the water density variation does not have a significant effect on the systematic character of the discrepancy occurring over burnup. In this case, the typically considered challenge of BWRs, the wide spread of the water density and the resulting massive change in the criticality of the fuel assemblies are captured by both codes very well.
We consider SCALE-6.2 as a reference solution for the further steps of our analysis since it is the currently used university standard tool. It is obvious that there is a good agreement between the industry standard code CASMO-5 and SCALE-6.2. In addition, the code is very well validated for different types of light water reactors and for burnup calculations against Monte-Carlo solutions as given in the user manual [8] and in an extensive report [15] as well as against experimental data published in Nuclear Engineering and Design [16].

Spectral Indexes Method
Spectral indexes (SI) are a method often used in the dosimetry and experimental analysis to get a deeper insight into the differences between model and experiment [17] without taking too much weight into the absolute value of specific reaction rates. SI is described as a ratio between averaged over the same spectrum cross-sections of the isotope of interest and reference one [18]. The method helps to reduce the effect of systematic biases, e.g., in the cross section set, and is often used for code quality and nuclear data library assessment. We will follow the idea of the SI method to identify the specific differences of the reaction rates produced by the different lattice codes to judge the consequence on the full core simulations. In [17], spectral index was defined as the ratio of two microscopic fission cross-sections averaged over the same neutron distribution: where i and r are the isotopes of interest and the reference one respectively; σ i and σ r are averaged microscopic cross-sections of the isotope i and r, respectively; σ i and σ r are microscopic cross-sections of the isotope i and r respectively; ϕ is the neutron flux; E is an energy. Since the depletion calculations were performed under a constant power approach, the neutron flux estimated by each code is not necessarily identical. Hence, the averaged microscopic cross-sections cannot be used for the SI evaluation since they contain integrated neutron flux in the denominator which will not necessarily cancel out in (2). For the current study, we will define SI through the ratio of reaction rates between the reference isotope and the isotope of interest as: where Σ i and Σ r are macroscopic cross-sections of the isotope i and r, respectively. The analysis was performed for SCALE-6.2 and WIMS-10A codes. SI were estimated for the both models over total neutron spectrum for fission cross-sections of: U-238 to U-235 or F28/F25; Pu-239 to U-235 (F39/F25); Pu-240 to Pu-239 (F40/F39); Pu-241 to Pu-239 (F41/F39); Pu-242 to Pu-239 (F42/F39). The difference between spectral indexes was calculated as: where SI WIMS (SI SCALE ) is the spectral index obtained with WIMS-10A and SCALE-6.2, respectively. Figure 7 depicts the difference between the spectral indexes obtained using WIMS-10A and SCALE-6.2 codes and calculated for both considered models. It should be noted that for both models, U-238 and Pu-240 fission cross-sections are overestimated by WIMS-10A in comparison with SCALE-6.2, while for Pu-242 the cross-sections are underestimated. For the BWR pin-cell model the difference in the F40/F39 spectral index varies from 1.8% to 3.7%, for the F28/F25, from 0.4% to 1.5% and for the F42/F39, from −1.9% to 0.2%. For the BWR 3 by 3 array model, the difference in the F40/F39, F28/F25 and the F42/F39 spectral indexes gradually rises as burnup and varies from 2.3% to 3.3%, 0.5% to 2% and −1.5% to 0.4% respectively.
Energies 2020, 13, x FOR PEER REVIEW 9 of 13 and the F42/F39 spectral indexes gradually rises as burnup and varies from 2.3% to 3.3%, 0.5% to 2% and −1.5% to 0.4% respectively. The difference in the F39/F25 spectral index rises with burnup in both models and varies roughly from −0.7% to 0.4%, while for the F41/F39 spectral index it has an opposite trend, and lies in the interval from approximately −0.4% to 0.4% for both models.
It should be noted that adding burnable poison to the system does not affect the main trends in spectral index difference as a function of burnup when compared to the UO2 fuel-containing system. Besides, the effect of the burnable poison is limited to the initial phase up to ~20 GWd/tU until the major part of the burnable poison is diminished. At this phase, Gd is competing in absorbing neutrons which affects Pu-239 accumulation and in addition the destruction of Pu-239 through the neutron absorption, leading to a reduced buildup of all higher plutonium isotopes. In the short term, it suppressed the deep in F39/F25 SI at the burnup 5 GWd/tU and, in the long term, decreased the difference between codes for this SI. Furthermore, in contrast with other SIs, F41/F39 sees almost no change for the model with burnable absorber in comparison with UO2 only.
The effect of spectral index discrepancy depends on the amount of isotopes appearing in the system. For example, inaccurate prediction of the Pu-240 fission cross section will lead to the incorrect estimation of its content, which will influence the capture reactions and thus the accumulation of isotopes along the breeding chain will not be preserved. It can significantly impact the final isotopic concentrations and the criticality of the spent fuel, especially in the case of modern high burnup fuel. Figure 7. Difference in % between spectral indexes obtained with the WIMS-10A and SCALE-6.2 codes for the BWR pin-cell (right) and BWR 3 by 3 array (left) models.

Isotopic Composition
The number densities of several main actinides for both models were calculated at the same set of burnup points using the SCALE-6.2 and WIMS-10A codes with Gd-155 for the BWR 3 by 3 array model. The difference between the number densities as a function of burnup for both models is shown in Figure 8. It was estimated as a subtraction of SCALE-6.2 number densities from WIMS-10A ones at every burnup point of interest.
Already at first glance, it is clear that in both cases the two leading fissile isotopes U-235 and Pu-239 show different trends, with the Pu-239 content lower through the whole burnup period and the U-235 content lower after a short initial peak. This can partly explain the observed differences in the infinite multiplication factor. The trend to the increasing negative discrepancy correlates to a slightly lower amount of the main fissile isotopes. However, it would be unwise to claim here that the burnup model has deficiencies, since the formation and destruction of the leading Pu isotope 239 is a complex balance between power, isotopic composition, capture cross sections, neutron flux and neutron spectrum. All these parameters can be easily influenced by other parameters, e.g., discretization and/or self-shielding [14]. The difference in the F39/F25 spectral index rises with burnup in both models and varies roughly from −0.7% to 0.4%, while for the F41/F39 spectral index it has an opposite trend, and lies in the interval from approximately −0.4% to 0.4% for both models.
It should be noted that adding burnable poison to the system does not affect the main trends in spectral index difference as a function of burnup when compared to the UO 2 fuel-containing system. Besides, the effect of the burnable poison is limited to the initial phase up to~20 GWd/tU until the major part of the burnable poison is diminished. At this phase, Gd is competing in absorbing neutrons which affects Pu-239 accumulation and in addition the destruction of Pu-239 through the neutron absorption, leading to a reduced buildup of all higher plutonium isotopes. In the short term, it suppressed the deep in F39/F25 SI at the burnup 5 GWd/tU and, in the long term, decreased the difference between codes for this SI. Furthermore, in contrast with other SIs, F41/F39 sees almost no change for the model with burnable absorber in comparison with UO 2 only.
The effect of spectral index discrepancy depends on the amount of isotopes appearing in the system. For example, inaccurate prediction of the Pu-240 fission cross section will lead to the incorrect estimation of its content, which will influence the capture reactions and thus the accumulation of Energies 2020, 13, 2573 10 of 14 isotopes along the breeding chain will not be preserved. It can significantly impact the final isotopic concentrations and the criticality of the spent fuel, especially in the case of modern high burnup fuel.

Isotopic Composition
The number densities of several main actinides for both models were calculated at the same set of burnup points using the SCALE-6.2 and WIMS-10A codes with Gd-155 for the BWR 3 by 3 array model. The difference between the number densities as a function of burnup for both models is shown in Figure 8. It was estimated as a subtraction of SCALE-6.2 number densities from WIMS-10A ones at every burnup point of interest.
Already at first glance, it is clear that in both cases the two leading fissile isotopes U-235 and Pu-239 show different trends, with the Pu-239 content lower through the whole burnup period and the U-235 content lower after a short initial peak. This can partly explain the observed differences in the infinite multiplication factor. The trend to the increasing negative discrepancy correlates to a slightly lower amount of the main fissile isotopes. However, it would be unwise to claim here that the burnup model has deficiencies, since the formation and destruction of the leading Pu isotope 239 is a complex balance between power, isotopic composition, capture cross sections, neutron flux and neutron spectrum. All these parameters can be easily influenced by other parameters, e.g., discretization and/or self-shielding [14].
The difference in number densities between WIMS-10A and SCALE-6.2 behaves similarly for U-235, Pu-240 and Pu-241 nuclides in both models. The behaviour changes for Pu-239 nuclide when instead of decreasing, the difference after 30 GWd/tU as in case of BWR pin-cell, it starts to gradually decrease slightly for the BWR 3 by 3 array model.
From initial time and up to 15 (10) GWd/tU in the case of the BWR pin-cell (3 by 3 array), U-235 is consumed more slowly in WIMS-10A than in SCALE-6.2; however, at around 30 GWd/tU the trend changes to the opposite effect and U-235 is being consumed faster.
As for the Pu-239 and Pu-241 isotopes, they accumulate constantly slower in WIMS-10A (underestimation) for both the pin-cell and 3 by 3 array models, while Pu-240 starts to build up faster from roughly 30 GWd/tU burnup.
Overall, the difference in U-235 number densities in WIMS-10A and SCALE-6.2 is almost twice as small for the model with the burnable absorber for burnups up to 20 GWd/tU. It can be related to the high absorption of neutrons by Gd at the beginning of the assembly life, which leads to a reduction in the thermal neutron flux and to the lower depletion of U-235 and, hence, a difference in the number densities.
As was mentioned in [2], 70% of the reactivity reduction in the storage cask system with BWR spent fuel comes from major actinides. Hence, it is important to have a reliable, robust estimate of their number densities for the high target burnups. The given results show that for fuel assemblies with typically discharged burnups (40-50 GWd/tU), the difference in number densities between SCALE-6.2 and WIMS-10A for major actinides becomes significant in both cases. For example, number densities for Pu-239 differ by around 1.5% at the EoL in the case of BWR pin-cells. This would account to a 180 pcm difference in k inf , while the change in number densities of the other observed isotopes would account for the additional 150 pcm.
Energies 2020, 13, 2573 11 of 14 their number densities for the high target burnups. The given results show that for fuel assemblies with typically discharged burnups (40-50 GWd/tU), the difference in number densities between SCALE-6.2 and WIMS-10A for major actinides becomes significant in both cases. For example, number densities for Pu-239 differ by around 1.5% at the EoL in the case of BWR pin-cells. This would account to a 180 pcm difference in kinf, while the change in number densities of the other observed isotopes would account for the additional 150 pcm.

Neutron Spectrum
The neutron spectrum in the form of neutron flux per unit lethargy was investigated in the SCALE-6.2 and WIMS-10A codes for both the BWR pin-cell and BWR 3 by 3 array models. Figure 9 and Figure 10 present the flux over energy distribution obtained in both codes for two burnup points-0 and 57.5 GWd/tU (beginning and end of life). The results show that neutron flux in the thermal energy region is in good agreement in SCALE-6.2 and WIMS-10A for all considered models and burnups. This is due to the comparable, almost identical, energy group number and structure in the thermal region used in both codes. The epithermal and fast energy regions have slight depressions in the neutron flux caused by absorption resonances which are much more pronounced in SCALE due to the higher number of energy groups and due to the different energy bin structure in these regions. SCALE-6.2 has approximately 1.5 times more energy groups in those regions in comparison

Neutron Spectrum
The neutron spectrum in the form of neutron flux per unit lethargy was investigated in the SCALE-6.2 and WIMS-10A codes for both the BWR pin-cell and BWR 3 by 3 array models. Figures 9 and 10 present the flux over energy distribution obtained in both codes for two burnup points-0 and 57.5 GWd/tU (beginning and end of life). The results show that neutron flux in the thermal energy region is in good agreement in SCALE-6.2 and WIMS-10A for all considered models and burnups. This is due to the comparable, almost identical, energy group number and structure in the thermal region used in both codes. The epithermal and fast energy regions have slight depressions in the neutron flux caused by absorption resonances which are much more pronounced in SCALE due to the higher number of energy groups and due to the different energy bin structure in these regions. SCALE-6.2 has approximately 1.5 times more energy groups in those regions in comparison to WIMS-10A leading to more detailed visibility of resonance structures. For the case with the burnable absorber, neutron flux in the thermal region at the BoL is lower than for the pin-cell configuration. to WIMS-10A leading to more detailed visibility of resonance structures. For the case with the burnable absorber, neutron flux in the thermal region at the BoL is lower than for the pin-cell configuration.

Conclusions
Development and application of burnup credit methods for Boiling Water Reactors require a deep understanding of the fuel assembly behaviour in the reactor core. Lattice physics computations are an essential part of this analysis and thus reliable and robust code results are an essential part of building trust in the future application of the method. Therefore, it is important to estimate the quality of the lattice codes and identify the level of discrepancy between them. For this purpose, the UK code system WIMS-10A is compared with the US code system SCALE-6.2, and the leading independent industrial code CASMO5 for BWR burnup calculations. Two different models were investigated in the present study for boiling water reactor systems: a pin-cell and a 3 by 3 array, with gadolinia

Conclusions
Development and application of burnup credit methods for Boiling Water Reactors require a deep understanding of the fuel assembly behaviour in the reactor core. Lattice physics computations are an essential part of this analysis and thus reliable and robust code results are an essential part of building trust in the future application of the method. Therefore, it is important to estimate the quality of the lattice codes and identify the level of discrepancy between them. For this purpose, the UK code system WIMS-10A is compared with the US code system SCALE-6.2, and the leading independent industrial code CASMO5 for BWR burnup calculations. Two different models were investigated in the present study for boiling water reactor systems: a pin-cell and a 3 by 3 array, with gadolinia poison rod in the centre. For both models, the SCALE-6.2 and CASMO5 results of burnup calculations demonstrated very good agreement. However, bias is observed for the infinite multiplication factor estimated with the ENDF/B-VII.0-based neutron library of SCALE-6.2 for both models, and with ENDF/B-VII.1, for 3 by 3 array model only.
A large discrepancy for higher burnups is detected between WIMS-10A and SCALE-6.2 along with CASMO5 k inf values. The difference between WIMS-10A and SCALE-6.2 (v7_252 library) increases with burnup and reaches a maximum of −850 and −1600 pcm for the BWR pin-cell and 3 by 3 array model, respectively. Thus, WIMS-10A generally underestimates the infinite multiplication factor in comparison with the SCALE-6.2 and CASMO5 codes. Furthermore, it was demonstrated that the discrepancy between the infinite multiplication factors calculated by WIMS-10A and SCALE-6.2 rises as the moderator density decreases, but these differences are almost negligible compared to the discrepancies over burnup. Since the SCALE-6.2 and CASMO5 results were in a good agreement with each other for k inf , further evaluation of WIMS-10A was based on the SCALE-6.2 code system to provide a deeper understanding of possible sources of the discrepancies.
The spectral indexes analysis showed that for both models, the U-238 and Pu-240 fission cross-sections are overestimated by WIMS-10A in comparison with SCALE-6.2, while the fission cross-section for Pu-242 is underestimated. The number density analysis detected that U-235 burns faster in the WIMS-10A model for burnups higher than 30 GWd/tU, while Pu-239 builds up slower on the whole range of considered burnups. In both cases, the trend in isotopic composition variance is similar, which may result from the difference in the burnup chains included in the solvers.
The neutron spectrum analysis revealed that the WIMS-10A and SCALE-6.2 neutron flux distribution over energy for both models have similar behaviour, with small depressions in the epithermal and fast energy regions due to the different energy group number and structure of those regions.
Based on the above, SCALE-6.2 provides more conservative results for k inf estimation in comparison with WIMS-10A. For burnups above 30 GWd/tU, WIMS-10A in both cases underestimates number densities of important fissile materials. Since CASMO5 and SCALE-6.2 codes produce close results, it makes SCALE-6.2 a more reliable tool for burnup computations for LWRs as required for burnup credit analysis. Hence, SCALE-6.2 will be used as a basic lattice code for further burnup credit development, and the deviations have been fed back to the code developer of WIMS.
The considered models and results were not compared against experiments, hence it is difficult to judge in the differences between WIMS-10A and SCALE-6.2 and CASMO5 in a form to decide which code would be closest to a real experiment. However, due to a very good validation base and the wide industrial use, there is high confidence in the CASMO5 results, and the code-to-code comparison offers access to very detailed results for comparison which would hardly be possible on the basis of experimental data. The possible reasons for the discrepancies can vary from the nuclear data inconsistencies-for instance, in the code specific cross-section or decay data libraries-to the differences in computational methods utilized by the considered codes. Even if the master libraries are based on the same nuclear data library-such as ENDF/B-VII.0 in a case of WIMS-10A and SCALE-6.2-the way they were processed to the code master library can vary, which may lead to the difference in the depletion results. The difference between the master libraries could be estimated by a comparison of specific nuclide reactions from those libraries.
With the BEIS nuclear innovation programme in mind, we recommend a deeper investigation into the discrepancies, working with the code developer to improve the trust in the UK lattice code system WIMS-10A for the Boiling Water Reactor application if this code system is foreseen for future applications.