Multiphysics Modeling and Validation of Spent Fuel Isotopics Using Coupled Neutronics/Thermal-Hydraulics Simulations

Multiphysics coupling of neutronics/thermal-hydraulics models is essential for accurate modeling of nuclear reactor systems with physics feedback. In this work, SCALE/TRACE coupling is used for neutronic analysis and spent fuel validation of BWR assemblies, which have strong coolant feedback. 3D axial power profiles with coolant feedback are captured in these advanced simulations. ,e methodology is applied to two BWR assemblies (2F2DN23/SF98 and 2F2D1/F6), discharged from the Fukushima Daini-2 unit. Coupling is performed externally, where the SCALE/T5-DEPL module transfers axial power data in all axial nodes to TRACE, which in turn calculates the coolant density and temperature for each of these nodes.Within a burnup step, the data exchange process is repeated until convergence of all coupling parameters (axial power, coolant density, and coolant temperature) is observed. Analysis of axial power, criticality, and coolant properties at the assembly level is used to verify the coupling process. ,e 2F2D1/F6 benchmark seems to have insignificant void feedback compared to 2F2DN23/SF98 case, which experiences large power changes during operation. Spent fuel isotopic data are used to validate the coupling methodology, which demonstrated good results for uranium isotopes and satisfactory results for other actinides. ,is work has a major challenge of lack of documented data to build the coupled models (boundary conditions, control rod history, spatial location in the core, etc.), which encourages more advanced methods to approximate such missing data to achieve better modeling and simulation results.


Introduction
e ability to accurately and efficiently predict end of life (EOL) fuel composition is essential in ensuring safe and economic management of spent nuclear fuel. In addition, verification, validation, and uncertainty quantification of neutronics/depletion codes with spent fuel data are important in order to evaluate the accuracy of such physics codes. Unfortunately, fuel depletion calculations are expensive to perform, making brute-force uncertainty analysis methods difficult to apply [1]. In single physics or decoupled simulations, it can often be difficult to directly simulate a particular experiment lacking detailed data. For example, time-dependent axial coolant density distributions are needed for light water reactor simulations, especially for boiling water reactors (BWRs) but data on this are rare. Approximations have been used in [2,3] where an analytical form for axial coolant temperature distribution for a pressurized water reactor (PWR) is used. For BWRs, even with reported experimental data, major assumptions were needed to be made about the behavior of axial coolant density distributions in order to create verification models in [4]. However, such assumptions can introduce uncertainties in the calculations because void fraction feedback is not captured [5,6]. e current study seeks to avoid this difficulty by coupling the thermal-hydraulics code TRACE to the T5-DEPL depletion sequence in SCALE for validation of spent fuel composition, thereby avoiding the need for detailed history of void distribution data.
Often, verification studies that use spent fuel compositions as a comparative metric rely on the concept of a ratio of computed to experimentally determined isotopic concentration or code bias (C/E) [2,4,7]. e results from these verification studies can be further used for uncertainty quantification (UQ) [8]. Many methods of uncertainty quantification based on Monte Carlo parametric sampling [3,9,10] exist; data on C/E can be helpful in approximating uncertainty bounds of the isotopic concentrations in spent fuel analysis, similar to what has been done in [11]. An uncertainty analysis of a spent nuclear fuel storage system was performed in [5]. POLARIS was used to calculate the composition of commercial BWR spent fuel for 77 samples.
is study [5] presented and compared BWR spent fuel isotopics for both actinides and fission products. en, the study compared k eff of a spent fuel cask based on codepredicted isotopic compositions with a cask based on experimentally reported isotopic compositions.
Neutronics and thermal-hydraulics coupling involves exchanging data on axial power distribution and coolant characteristics across physics codes [12]. is process can be done with external coupling [13], which is the method used in this study for validation of spent fuel isotopic composition (see Section 3). A variety of neutronics and thermal-hydraulics coupling studies have been conducted previously, which include coupling to subchannel codes (e.g., CTF), coupling to system codes (e.g., TRACE), and coupling to computational fluid dynamics codes (e.g., ANSYS/FLU-ENT). Internal coupling between Serpent 2-SUBCHAN-FLOW was carried out in [14], and the coupling scheme was benchmarked against TRIPOLI4-SUBCHANFLOW and MCNP5-SUBCHANFLOW. In addition, MCS neutronics was coupled to the CTF subchannel code to perform 3D pinby-pin analysis of the BEAVRS benchmark [15]. Applications of coupling involving system codes have been used to simulate BWR accidents such as anticipated transient without scram using Simulate-3K/RELAP5 [16] and BWR instability accidents with TRACE/PARCS [17]. Lastly, nuclear applications that involve coupling to computational fluid dynamics codes were performed in [18,19].
Overall, the present study seeks to develop advanced depletion models which include coupling to the thermalhydraulics system code TRACE to aid in verification and validation of neutronic analysis. e spent fuel isotopic data will be used for validation. 3D neutronic models are developed using SCALE/T5-DEPL in order to capture axial effects in the lattice (e.g., axial power and nonuniform fuel composition). Also, TRACE is coupled to these neutronic models to provide a burnup-dependent axial void distribution. e coupling methodology is applied to two BWR benchmarks with reported isotopic measurements at end of cycle. e remaining sections of this paper are organized as follows. Section 2 presents brief descriptions of the BWR benchmarks used in the validation process. Section 3 presents the SCALE/TRACE coupling methodology. e results of this study and their discussion are presented in Section 4, while the conclusions drawn from this work are given in Section 5.

Fukushima Daini-2 2F2D1 Assembly.
e first benchmark used in analysis is the 2F2D1 BWR assembly from the Fukushima Daini-2 reactor. e design is 8 × 8-4 with 60 fuel rods and one water rod [20]. e top portion of Figure 1 presents the radial layout of the assembly. In this figure, the location of a special rod is marked, labeled "F6," where samples have been taken for isotopic analysis [21]. e F6 rod has an active height of 371 cm with 4.5 wt.% enrichment. Isotopics are measured at three different axial elevations shown in Figure 1(a). Within the dataset, an average void fraction and discharged burnup are reported for each sample. Sample power history is determined by scaling the reported total assembly power to enforce the given burnup to each sample.

Fukushima Daini-2 2F2DN23 Assembly.
e second benchmark used in this study is taken from the same unit as the previous.
is is the 2F2DN23 assembly from the Fukushima Daini-2 unit [22]. Similarly, this is an 8 × 8-2 design but with two smaller water rods. e bottom portion of Figure 1 presents the radial layout with pin-by-pin enrichment of this assembly. Again, the eight locations in which samples have been removed for analysis from a selected rod, the "SF98" rod, are labeled. For each of these samples, time-dependent power, constant and average void history, and discharged burnup are reported. e sample sections taken from this rod are all 0.5 mm in thickness. Sample 1 is taken from the naturally enriched (0.71 wt.%) region of the rod.

Methodology
TRACE [23] is a two-phase flow solver and best-estimate reactor system code developed by U.S. Nuclear Regulatory Commission (NRC) to analyze loss of coolant accident (LOCA), operational transients, and other accident scenarios within light water reactors (LWRs) such as PWR and BWR. TRACE is capable of simulating thermal-hydraulic phenomena in LWRs with high accuracy. It is a 1D (z-direction) code that solves the time-averaged and area-averaged two-phase two-fluid model equations. KENO-V.a is a 3D Monte Carlo neutronic solver in the SCALE code system developed by Oak Ridge National Laboratory [24,25]. KENO-V.a performs 3D Monte Carlo criticality calculations for arbitrary geometries, and it operates in multigroup or continuous energy mode. KENO-V.a is able to determine k eff , neutron lifetime, power and flux densities, and many other important reactor physics quantities.

Multiphysics Coupling.
In this study, the nuclear fuel depletion sequence (T5-DEPL) is coupled to the thermal hydraulic code TRACE. By coupling a neutronics and thermal hydraulics code, there is no longer a need to rely on axial void fraction distributions-as these data can be missing or imprecise. A parametric description of the coupling scheme between these two simulated systems is shown in Figure 2. Internal to the T5-DEPL depletion sequence, KENO-V.a uses Monte Carlo neutron transport to calculate the axial power profile in all nodes, while ORIGEN is used to deplete the fuel in each node based on its respective power level. Externally, TRACE uses the nodal powers to calculate coolant density and temperature and return them to the T5-DEPL sequence for a new neutronics calculation.
e coupling procedure is shown in Figure 3 and can be summarized as follows: (1) e TRACE model is executed first with uniform power in all axial nodes. After completion, coolant density and temperature at each axial node are extracted from the TRACE output.
(2) e coolant density and temperature in all axial nodes are used to update the SCALE/T5-DEPL model. e updated neutronics input is executed and the axial power profile is extracted and normalized.  Figure 1: Pin layout and sample location for 2F2D1 assembly (a) and 2F2DN23 assembly (b) in the Fukushima Daini-2 unit. Sample thickness in SF98 rod is too thin for display ( ∼ 0.5 mm).
Science and Technology of Nuclear Installations 3 (3) e updated axial power distribution is used in a new TRACE calculation, which yields updated coolant temperature and density. (4) Steps 2-3 are repeated until the convergence criterion is met. e convergence criterion uses the average relative difference between two successive iterations in all axial nodes. When this quantity falls below a threshold value (e.g., 6%) for axial power, coolant density, and coolant temperature, the convergence criterion is satisfied and the loop is terminated. (5) Steps 2-4 are repeated until all burnup steps are completed. e next burnup step is initiated upon completion of step 4.

Multiphysics
Modeling. SCALE/T5-DEPL is used to model both the 2F2D1 and 2F2DN23 fuel assemblies in 3D with a reflective boundary condition radially and a vacuum boundary condition axially. e assembly is modeled in a square channel without rounded corners. e Monte Carlo parameters used are 1100 cycles with 20,000 neutrons per cycle, with an initial 100 cycles skipped. ese parameters in KENO-V.a ensure about ∼10 pcm uncertainty in k eff . For 2F2D1, twenty three nodes are used in the axial direction, which include three 3.5 cm thick axial slices to represent the samples (TU101, TU102, and TU106) in the F6 rod. For 2F2DN23, twenty one axial nodes are used, and they include eight 0.5 mm thick axial slices to represent the samples in the SF98 rod. For each sample tested, an irradiation/power history must be given to SCALE/T5-DEPL. Fortunately, there is an option in this code to ensure that the sample power can be specified. e option "power basis" in SCALE/T5-DEPL is used for each measured sample in the rods SF98/F6. To explain, the total power for all materials in the problem will be normalized such that the power in the sample material matches the sample power history specified. e burnup of each sample is measured during experiment and is reported in Figure 1. Matching this burnup value is essential in SCALE/T5-DEPL to achieve accurate validation of fuel isotopics. Fortunately, for the SF98 rod, the sample power is provided directly from the benchmark for each of the eight samples, and it is shown in Figures 4 However, only the total assembly power history is given for rod F6. erefore, the assembly power provided by the benchmark is scaled such that the final burnup matches the sample burnup (i.e., 18.2, 16.1, and 14.0 GWD/MTU in Figure 1). Only the specific power values (kW/kg) are changed, while the total cycle time (418 days for 2F2D1) and time step length are maintained as in the benchmark. e adjusted sample power history is shown in Figure 4(a). e original power history and other neutronic parameters of this benchmark can be found in the SFCOMPO database [21].
TRACE is used to model the two assemblies as a BWR channel (CHAN) component with inlet boundary condition as a FILL component and outlet boundary condition as a BREAK component. e inlet flow rate is set to 17.56 kg/s, while the inlet and outlet coolant temperatures are set to 550 K and 559 K, respectively. e total number and size of the axial nodes in TRACE match those in KENO-V.a. We ensured that mesh sizes in both KENO-V.a and TRACE are appropriate to avoid large numerical uncertainties in the calculations. e extraction of the coolant temperature and density from TRACE is done using the PyPost postprocessor. Notice that since we model the assembly in TRACE as a single component, then the coolant properties calculated by TRACE are generalised for all rods including the SF98/F6 rods in the assembly. To have pin-by-pin resolution of the coolant distribution, subchannel codes (e.g., COBRA-TF) need to be used. Since this problem involves depletion for multiple cycles, burnup-dependent simulation increases the complexity of the coupling process. Based on preliminary tests, we found that TRACE parameters tend to converge after 300-500 seconds following any power changes from KENO-V.a. erefore, to reduce TRACE computational time, only a small time is simulated for coolant properties. Moving on, restart capabilities in SCALE are activated to keep track of isotopic composition. e converged values of material isotopic composition of the previous burnup step are used to initiate the calculations of the next time step. en, on the last burnup step, the converged values of isotopic composition are used for validation. Lastly, it is worth acknowledging that the multiphysics feedback between the two physics is performed only when the power level changes. For 2F2D1, 8 power changes occur, while 2F2DN23 experiences 13 power changes; all are shown in Figure 4. erefore, the convergence of all coupling parameters is ensured 8 times for 2F2D1 and 13 times for 2F2DN23, over the whole depletion time.

Results and Discussion
e results of this study are presented in two sections. First, verification of the coupling process is performed by analyzing the convergence and axial power behavior at the assembly level (i.e., the power basis option is deactivated and all assembly materials are depleted by specifying assembly power). In the second section, validation of the coupling process is performed by analyzing fuel isotopics at the sample level (i.e., the power basis option is activated by specifying sample power). (Assembly-Based). Convergence results at three burnup steps (beginning, middle, and end of irradiation time) are shown in Figure 5 for the 2F2DN23/ SF98 benchmark and in Figure 6 for the 2F2D1/F6 benchmark. e convergence of the three coupling parameters: axial power, coolant density, and coolant temperature, is monitored. e convergence criterion is set to ensure that the average relative difference of all axial nodes falls below a threshold of 6% with at least 3 iterations executed. e threshold value is expressed in relative form regardless of the coupling parameter unit and it balances between accuracy and computational costs. If the average relative difference of all axial nodes falls below 6% on the second iteration, another iteration will still be performed to reach 3 in total. e third iteration is used to identify the convergence, i.e., convergence cannot be ensured with two iterations only. However, if three iterations are not sufficient to converge, additional iterations are executed until all average relative differences fall below 6%.

Neutronic Analysis
For 2F2D1/F6, the results show fast convergence for all parameters as they all converge after 3 iterations (the minimum limit) at the three selected burnup steps. It is clear that the coupling parameters seem to vary little between iterations, implying that the feedback from coolant is not significant for this benchmark. is can be justified by the almost steady power level in Figure 4(a), which causes the void fraction to be almost constant with burnup. erefore, for 2F2D1/F6, the evolution in the axial power profile captured by the SCALE/T5-DEPL 3D model is more significant than the evolution of the coolant properties in TRACE when considering reactor behavior over burnup. On  e notations are as follows: C fuel gives the nodal fuel composition, P is the nodal power distribution, and T cool and ρ cool are temperature and density of the coolant, respectively. P initial is the converged power profile of the previous burnup step. For the first burnup step, P initial is a uniform power density to begin the iterative process.
the other hand, for 2F2DN23/SF98, the results behave differently, as the system may converge within 3, 4, or 5 iterations, implying strong feedback from the coolant. is can be attributed to the higher variability in the power history associated with 2F2DN23/SF98, as can be observed in Figures 4(b)-4(d).
After confirming the convergence behavior of the coupling process, it is good to verify the change of some neutronic parameters with burnup. Figure 7 shows the variation of assembly k ∞ with burnup. For 2F2DN23, we can notice that k ∞ decreases significantly at beginning of cycle due to the xenon effect. Afterward, k ∞ starts to increase due to the depletion of gadolinium poison, until it reaches the reactivity peak, which is close to 10 GWD/MTU. Next, k ∞ continues to decrease toward the end of cycle due to the depletion of U-235.
is is a typical behavior for BWR assemblies that are loaded with gadolinium rods. is behavior can be observed also for the 2F2D1 assembly, except that the simulated cycle length is shorter. We can see in Figure 7(b) that the reactivity peak is very close to 12 GWD/MTU, which is the end of cycle for 2F2D1. It is worth noticing that for some depletion steps in Figure 7(a), k ∞ seems not to converge (e.g., between 15 and 20 GWD/ MTU), since we did not assign a constraint on converging k ∞ as a condition to exit the coupling loop, as described in Figure 2. Doing so will increase the computational costs significantly as k ∞ is a sensitive parameter. Nevertheless, the depletion trend is still clearly seen for this problem in spite of some incompletely converged k ∞ values. e evolution of the axial power for the two sampled rods is plotted in Figure 8 for different burnup steps. For the SF98 rod, the samples span the axial direction. erefore, the locations of the eight samples are used in Figure 8   e power profile for this benchmark did not experience a large change due to the short irradiation time. e bottompeaked shape of the axial power is caused by the void distribution feedback of TRACE. Only small changes in the power peak at the channel bottom occur during burnup.

Isotopic Analysis (Sample-Based).
To begin, all actinides and fission products are tracked, but six major actinides are plotted as a function of burnup at the sample level. e plots for these results are shown in Figure 9 for 2F2D1/F6 (three samples) and Figure 10 for 2F2DN23/SF98 (eight samples). As expected, the U-235 and U-238 concentrations decrease over burnup, while the remaining 4 actinides tend to build up as burnup increases. e differences in burnup for the last   reported concentration in each sample are due to the differences in the sample burnup (see Figure 1). In general, uranium isotopes for all samples tend to change with the same slope. However, it is clear that plutonium isotopic concentrations are more sensitive to depletion conditions, as these isotopes build differently for each sample. e decay time before the rod measurement ( ∼ 1500 days after rod discharge) causes a sharp decrease in Pu-241 concentration at end of life. is is also the case for the small increase in Pu-239 at end of cycle, due to the decay of other fission products into Pu-239. In addition, Figure 10 shows that sample 1 discharge burnup is much smaller than other samples, due to its location in the natural uranium region (0.71 wt.%). e validation results of fuel isotopics are demonstrated in Figure 11 for selected isotopes in the two benchmarks. In addition, numerical values of the C/E ratio are listed in Table 1 for 2F2DN23/SF98 and in Table 2 for 2F2D1/F6. For completeness, the experimentally determined isotopic data are reported in Tables 3 and 4 of Appendix A, which can be used in conjunction with Tables 1 and 2 to obtain the calculated results. In the scatter plot of experimental versus calculated, the points that are closer to the diagonal line indicate better agreement with experiment. In this plot, the uranium isotopes as well as Nd-148 tend to have good agreement with data as can be inferred from their good C/E values. e good agreement for uranium isotopes can be attributed to their large content in the fuel and their low measurement error. Also, it was seen early that these isotopes are relatively insensitive to depletion conditions. e plutonium isotopes tend to have fair agreement with data, as their relative error can be as low as 2%, but as high as 64%.
In terms of uncertainty sources associated with this work, we believe uncertainty comes from three major sources: first, SCALE/T5-DEPL and the nuclear data libraries. Uncertainty in nuclear data, especially for minor actinides and fission products, can be a large factor particularly because the 56 multigroup energy library was used, which could be less accurate as compared to using continuous energy cross section library. e 56-group energy library was used because SCALE/T5-DEPL depletion calculations tend to be very slow with continuous energy cross sections. is coupling process requires multiple SCALE/ T5-DEPL runs, so using continuous cross section libraries becomes impractical in terms of computational costs. e second uncertainty source is associated with the TRACE model, which is directly connected to the third source about lack of reported data. As this is a neutronics benchmark, very little information is provided to aid in building an accurate thermal-hydraulics model representing the fuel channel.
erefore, some assumptions about boundary conditions based on commercial BWRs (inlet flow rate, outlet temperature, friction coefficient, etc.) had to be made. ese assumptions would affect the quality of coolant density and temperature predictions by TRACE. Additional missing pieces associated with neutronics include configuration of the assemblies that surround 2F2D1/2F2DN23 in the core, control blade history, SF98/F6 rod burnup, and fine-scale power history compared to the approximated profile provided by the benchmark. ese factors are very important for realistic 3D modeling. Despite the previous challenges, satisfactory results are achieved at some axial elevations as indicated by the last row of Tables 1 and 2. e average relative difference for all isotopes can be less than 15% for some samples but can also be higher than that for other samples. It is worth mentioning that the current validation results are similar to previous efforts such as [4,5,8]. Finally, since lack of documentation was a major challenge for this work, future efforts will focus on incorporating machine learning and data science to aid in identifying the missing input parameters such that they can yield the desired output. Techniques such as Gaussian process optimization and deep reinforcement learning can be efficiently used for this purpose. ese efforts can help in completing the current benchmarks to achieve accurate modeling and simulation.

Conclusions
Multiphysics coupling of SCALE/TRACE is used in this work for neutronic analysis and spent fuel validation of BWR spent    e coupling is performed externally where SCALE/T5-DEPL module sends axial power profile to TRACE, which returns coolant density and temperature. is process is repeated until convergence of all quantities at each time step. e methodology is applied to two BWR assemblies, discharged from the Fukushima Daini-2 unit. e neutronic analysis on assembly level shows that the coupling process is verified, where the 2F2D1/F6 benchmark seems to have insignificant void feedback. Validation of the coupling process is done through spent fuel isotopics data, which demonstrated good results for uranium isotopes and satisfactory results for other isotopes, given the significant lack of documented data faced in this work. Future work will include applying machine learning methods to help in identifying the missing data such that more accurate modeling and results can be achieved.