CTF-PARCS Core Multi-Physics Computational Framework for Efﬁcient LWR Steady-State, Depletion and Transient Uncertainty Quantiﬁcation

: Best Estimate Plus Uncertainty (BEPU) approaches for nuclear reactor applications have been extensively developed in recent years. The challenge for BEPU approaches is to achieve multiphysics modeling with an acceptable computational cost while preserving a reasonable ﬁdelity of the physics modeled. In this work, we present the core multi-physics computational framework developed for the efﬁcient computation of uncertainties in Light Water Reactor (LWR) simulations. The subchannel thermal-hydraulic code CTF and the nodal expansion neutronic code PARCS are coupled for the multi-physics modeling (CTF-PARCS). The computational framework is discussed in detail from the Polaris lattice calculations up to the CTF-PARCS coupling approaches. Sampler is used to perturb the multi-group microscopic cross-sections, ﬁssion yields and manufacturing parameters, while Dakota is used to sample the CTF input parameters and the boundary conditions. Python scripts were developed to automatize and modularize both pre- and post-processing. The current state of the framework allows the consistent perturbation of inputs across neutronics and thermal-hydraulics modeling. Improvements to the standard thermal-hydraulics modeling for such coupling approaches have been implemented in CTF to allow the usage of 3D burnup distribution, calculation of the radial power and the burnup proﬁle, and the usage of Santamarina effective Doppler temperature. The uncertainty quantiﬁcation approach allows the treatment of both scalar and functional quantities and can estimate correlation between the multi-physics outputs of interest and up to the originally perturbed microscopic cross-sections and yields. The computational framework is applied to three exercises of the LWR Uncertainty Analysis in Modeling Phase III benchmark. The exercises cover steady-state, depletion and transient calculations. The results show that the maximum fuel centerline temperature across all exercises is 2474K with 1.7% uncertainty and that the most correlated inputs are the 238 U inelastic and elastic cross-sections above 1 MeV.


Introduction
Best Estimate Plus Uncertainty (BEPU) [1] approaches have received a lot of interest during the last few years in the context of nuclear reactor safety analysis. Efforts are being made to improve the multi-physics and multi-scale computational modeling of nuclear reactors. Various coupling approaches have been used to establish robust computational models under different conditions from steady state to transient with reasonable accuracy and computational efficiency. Typically, such modeling approaches restrain the multi-physics analysis to a coupling between thermal-hydraulics and neutronics. The fuel performance is simplified to capture the main feedbacks through the thermal solver inside the thermal-hydraulics modeling. This is performed because fuel performance is computationally very expensive due to the complexity of the multi-scale phenomena involved. Various computational frameworks have been developed and established around the world to develop such multi-physics calculations. In the U.S., VERA [2], NEAMS [3] and MOOSE [4] are the main examples that can have an even higher multi-physics fidelity. In France, CEA has developed its own platform using SALOME [5], while EPFL and PSI in Switzerland have explored the OpenFOAM capabilities to establish their platform [6]. In Finland, VTT has recently developed a new multi-physics platform called Kraken [7] and in Korea, UNIST has developed MPCORE [8], a multi-physics coupling that can even model the core at a pin-by-pin scale. In this work, we present the core multi-physics computational framework for efficient Light Water Reactor (LWR) steady state, depletion and transient uncertainty quantification established at North Carolina State University (NCSU) [9]. In this framework, the subchannel thermal-hydraulic code CTF [10] and the nodal expansion neutronic code PARCS [11] are coupled at the reactor core level. This CTF-PARCS coupling leverages the ongoing developments of CTF at NCSU.
A wide variety of methods are available for coupling the different physics. A first classification of these coupling methods is usually made by dividing them into tight and loose coupling [12]. In tight coupling, the different physics are embedded in a unique code that solves all the constituent equations. On the other hand, in the loose coupling, the different physics are solved separately using an external iteration scheme to track the convergence of the coupled solution. Furthermore, loose coupling methods can be divided into internal and external coupling schemes [13]. The internal coupled codes are compiled in the same project, and the communication and data exchanges take place by sharing their internal memory. On the other hand, the externally coupled codes remain isolated in their own compilation projects, and the communication and synchronization relies on a coupling interface based on any of the available inter-process communication protocols. Different solution strategies are available for solving the coupling of the different physics and their effectiveness depends on the above-mentioned coupling methods. For example, for tightly coupled codes, the implementation of the implicit Jacobian-free Newton-Krylov (JFNK) method [12] is desirable because of its high convergence rate. Intrusive and extended code modifications are needed for the implementation of the JFNK method. For this reason, although conceptually it can be applied in loosely coupled codes, it is usually reserved for tight coupling. For loose coupling methods, the Picard iteration method [14] is widely used due to its simplicity. This method has the associated problem of poor convergence rate, but this can be remedied by applying acceleration methods such as the Alternating Nonlinear method, the Residual Balance, or the more convenient Adaptive Residual Balance [15]. The Semi-Implicit Operator Splitting [12] is another well-developed and optimized solution strategy to consider in loosely coupled codes. Applications of these methods and solutions strategies to the neutron kinetics and thermal-hydrualics coupling can be found in these references [16][17][18][19]. From the presented methods and solution strategies the loose coupling and the Picard iteration have been selected for the CTF-PARCS coupling. No acceleration methods are currently implemented apart from the relaxation of the relative tolerance of CTF.
The developed CTF-PARCS coupling is embedded in an uncertainty quantification framework that allows efficient and robust statistical analysis of various input and output quantities of interest. This allows the framework to be applicable to BEPU approaches that alleviate some of the conservatism involved in safety studies. In this context, the OECD/NEA Light Water Reactor Uncertainty Analysis in Modeling (LWR-UAM) Benchmark [20] was initiated to facilitate the development of uncertainty quantification methodologies and propagate consistently uncertainties across multi-scale and multi-physics coupled calculations. The application of BEPU approaches to transient core calculations is challenging due to the increased computational cost. In [21], the CEA multi-physics platform was used to perform such studies, while in [22], a coupling between PARCS and system thermal-hydraulics code TRACE [23] was used for I-2a, I-2b and I-2c Pressurized Water Reactor (PWR) exercises of the LWR-UAM Phase III [24]. In this work, we perform a similar uncertainty quantification study using the developed CTF-PARCS computational framework for the same exercises of the LWR-UAM Phase III. These exercises involve the modeling of the TMI-1 core under various conditions. For Exercise I-2a, we model the steady-state Hot Full Power (HFP) at the Beginning of Cycle (BOC) and End of Cycle (EOC). For Exercise I-2b, we perform the depletion as a multi-state coupled calculation from BOC to EOC. For Exercise I-2c, we model a rod ejection accident (REA) initiated from HFP conditions at both BOC and EOC.
This article is structured in five sections; Section 1 was the above presented introduction. Section 2 discusses in detail the developed multi-physics computational framework. The different codes that are being coupled are first presented followed by a description of the Python tool used to streamline the transition from lattice to core neutronics. The implemented CTF-PARCS coupling approaches for steady-state, depletion and transient calculations are discussed, and the uncertainty quantification approach is presented. Section 3 describes the LWR-UAM Phase III uncertainty quantification exercises selected to apply the developed computational framework. Section 4 presents and analyzes the obtained results for all the studied exercises. Finally, in Section 5, the general conclusions from this work and future intended developments are discussed.

Computational Framework
This section describes the computational framework used to simulate the proposed steady-state, depletion and transient calculations. It focuses on the capabilities and roles of the standalone codes involved in the framework, as well as the interaction of the codes embedded in the framework to perform the calculations and the adopted uncertainty quantification approach.

Simulation Tools
This subsection contains a short description about each one of the simulation tools involved in the developed computational framework. This includes the main capabilities and the most important features of these codes used in our current activities.

Polaris
Polaris [25] is a two-dimensional (2D) deterministic neutron transport code from the SCALE 6.2.4 suite [26]. Polaris is used to produce multi-history, multi-branch cross-section libraries at different burnup. The cross-sections are condensed from one of SCALE 6.2.4's predefined fine-group structure into a coarse-group structure. The generation of the crosssections follows a traditional approach, where each fuel assembly representative of the LWR core of interest is modeled individually in Polaris to generate a set of homogenized fewgroup cross-sections in the SCALE format. Reflective boundary conditions and a critical spectrum are applied in these lattice models. Polaris can also generate homogenized fewgroup cross-section libraries for the reflector regions of the core. The models implemented to generate the nominal set of cross-sections are described in Section 3.2.
The generation of the two-group cross-sections with Polaris results in a set of SCALEformatted libraries (hereafter referred to as t16 format). Each library contains cross-sections representative of one fuel or reflector assembly constitutive of the core. The Polaris libraries contain information at different burnup and branch points. The burnup points are required in the cross-section modeling to account for the burning, breeding, decaying and activation of actinides and fission products found in the materials depleted in the Polaris models.
Branches are used to model instantaneous state variations (temperature, density, soluble boron concentration, etc.) during operations and transients. Histories account for the longterm effects of the neutron spectrum on the burned materials. The historical parameter is modeled separately in Polaris, which means that modeling H different historical points for a given fuel assembly requires H Polaris inputs. Overall, a core containing A F types of fuel assemblies and A R types of reflector assemblies with H F fuel assembly histories and H R reflector assembly histories requires T = A F × H F + A R × H R Polaris inputs and produces T libraries in t16 format.

Sampler
Sampler is a module of SCALE 6.2.4 [26] used in association with Polaris to generate perturbed cross-section libraries. Sampler is used in the framework to perturb the multi-group cross-sections, the fission yields of fissionable actinides, and manufacturing parameters of the lattice cell. Perturbing consists of modifying simultaneously the uncertain inputs (cross-section, yields, manufacturing parameters) in the assembly models to obtain sets of cross-sections libraries representative of the input perturbations. The perturbed cross-section libraries are used in the CTF-PARCS simulations to assess the propagation of uncertainties (cross-sections, yields, manufacturing parameters) into key core outputs of interest.
In Polaris-Sampler, the range and nature of the sampling distribution for the manufacturing perturbations are user-defined and are described for our application in Section 3.2. The number of degrees of freedom is limited to the number of samplings for the crosssections and fission yields perturbations, because perturbation factors have been precomputed in the SCALE libraries. The perturbation factors were generated in SCALE from the covariance data available in ENDF/B-VII.1 [27] along with low-fidelity covariance data supplemented for minor fission products and actinides.

GenPMAXS
The Polaris generated libraries in t16 format are not readable by the core simulator, PARCS (Section 2.1.4), and need to be converted into a PMAXS format using the GenPMAXS sequence [28].
The GenPMAXS program combines into a PMAXS library the history, branch and burnup points of a given assembly. The construction of the GenPAMXS inputs relies on:

•
The configuration of the fuel or reflector assembly defined by the assembly specifications, such as the pin pitch, number of pin per assembly. Those specifications are required in the PMAXS for pin power reconstructions. • The cross-section parameterization (branch, history, burnup). GenPMAXS also requires the definition of a reference history and branch. • Pre-defined user inputs. For example, the name given to the PMAXS file, boolean flags such as inclusion of assembly and axial discontinuity factors.

PARCS
PARCS is a 3D reactor core neutron kinetics simulator code that solves the steady-state and time-dependent multi-group neutron diffusion and low-order transport equations (SP3) in Cartesian and hexagonal geometries [11]. Two additional important features of PARCS are the capability of reading burnup-dependent cross-sections in PMAXS format and making pin power reconstruction when group-dependent flux form functions are provided within the cross-sections library. PARCS exists as a standalone code and as a coupled version to TRACE. The TRACE-coupled internal version is bundled directly into the TRACE distribution and has been used in this project as it facilitates the multi-physics multi-scale coupling between TRACE, CTF, and PARCS.

CTF
CTF is the shortened name given to the version of COBRA-TF being developed and improved by the Reactor Dynamics and Fuel Modeling Group (RDFMG) at North Carolina State University (NCSU) [10]. CTF is a subchannel thermal-hydraulic simulation code designed for Light Water Reactor (LWR) vessel and core analysis. It uses a two-fluid three-field modeling approach for the two-phase flow. In the last decade, CTF has been extensively developed and validated for PWR, BWR, VVER, Small Modular Reactor (SMR), Fast Breeder Reactor (FBR), and research reactor applications. CTF is a state-of-the-art subchannel code for reactor thermal-hydraulics bundle and core analysis and is a part of the U.S. DOE CASL [29] and E.C. NURESAFE [30] projects.

Dakota
Dakota [31] is a general toolkit to perform uncertainty quantification, optimization and sensitivity analysis for any code that can be treated as a black box. In this work, the code is the multi-physics coupling between CTF and PARCS, and Dakota is used for performing the uncertainty propagation and more specifically the sampling of CTF input parameters and boundary conditions and the parallel distribution and execution of the different calculations. Custom Python scripts are used for providing an interface between the code and Dakota and for post-processing and analyzing the obtained results.

Polaris Pre-Processing
In general, 50 to 100 Polaris inputs are required to construct a full-core LWR model in PARCS to account for the core history of each assembly type in nominal conditions. The construction of multi-history inputs of different assemblies requires the modification of the following parameters in the Polaris inputs: • The pin layout. • The material definitions. • The reference state. • The core plate thickness or shroud thickness, if any, for reflector assemblies.
The historical, branch and burnup points can remain identical for each fuel assembly. The same statement applies for reflector assemblies, although the reflector parameterization does not have to align with the fuel assemblies' parameterization. After establishing the parameterization of the nominal configuration of each assembly, Sampler creates N perturbed cases per assembly history. The execution of Polaris/Sampler produces N t16 libraries. These libraries must be converted into the PMAXS format; therefore, the GenPMAXS inputs must be implemented for each assembly required by PARCS for each perturbation computed used for the uncertainties analysis. The manual implementation of the Polaris and GenPMAXS inputs is error-prone, redundant and time-consuming; therefore, the generation of these inputs as well as the code executions were streamlined into a unique data management sequence. The sequence accelerates the generation of the PMAXS libraries, which enhances the consistency and robustness of the libraries. It relies on two Polaris templates files: one for the fuel assemblies and one for the reflector assemblies. Both templates are formatted in the Sampler formulation to avail the generation of the perturbed assemblies. One of the modules of the sequence contains the properties of the assemblies, such as the pin layouts of fuel assemblies and the historical points chosen for the cross-section parameterization. The branch points and the burnup points are hardcoded into the templates because they are uniform across all the fuel and reflector assembly models. Optionally, the dimension and composition of the pins can be implemented in the sequence with corresponding uncertainties, if manufacturing uncertainties are included in the model. The construction of the Polaris assemblies is illustrated in Figure 1.
After execution of the sequence, the creation of the Polaris inputs is separated into two phases: (a) the nominal values are used to produce the nominal Polaris assemblies, and (b) the Sampler inputs are built; therefore, the Sampler aliases are preserved from the template, and a distribution block is added at the end of the Polaris inputs using the nominal values and the associated standard deviations. After the construction of the nominal and perturbed inputs, the nominal cases are executed, which are followed by the perturbed cases. Typically, 100 libraries per assembly and per history leads to about 10,000 perturbed Polaris inputs/outputs. If the computational power available is sufficient to run the entire set of nominal and perturbed Polaris simulations, the sequence will submit the jobs and terminate when the t16 libraries exist. Otherwise, the sequence can be stopped while the Polaris simulations are running and be resumed at a later time. The sequence is set to submit the simulations on a Linux High-Performance Computing with a Slurm workload; the portability of the sequence depends on the Operating System and the workload. Additionally, the perturbation of manufacturing parameters involves the perturbation of the pin dimensions or reflector dimensions in the assembly, which can induce ray spacing errors in SCALE 6.2.4. It was found that~8% of the Polaris inputs failed due to ray spacing error with the default spatial discretization; this number tends to increase if the spatial discretization is refined. To handle robustly the situations of (a) resuming the submission of the simulation on the HPC or (b) fixing ray spacing issues, the sequence reads the existing Polaris outputs to re-run or re-index cases. Therefore, a set of M perturbed inputs are generated, while N perturbed are actually run, to have E inputs cases ready to substitute the failed input configurations (see Equation (1)).
At a given perturbation n, 1 ≤ n ≤ N, the sequence restarts are handled as follows: • If a SCALE execution file exists (exec file), the job is currently running, and the sequence moves on to the next perturbation, n → n + 1. • If no t16 file and no Polaris output (.out) exist, the perturbation case has not been run yet: the sequence runs the simulation normally and moves on to the next perturbation, n → n + 1.

•
If the t16 file and Polaris output file both exist, the number of state points computed in the t16 file (S t16 ) is cross-compared to the expected number of state points predicted by the Polaris output (S out ).
-If S t16 = S out , the perturbation case has already run normally, and the sequence moves on to the next perturbations, n → n + 1.

-
Otherwise, the sequence raises an error, and no diagnostic is predicted.
• If the Polaris output exists, but the t16 does not, the Polaris output is read by the sequence to find out if the simulation stopped due to a ray spacing error.
-If a ray spacing error occurred, another perturbation is drawn from the E supplemental inputs. The supplemental case indexed m is re-indexed into n, the case is re-run, and then the sequence moves on to the next perturbations, n → n + 1.

-
Otherwise, a job crash is assumed (exceeded wall time, reboot, etc.) and the perturbation n is re-run without modification and then moves onto the next perturbation, n → n + 1. A summary of the execution and verification phase is provided in Figure 2. Once all the nominal and perturbed t16 libraries are generated, the sequence advances to the PMAXS generation phase. The GenPMAXS files are constructed immediately from the information contained in the Polaris input files without using templates. The type of reactor, the name of the assembly, the number of history points, the number of energy groups, the types of state parameters (boron concentration, fuel temperature etc.), the values of the historical points, the number of historical points, the number of branches, the values of the branch points, the number of burnup points, the values of the burnup points and the names of the t16 libraries are extracted from the Polaris input. The sequence selects a reference history point based on the type of reactor by calculating the Euclidian distance between each historical points to the predefined reference conditions of the reactor. The reactor conditions involve the fuel temperature (T F ), the moderator temperature (T M ), the moderator density (ρ M ) and the boron concentration (C B ). For example, the PWR reference conditions are defined in the sequence with T F = 900 K, T M = 573 K, ρ M = 0.726 g · cm −3 , C B = 900 ppm. This approach is repeated to determine the reference branch point. The same approach also operates for the reflector, with the exception that PMAXS libraries are generated for corner reflectors in addition to the regular reflector libraries. Corner reflector libraries are computed from regular reflector libraries, using an additional parameter r 2D in GenPMAXS computed from the assembly pitch P FA and the shroud thickness d, as shown in Equation (2).
After the generation of the GenPMAXS inputs, the sequence executes automatically the GenPMAXS inputs to generate the PMAXS data for the both the nominal case and the perturbed case.

CTF-PARCS Coupling
The subchannel thermal-hydraulic code CTF and the core neutron kinetics simulator PARCS have been recently coupled at NCSU. The capabilities of this coupled code are being expanded, and the code can be currently used as a three-dimensional multi-physics core simulator for steady state, depletion and transient. This core simulator, referred to as CTF-PARCS in this manuscript, will be used in the multi-physics simulation framework for uncertainty quantification proposed here. A flexible Message Passing Interface (MPI) communication protocol based on a server-client coupling algorithm has been developed as a kernel of our multi-physics multi-scale simulation platform. CTF-PARCS is the core simulator in the general simulation platform developed at NCSU and presented in Figure 3. This coupling algorithm facilitates the implementation of loose coupling methods between different simulation tools without the necessity of extended modification in their sources and compilation projects [32]. In the developed CTF-PARCS multi-physics core simulator, CTF acts as the server while PARCS is the client, as indicated in Figure 4. The coupled simulation is activated by using a command line argument in CTF, being as a server able to launch additional MPI processes acting as clients, as will be PARCS in this case. A coupling interface input file (input-name.ci) handles the information about the different clients launched by the server. The neutron kinetics coupling class developed in CTF contains protocols to control the run-time advance of PARCS, such as initialization, steady-state iteration, depletion step, transient time step, convergence, finalization, and error. The synchronization is managed through tags communicated using the MPI standard functions MPI_SEND and MPI_RECV. The feedback parameters (included in Figure 5a) are exchanged during each global coupled code iteration using the MPI communication PORT opened between both codes. Mapping and auto mapping procedures have been developed to provide the information about the different physical domains, spatial discretization, and spatial mesh overlays. It must be noted that CTF is a MPI parallel code, and different approaches can be taken regarding the subdomain decomposition. Currently, the MPI process with rank 0 in CTF is making the full domain communication in the coupling algorithm following in fact a global mapping scheme provided as input to CTF.  To analyze the global convergence of the core simulator, the three developed simulation modes, steady state, transient, and depletion, must be analyzed separately.

•
In steady-state mode, a simple Picard iteration scheme [14] has been adopted. This method implies the construction of an outer convergence loop to check for global convergence of the coupled code and as many local convergence loops as physics are being coupled together. These nested loops represent the global coupled convergence. In this method, the convergence of the thermal-hydraulic code has to be fully met before providing feedback to the neutronic code. The neutronic code then uses the thermal-hydraulic feedback to resolve the neutron flux distribution until its own convergence criteria are met. The global coupled code convergence is tracked by the residuals of the feedback parameters using the L 2 and L ∞ norms of each parameter until the defined tolerance is reached. The L 2 and L ∞ norms are defined in Equations (3) and (4), respectively, where i is the cell number, nc is the total number of cells, j is the index of the residual of a given feedback parameter, and nr is the total number of residuals. R i,j is the residual of the jth feedback parameter in the ith cell for the current iteration, which is computed as R i,j = R n i,j − R n−1 i,j , among which n is the current iteration number. The coupled code steady-state convergence behavior is illustrated in Figure 6.
• For transient simulations, a temporally explicit coupling scheme has been adopted. Therefore, the neutronic and thermal-hydraulic conservation equations are not solved altogether, but each code solves its own system of transient equations using the most limiting time step size of both codes. This explicit client/server coupling scheme is presented in Figure 5b. The temporally explicit methods are usually adopted when coupling thermal-hydraulics with neutron kinetics because the numerical stability is always preserved. Although it is widely adopted, the temporally explicit coupled methods have two main inconveniences. The first issue is the numerical diffusion, as each of the coupled physics relies on the solution of the other physics at the previous time step. The second inconvenience is the limitation of the time step size, as all the coupled codes must adopt the smaller time step size between all the involved physics. • In depletion mode, a multi-state simulation procedure is initialized in the core simulator. The number of states of this multi-state mode is defined by the number of depletion steps set in PARCS input. At the beginning of each state, the 3D burnup distribution is sent from PARCS to CTF to initialize the burnup-dependent models and material properties in the fuel rod objects. After this initialization, the Picard iteration global convergence scheme is used in the same way as in steady-state simulations. It must be noted that the boron concentration is not feedbacked from CTF to PARCS in depletion mode, as the critical boron is the parameter varied by the latter to find the criticality of the core. This multi-state iterative procedure is represented in the scheme of Figure 7.
The global convergence of the developed core simulator in depletion mode can be analyzed by looking at the L ∞ norms presented in Figure 8. In Figure 8, it can be observed that the global convergence is longer in the first iteration because the initial conditions are far from the real solution. All the remaining states converge with a similar convergence ratio. Each peak on the residual represents the beginning of each state. The peaks are large due to the change in the 3D burnup distribution in the core.
Some additional capabilities have been implemented in CTF during the development of CTF-PARCS coupling for improving the prediction of the effective fuel temperature: • The first development is the possibility of defining customized ways to compute the Doppler fuel temperature. This effective fuel temperature can be defined as a linear combination of the pellet radial average temperature (T avg ), the fuel surface temperature (T surf ), and the fuel centerline temperature (T cen ), as shown in Equation (5).
The parameters A , B, and C are defined by the user in the coupling interface input file.
This allows the user complete flexibility on the definition of the effective temperature. Typical effective fuel temperature formulas found in the literature that can be reproduced with this development are the Rowlands [33], as shown in Equation (6), and the Santamarina [34], as shown in Equation (7).
• The second added feature is the initialization of the pellet power and exposure radial profiles in CTF. This capability has been implemented by linking an external 1D depletion module to CTF. This external depletion module can be called from the fuel rod class of CTF to make a depletion and obtain the radial power and exposure profiles for each nuclear rod and axial node of the model.
In this project, the effective fuel temperature is calculated using the Santamarina formula [34] of Equation (7), and the radial power and exposure profiles within fuel pellets are calculated using the above-mentioned capability. For the latter, we have selected the TUBRNP model developed in TRANSURANUS code [35] from the publicly available model to compute the radial dependence of the power and burnup in the fuel pellets. This model determines the power profile based on the fuel rod geometry and the initial fissile material concentration. Single group thermal flux diffusion theory is applied to calculate the concentration of 235 U, 238 U, 239 Pu, 240 Pu, 241 Pu, and 242 Pu isotopes in each pellet ring. The microscopic cross-sections are spatially constant. Bessel functions are used to solve the 1D single group diffusion equation, obtaining the flux shape within the pellet as a function of the isotopic composition of each ring. For further details, the reference [35] can be consulted. This model has been coded and compiled as an external library and is used to initialize the fuel rod power and burnup profile at the beginning of the CTF-PARCS steady state and before every one of the depletion steps. PARCS sends the 3D burnup distribution to CTF at the beginning of the simulation and at the end of each depletion step.

Uncertainty Quantification
The previously discussed multi-physics computational model is embedded into an uncertainty quantification framework, as illustrated in Figure 9. The framework combines different custom Python scripts with Dakota functionalities to perform a consistent and efficient uncertainty quantification study. The consistency of the framework requires applying all the multi-physics inputs perturbations to all the relevant physical domains. For example, if the fuel rod radius is perturbed, this perturbation needs to be applied to Polaris but also to be considered in the fuel rod geometry, in the calculations of the flow areas, and in the wetted perimeters in CTF. The custom-made Python scripts allow an easy and straightforward application of various multi-physics inputs. Once the perturbations are applied, Dakota distributes the calculations in parallel and executes them until all the outputs have been obtained. The sampled inputs and outputs are used for the statistical analysis, where different statistical quantities such as mean, standard deviations and correlations are calculated.
The work emphasizes estimating the correlations between multi-physics outputs and the 56-group microscopic cross-sections available in SCALE and used by Sampler and Polaris. In state-of-the-art applications, the two-group macroscopic cross-sections are one of the figures of merit correlated to the downstream outputs. The objective here is to (a) retrieve a finer level of information than two-group macroscopic data, and (b) correlate the uncertain inputs from a more upstream level (i.e., the 56-group cross-sections before the lattice calculations) to the multi-physics outputs. The fundamental difficulty in such conventional two-step multi-physics calculations is that the computational cost does not allow any detailed sensitivity analysis to determine the important inputs. Recent studies [22] addressed this issue by perturbing independent separate groups of inputs and evaluating the standard deviations obtained for each group with the standard deviation obtained when all the inputs are perturbed altogether. The issue with such approaches is that (a) it requires larger (but feasible) number of calculations, and (b) the grouping of inputs is done at a very coarse level. For example, the inputs are separated by physics domains (e.g., neutronics, thermal-hydraulics). The large number of inputs that are highly correlated constitutes the major challenge for going to a finer level of sensitivity analysis. This is specifically true for neutronics due to the large number of microscopic cross-sections and yields stored in the continuous energy and multi-group libraries. The other sources of input uncertainties (e.g., manufacturing parameters, boundary conditions) can be assumed independent with reasonable assumptions and do not raise problems for the statistical analysis. In the developed computational framework, we exploit the existence in SCALE of pre-calculated perturbation factors for all the 56-group microscopic cross-sections and yields. The PALEALE utility of the SCALE suite was used to pre-process and extract all these perturbation factors. This pre-processing needs to be done only once, since the perturbations are constant and they apply to every uncertainty quantification study performed with Sampler. The obtained microscopic cross-sections and yields perturbations can then be used to estimate their correlations with the multi-physics outputs of interest. Due to the large quantity of data, only a selection of the most important isotopes is analyzed: 235 U, 238 U, 239 Pu, 1 H, and 16 O. For these isotopes, the reactions listed in Table 1    Although the yields have been studied, their impact in this work is found negligible for the analyzed cases and thus is not discussed. The Spearman rank correlation coefficients are used for estimating the correlations. Spearman coefficients are computed based on the ranking of the variables and thus capture monotonic trends. If we assume a sample of size N where R x and R y are the vectors of ranks for each sample of two quantities x and y, respectively, then the Spearman coefficient r xy can be computed through Equation (8). The Spearman coefficients are preferred over the typical Pearson correlation that captures only linear relationships. In the future, more advanced methods such as the Hilbert-Schmidt Independence Criterion (HSIC) indices could be considered [36].
The developed statistical framework allows the calculations of correlations between multi-physics outputs and the 56-group SCALE microscopic cross-sections. Before presenting the different studies and their results, it is important to discuss three aspects relevant to such an approach: • The first aspect is the statistical noise. Due to the limited number of samples, even independent quantities can be correlated up to some degree. For this reason, a cutoff threshold needs to be determined in order to screen out all the correlations that are within the statistical noise. To calculate the cutoff threshold, a very simple brute force approach is performed. For a predefined number of iterations N o = 10 6 , in each iteration, two standard normal variables are sampled independently with the intended number of samples N. The correlation between these two variables is then estimated and stored. This results in a vector of correlations of size N o that represents a distribution of empirically estimated correlations of independent variables and thus can be used to derive a cutoff based on some statistical measure. The 99.9% percentile is considered in this work as the cutoff. As we are going to see later, N = 100 samples were used, and the calculated cutoff threshold was found to be 0.35. This means that any correlation that, in absolute value, is larger than 0.35 is considered as statistically significant. The cross-sections that show at least one group to be statistically significant are analyzed at all their energy groups, as shown in the uncertainty quantification results section, in order to better understand the energy dependence of the estimated correlations. • The second aspect is the fact that if a very large number of inputs contribute approximately equally to the output, then it will be difficult to even distinguish them from the noise. This can be remedied by increasing the number of samples, since the cutoff threshold will decrease. Additionally, in most of the applications, few inputs contribute significantly more than the others, which facilitates their identification. The latter of course cannot be known in advance but only after performing the multi-physics calculations for the specific application. • The third and most important aspect is that the estimated correlations do not imply causation; in other terms, if an input shows strong correlation with an output, it does not mean that it is actually responsible for the output's uncertainty. For this reason, a subjective analysis of the correlation matrices together with the physics understanding and experience can help deduce whether an input is actually important or it has high correlation, because it is correlated to another input, which is the one that is actually important. In any case, even if there is no objective way for such a simple uncertainty quantification approach to assign sensitivities in the multi-physics context, the resulting correlations are still very valuable and can inform future experiments to reduce the corresponding input uncertainties or correlations with other inputs. Finally, an ultimate test for this approach would be to compare the identified microscopic cross-sections for a steady-state standalone neutronics study with the ones identified using a perturbation-based approach. The validation of this proposed approach is left for future work.

Case Studies
The exercises I-2a, I-2b and I-2c of LWR-UAM Phase III [24] consist of various transient and steady-state calculations for the TMI-1 core. The TMI-1 core has a nominal power of 2772 MW th and consists of 177 UO 2 fuel assemblies with 11 different designs. The different designs consist of variations in UO 2 enrichment, number of gadolinium rods and number of discrete burnable poisons.
Exercise I-2a includes two steady-state calculations for BOC at HFP and Hot Zero Power (HZP) conditions. The main goal of this work is to apply and test the robustness of the multi-physics framework, so we focus on the HFP condition at both BOC and EOC. The 3D burnup distribution at BOC and EOC is provided in the benchmark specifications, and its axially averaged radial distribution is shown in Figure 10. The average core burnup is 18.08 GWD/MTU at BOC and 42.06 GWD/MTU at EOC. Exercise I-2b consists of a cycle depletion calculation from BOC to EOC. The cycle length is 664 Effective Full Power Days (EFPD). Exercise I-2c includes a series of REA from different initial conditions. The control rod with the highest worth, which is fully inserted into the core, is ejected at a constant rate of 2380 cm/s. No SCRAM is simulated to allow for a better understanding of the feedback effects. The location of the ejected control rod is the D12 in Figure 10. During the transient, the boron concentration and the position of the other control rods are assumed to be constant. Four different initial states are considered for the REA in the LWR-UAM benchmark: HZP BOC, HFP BOC, HZP EOC and HFP EOC. In this work, only the HFP BOC and EOC are modeled, because they lead to the most limiting condition in terms of the maximum fuel temperature and fuel stored enthalpy. The inserted reactivity in all cases is relatively small, and the HZP cases do not lead to large temperature variations. The input uncertainty quantification for this work was performed based on the LWR-UAM Phase II [37] and Phase III [24] specifications and recommendations. It is important to mention that for the fuel and cladding thermal conductivities, the reduced uncertainty recommended in LWR-UAM Phase III and obtained after a Bayesian calibration is used in this work. A total of N = 100 macroscopic cross-sections samples (Section 3.2) for each unique assembly and reflector were generated with Sampler (Section 2.1.2). The Wilk's formula was used to calculate that 93 samples were required to obtain a 95% tolerance interval with 95% confidence. A maximum of 1000 samples in SCALE can currently be generated from the SCALE libraries. Due to the computational cost for such uncertainty analysis, a sample size of N = 100 was deemed sufficient to test the computational framework. This sampling size choice was also justified in [38], where some of the LWR-UAM Phase III exercises were studied with a TRACE-PARCS coupling. More details about the Wilk's formula and its limitations can be found in [39]. The sampled cross-sections were used together with five additional CTF thermal-hydraulics parameters, three boundary conditions and nine manufacturing parameters that are defined in Tables 2 and 3 together with their probability density functions (PDF). In Table 2, the PDF are presented as multipliers on the nominal values, while in Table 3, the PDF are presented with the nominal values of the parameters along with their absolute standard deviations. The same standard deviation was used for all 235 U enrichments that vary from 4.0 to 5.0 w/o%. The manufacturing parameters were originally sampled with Sampler using a random sampling size of N = 100 and were applied in the Polaris lattice calculations, while in CTF, a consistent perturbation of relevant impacted parameters was performed for each sample, as explained in Section 2.4. The perturbations for all manufacturing parameters were applied to all fuel rods in the assembly. The CTF parameters and the boundary conditions were sampled with Dakota, using a random sampling of size N = 100.

Modeling
The core modeling is performed using a conventional two-step approach employed in LWR calculations. In the first step, lattice neutronics calculations are performed with Polaris. The ENDF/B-VII.1 cross-section library with SCALE's 56-group structure is used to generate the two-group macroscopic cross-sections with a 0.625 eV energy cutoff. SCALE's perturbation factors are used to generate the perturbed homoginized cross-section libraries. The perturbation factors are derived from the SCALE covariance library in 56 groups. The lattice calculations are performed for each unique fuel assembly design at 22 burnup points ranging from 0.0 to 65.0 GWD/MTU, 45 branch points and 8 history points. Fuel temperature, moderator density, boron concentration, and control insertion are the state parameters selected to define the branch and history points. The parameter coverage is given in Table 4 for the fuel assemblies. The four types of assembly layout are given in Figure 11.
Radial and axial reflector libraries were generated based on the LWR-UAM-Phase I specifications [40]. The radial reflector is constructed with three slabs neighboring a fuel assembly (Figure 12a). From left to right in Figure 12a is represented the fuel assembly, a 0.2 cm layer of coolant (first slab, blue), the shroud (second slab, yellow) and another layer of coolant (third slab, blue). The axial reflector is modeled with two slabs adjacent to a fuel assembly Figure 12b. From left to right in Figure 12b is represented the fuel assembly, the reflector (coolant, blue) and the core plate (second slab, yellow). The sub-regions in each slab represent the mesh grid used for the transport calculations in Polaris. The fuel assemblies in the radial and axial reflector models only provide a flux background using a critical spectrum to compute the homogenized cross-sections, so the reflector libraries only account for the reflector region. Reflective boundary conditions are implemented on the west side of the model, and vacuum boundary conditions are implemented on the east side of the model. The reflector libraries contain one history and 12 branches. The branch parameters for reflector libraries are the moderator density and boron concentration, and the same coverage as the fuel assemblies is used.  In the second step, the recently developed CTF-PARCS coupling was applied in an assembly-resolved model, leading to a total of 178 radial cells (177 fuel assemblies + 1 for the radial reflector). Axially, the core was divided into 28 meshes, 26 for the fuel active length and two for the top and bottom reflector. This level of nodalization is in line with LWR-UAM Phase III specifications and is visualized in Figure 13. In future uncertainty quantification studies, the uncertainties due to nodalization should also be considered using the Richardson extrapolation, as described in [41]. One representative fuel pin is modeled for each assembly with 10 radial nodes in the fuel and 2 radial nodes in the cladding. For each transient calculation, a steady-state coupled calculation is performed initially, and the core is brought to criticality by an eigenvalue fission source normalization. It is important to note that the normalization approach can have a non-negligible impact on the uncertainty quantification [42,43]. For the depletion calculation, a total of 24 depletion steps are performed. The first four steps are of 1 EFPD, and the next 20 steps are of 30 EFPD. The discussed modeling so far represents the standard fidelity of multi-physics calculations applied to uncertainty quantification studies, such as the one used in [22]. Compared to this standard fidelity and especially with regard to the thermal-hydraulics modeling in such multi-physics coupled calculations, some improvements in the thermal modeling have been implemented in this work, as discussed also in Section 2.3. The following improvements were implemented in CTF and used in the simulations presented in Section 4:

Results
In this section, we present the obtained results and statistical analysis for the different cases, using the previously discussed computational framework. The results are presented into three subsections. The first subsection is for the HFP BOC REA and includes two subdivisions: one for the steady-state results and one for the transient. The second subsection is for the HFP EOC REA and includes the same two subdivisions. The third and last subsection discusses the results for the cycle depletion from BOC to EOC. This creates a total of five distinctive case studies. For each case study, the results are analyzed in a hierarchy of functional complexity. At a first level, scalar quantities such as the maximum temperature during the transient are investigated with results about their Empirically estimated Cumulative Distributions (ECDF) and their first two statistical moments. For these quantities, the Spearman rank correlations are calculated, and as mentioned previously, only the correlations above the predefined absolute threshold of 0.35 are selected as statistically significant. It is important to stress that (a) these correlations do not necessarily imply causation, and (b) to derive any meaningful conclusion about sensitivities, a subjective judgment is necessary by understanding the underlying input uncertainties (e.g., correlations) and the physical phenomena modeled in each calculation. The second level of output quantities of interest are functional 1D quantities representing axial or temporal evolutions, such as the spatial maximum temperature trend during the transient, or the radially averaged axial power profile at HFP BOC. For these quantities, the mean is estimated together with the 95% percentile as the upper bound and the 5% percentile as the lower bound. The third level involves functional 2D radial output quantities such as the axially averaged power profile at HFP BOC. For these quantities, the mean and standard deviation are presented. All the selected scalar, functional 1D and functional 2D quantities for the steady state, transient and depletion cases are shown in Tables 5-7.   The steady-state calculation is performed for the HFP conditions as described in the specifications of LWR-UAM Phase III [24]. It is important to mention that an equilibrium Xenon and Samarium concentration is used. The results for the scalar outputs of interest are shown in Figure 14. The k eff shows an uncertainty of 0.44%, which is a result in the same range of that found in the literature for typical LWR uncertainty quantification studies (0.5%) such as in [40]. Small uncertainties of 0.81% and 0.37% are obtained for P max lin and D min cool , respectively. The T max fc shows a larger uncertainty of 3.14%. In order to understand the sources of these uncertainties, we estimate the Spearman rank correlation coefficients. The 56-group cross-sections that exhibit strong correlations for the k eff and P max lin are shown in Figure 15. For the k eff , we see that the 235 Uν p at low energies is the most significant input. This is in line with results from the LWR-UAM Phase I benchmark for steady state [40], since theν p has a direct contribution to the produced neutrons and thus to the k eff . The positive correlation is reasonable, since an increase of theν p will lead to an increase of the neutron population and then a higher k eff . As mentioned, correlation does not mean necessarily causation. This is the reason for example we cannot be sure that theν p at all energies below 1 keV is important since they are highly correlated. We can, however, deduce that since among all the independent isotopic crosssections only theν p seems to be highly correlated and since that this is in accordance with our understanding of the physical phenomena, then we are quite confident thatν p is the most important input at low energies for the k eff . With a similar reasoning, we analyze the cross-section correlations for P max lin , where we can see a very strong correlation for the 238 U elastic and inelastic scattering at high energies (above 1 MeV). From these two cross-sections, the inelastic scattering of 238 U is a quantity that has been shown to be one of the dominant inputs in uncertainty quantification studies such as in the LWR-UAM Phase I benchmark [40]. For the elastic cross-section, we have little information, but as we can see in Figure 16, it is strongly negatively correlated with the inelastic scattering of 238 U for this specific high-energy group. This is probably an indicator that the effect of both crosssections is essentially one effect, but since they are highly correlated to each other, both appear to be statistically significant. To better understand the impact of these scattering cross-sections, in Figure 17, we show the estimated Spearman correlation between the axially maximum linear power in every assembly and the 238 U elastic and inelastic crosssections. In this figure, only the maximum correlation in absolute value is shown for each assembly. What can be observed is that the negative correlation is also evident in these radial results and that two distinctive regions are present. One region is in the center with positive/negative contributions for the elastic/inelastic cross-sections and one region is in the periphery with negative/positive contributions for the elastic/inelastic cross-sections. The highly negative correlations between the elastic and inelastic cross-sections do not allow a separation of their impact on the two regions. What can be deduced is that based on the location of the maximum power, the impact of these cross-sections can be reversed. This is indeed what we observe in all the next studies, because the maximum local power will move toward an assembly in the periphery of the core. We attribute the presence of these two regions to the fact that the flux is normalized to 1 and thus to the differential impact of the scattering on the distribution of neutrons in the core, leading essentially to more neutrons pushed toward the periphery. For both k eff and P max lin , no other source of uncertainty from the manufacturing, boundary conditions or CTF inputs was found to be statistically significant.  6  8  10  12  14  16  18  20  22  24  26  28  30  32  34  36  38  40  42  44  46  48  50  52  54  56 Energy group For the other outputs, no cross-section showed a statistically significant correlation, but instead for the T max fc , the following inputs and their respective Spearman correlations were found: The most important input seems to be the cladding inner radius, which defines the gap thickness and thus gap conductance. A positive correlation is found because an increase of the gap width leads to a decrease of the gap conductance, less heat being extracted from the fuel resulting in a higher T max fc . For the same reason, the gap conductance is found to be the second significant input and with a negative correlation. Close to the cutoff threshold, we also find significant contributions from the fuel thermal conductivity and the fuel density, both with a negative sign. The former is simple to understand, because an increase in the thermal conductivity increases the heat extracted from the fuel and thus decreases the T max fc . For the latter, the effect is more difficult to explain, because it is an input to both Polaris and CTF. One explanation can be the contribution of the fuel density through the thermal conductivity, since an increase of the density will lead to an increase of the thermal conductivity. Another explanation could be the self-shielding on the fuel pellet outer layers and the effect on the neutron utilization: the change in fuel density in Polaris increases both the number density of 235 U and 238 U. The two isotopes have competitive effects on the neutron utilization, i.e., 235 U will tend to increase the fission rate and 238 U will increase the absorption rate due to resonances. The increase of 238 U density may be dominant as compared to the 235 U density. The two effects (conducitivity and neutronics) may be both contributing, or one of the two effects predominates. To verify those postulates, two tests can be performed for future work: (a) the density of the fuel can be modified in Polaris in a few samples, and the effect on T max fc can be evaluated without modifying the conductivity in the thermal-hydraulics calculations, and vice versa, (b) the fuel density in Polaris can be kept constant, but the fuel conductivity can be sampled in the thermal hydraulics calculations to evaluate the effect of the conductivity alone.
Regarding the D min cool , only the mass flux and the inlet temperature were found to be statistically significant. The impact of both is obvious for the D min cool and their Spearman correlations are the following:
The functional 1D results for the different axial profiles are shown in Figure 18. The mean is estimated together with the upper and lower bounds defined as the 95% and 5% percentiles, respectively. We can see that the linear power P ave lin,ax has a shape close to a cosine slightly shifted toward the bottom with a very small uncertainty. A similar trend is observed in [22]. The T ave fc,ax follows the same trend as P ave lin,ax with the upper and lower bounds ranging at ∼50 K. For the T ave cool,ax and D ave cool,ax , the expected trends are found with relatively small uncertainty bounds. The functional 2D results are presented in Figures 19 and 20 where the radial 2D mean and standard deviation of P ave lin,rad and T ave fc,rad are shown, respectively. We can see that the highest linear power is in the location E10 and all its symmetrical locations, while the largest uncertainty is ∼2.7% located in the central assembly. The centerline temperature profiles follow the linear power profile with the hot spot being at the location of E10. The highest uncertainty is ∼ 3.2% and is located in the fuel assemblies neighboring the central assembly.

Transient
The REA calculation is performed for the HFP conditions described in the specifications of LWR-UAM Phase III [24] starting from critical initial conditions reached by using fission source normalization. As for the steady state, an equilibrium Xenon and Samarium concentration is used. The results for the scalar outputs of interest are shown in Figure 21.
The ρ in has a mean value of 0.391$ and an uncertainty of 12.5%. The P max lin mean has increased by a factor of 3 compared to the steady state with a significant larger uncertainty of 12.8%. The T max fc follows the same behavior with an increase of ∼700 K over the steady-state average but with a smaller uncertainty of 1.7%. This decrease is attributed to the fact that at steady state, the gap conductance could be open or closed depending on the initial gap width, while in the transient, the expansion of the fuel pin always leads to a gap closure and consequently decreases the impact of the gap width and the gap conductance. The H max f shows a mean value of 91 cal/g with an uncertainty of 2.0% similar to the T max fc . Investigating now the sources of uncertainties for these scalar output quantities, the Spearman rank correlation coefficients are calculated. Two clusters of trends can be observed coming from the cross-sections. The first cluster corresponds to ρ in and P max lin where the 56-group cross-sections show the same strong correlations for both of them as can be seen in Figure 22a,b. Similar to the steady-state results, the 238 U elastic and inelastic scattering cross-sections at high energies are statistically significant. Interestingly enough, they show an inverse sign for the correlation with regard to the steady-state results. This is attributed to the different location of the maximum linear power between the steady-state and transient results. The location moves from the center toward the periphery, and as can be seen in Figure 17, this leads to an inverse correlation. Although these cross-sections show strong correlations, the 238 Uν d andν p at high energies show even stronger correlations. The opposing signs between theν d andν p are attributed to their strong negative correlation in these energies, as can be seen in Figure 23. Further analysis of the 238 Uν d andν p uncertainties leads to the conclusion that the effective delayed neutron fraction is strongly correlated with these cross-sections and specifically with opposing correlations signs compared to ρ in . This means that probably the 238 Uν d andν p crosssections determine the uncertainty of the effective delayed neutron fraction, which is in the denominator of the ρ in . What we can deduce at the end from this analysis is that the uncertainty of the ρ in and P max lin is determined by the ejected control rod worth and by the effective delayed neutron fraction. The former uncertainty is mainly determined by 238 U inelastic and elastic scattering uncertainties, while the latter is mainly determined by the 238 Uν d andν p uncertainties. The remaining two quantities T max fc and H max f belong to the second cluster of trends and the cross-sections with stronger correlations are only the 238 U elastic and inelastic scattering cross-sections, as can be seen in Figure 22c,d. The functional 1D results over time for P max lin,t , T max fc,t and H max f,t are presented in Figure 24. It can be seen that P max lin,t shows large uncertainties around the maximum, as seen also in the scalar results, but then the uncertainty decreases rapidly. The location of the maximum linear power is located during the whole transient in the assemblies D13 and E12 that are neighboring the control rod ejection assembly (D12). The T max fc,t behavior shows that the maximum fuel centerline temperature stabilizes after 35 s. Initially, the hot spot location coincides with the location of the maximum linear power, but at around 10 s, as indicated by the non-linearity in the behavior, the hot spot location moves to the assembly where the rod was ejected, D12. This happens because the D12 assembly has higher burnup, as can be seen in Figure 10, and thus, the deteriorating thermal conductivity leads to larger temperatures later on in the transient. The uncertainty of T max fc,t creates bounds of ∼80 K. The H max f,t behavior follows closely the behavior of T max fc,t .

Steady State
The steady-state calculation is performed for the HFP conditions described in the specifications of LWR-UAM Phase III [24] at the EOC. The results for the scalar outputs of interest are shown in Figure 27. The k eff shows a slight increase compared to the BOC case, while the P max lin shows a significant increase. The uncertainty of D min cool remains similar, while the T max fc decreases almost by half. The underlying sources of these uncertainties are investigated through the estimation of the Spearman rank correlation coefficients. Some of the 56-group cross-sections show strong correlations for the k eff , P max lin and T max fc as shown in Figure 28, while no cross-section is identified as significant for D min cool . For k eff , a very interesting observation can be made since we do not see anymore the 235 Uν p being the dominant input, but instead, we see contributions from the 238 U elastic and inelastic scattering at high energies and 239 Pu fission and capture at low energies. This is attributed to the increased average burnup of the core and thus to stronger 239 Pu buildup that contributes a significant part of the neutrons production. This explains also the positive sign for the fission cross-section and the negative sign for the capture cross-section. The results for P max lin are similar to the BOC case with the 238 U elastic and inelastic scattering at high energies being the statistically significant inputs. The only difference is the inverse sign of the correlations, which is something attributed to the same reason as for the BOC transient results. The location of the maximum moves from the center region (E10) at BOC to the periphery (C10) at EOC. These regions, as shown in Figure 17, have opposite correlations for the 238 U elastic and inelastic scattering. For T max fc , we see the same cross-sections as important, but in the BOC calculations, these were not identified. This is because probably at BOC, their correlation was below the cutoff threshold, while at EOC, it increased slightly above the threshold.
From the other inputs, none was found statistically significant for the k eff and P max lin . For T max fc , the following were found to be significant: • Fuel density: −0.45. • Fuel thermal conductivity: −0.46 The fuel density and thermal conductivity correlation increased for T max fc compared to BOC, while we do not see as important anymore the gap heat transfer-related inputs. This is attributed to the increased average burnup at the EOC. Regarding the D min cool , the same inputs with the BOC case are found significant with their correlation being: • Mass flow rate: 0.41. • Inlet temperature: −0.79.
The functional 1D results for the different axial profiles are shown in Figure 29. The main impact is on the behavior of P ave lin,ax and T ave fc,ax , where a strong bottom peaking is observed. A similar trend is observed in [22]. The uncertainties on all the functional quantities are small, as can be seen from the tight upper and lower bounds in Figure 29.  The functional 2D results are presented in Figures 30 and 31 where the radial 2D mean and standard deviation of P ave lin,rad and T ave fc,rad are shown, respectively. We can see that the peaking of the linear power is in the location E10 similar to the BOC case but also in some peripheral assemblies such as C10 assembly. The P ave lin,rad uncertainty is up to ∼2.5% and in the central assembly. The T ave fc,rad follows the mean profile of P ave lin,rad and exhibits uncertainties up to ∼1.4% located in the fuel assemblies in the periphery of the core.

Transient
The REA calculation is performed for the HFP conditions described in the specifications of LWR-UAM Phase III [24] at EOC and starting from critical initial conditions using fission source normalization. The results for the scalar outputs of interest are shown in Figure 32. The ρ in mean value has decreased to 0.357% but its uncertainty increased to 16.4%. The P max lin mean has decreased compared to the BOC REA due to the smaller inserted reactivity, and its uncertainty remained similar at 11.9%. The T max fc shows a significant decrease of ∼400 K compared to the BOC REA with a similar uncertainty of 1.7%. This decrease is attributed to a less extreme REA transient at EOC conditions. The H max f decreases significantly its mean to 76 cal/g following the T max fc trend with a similar uncertainty of 1.6%.
In order to understand the change in the obtained uncertainties, the Spearman rank correlation coefficients are calculated. Similar to the BOC REA, two clusters of trends can be observed concerning the cross-sections. The first cluster corresponds to ρ in and P max lin , where the 56-group cross-sections that exhibit strong correlations are shown in Figure 22a,b. The difference with respect to BOC REA is that the 238 U elastic and inelastic cross-sections at high energies contributions have decreased significantly, while the 239 Puν d at low energies shows a strong statistical significance. Apart from that, the dominant cross-sections are still the 238 Uν d andν p at high energies. For T max fc and H max f , the same cross-sections as for the BOC REA case appear, namely, the 238 U elastic and inelastic cross-section at high energies, as can be seen in Figure 33c,d. From all the other sources of uncertainties, only the fuel density is statistically significant for the T max fc with a Spearman rank correlation of −0.35. As explained in the HFP BOC case, this correlation is attributed to the impact of the fuel density on the thermal conductivity.
The functional 1D results over time for P max lin,t , T max fc,t and H max f,t are presented in Figure 34, and similar conclusions to the BOC REA can be reached. The P max lin,t shows large uncertainties around the maximum with the uncertainty decreasing rapidly afterwards. The location of the maximum linear power is located during the whole transient in the assemblies D13 and E12 that are neighboring the control rod ejection assembly (D12). The T max fc,t behavior shows that the maximum fuel centerline temperature stabilizes after 30 s and is now located at the assembly E12. The uncertainty of T max fc,t creates bounds of ∼40 K. The H max f,t behavior follows closely the behavior of T max fc,t . The functional 2D results at the end of the transient for P ave lin,end and T ave fc,end are presented in Figures 35 and 36. The P ave lin,end peaking moves toward the center of the core compared to the BOC REA. This applies to the T ave fc,end as well that now follows more closely the P ave lin,end behavior. Concerning the uncertainties, the P ave lin,end shows a similar standard deviation up to 2.5% in the central assembly, while the T ave fc,end decreased by half a standard deviation of up to 1.5% in both central and peripheral assemblies.

Depletion
The depletion from BOC to EOC calculations is performed for the HFP conditions described in the specifications of LWR-UAM Phase III [24]. The results for the scalar outputs of interest are shown in Figure 37. The obtained Bu ave mean is 40.381 GWD/MTU, which is a value relatively close to the reference value of 42.06 GWD/MTU found in the LWR-UAM specifications. The Bu max mean is 62.78 GWD/MTU. The uncertainties for Bu ave and Bu max are 0.3% and 0.5%, respectively. The Spearman rank correlation coefficients between the 56-group cross-sections and Bu ave do not show any statistical significant cross-section. For the Bu max , however, four cross-sections are found significant and are shown in Figure 38. The main cross-sections are the 238 U elastic and inelastic scattering at high energies. The sign of these correlations is similar to the BOC steady-state results and reversed compared to the BOC transient and EOC steady state and transient. This, as explained in Figure 17, is due to the location of the maximum that occurs in the central region of the core. At a second order, the 235 U elastic and fission at low energies are identified, but they are very close to the cutoff threshold. From the remaining inputs, only the fuel density is significant, with a dominant role for Bu ave showing a correlation of −0.91, while for the Bu max , it has a correlation of −0.41. The impact of the fuel density and its negative sign is obvious, because the burnup step calculated by PARCS has the fuel density in the denominator. The functional 1D results are split into axial and temporal. The results for the temporal quantities are shown in Figure 39. As we can see, the Bu max t is increasing linearly for most of the depletion with small uncertainty bounds. The P max lin,t decreases during the depletion with small uncertainties indicating the decrease of the axial offset. The T max fc,t evolution over depletion indicates uncertain bounds of ∼80 K. The T max fc,t trend shows a decrease from 1700 K to 1540 K during half of the depletion and then a stabilization around this value. The results for the axial quantities are shown in Figure 40. We can see that the linear power P ave lin,BOC has a shape close to a cosinus with a slight shift toward the bottom and exhibits a very small uncertainty. At the EOC, the P ave lin,EOC shows a double peaking profile with one peak on the bottom and one on the top with very small uncertainties. The obtained radially averaged axial burnup profile Bu ave ax,EOC is compared in Figure 40c with the reference results of the LWR-UAM specifications. The comparison indicates a relatively good agreement for most of the fuel active height and larger discrepancies close to the top and bottom reflectors.
The functional 2D radial results for the obtained axially averaged burnup distribution Bu ave rad,EOC are presented in Figure 41. The estimated standard deviation reaches values up to 1.2% in locations close to the periphery. In Figure 42, the discrepancy between the nominal Bu ave rad,EOC and the reference from the LWR-UAM specifications is shown. The errors span from −3.5% to 1.5%, indicating a reasonable agreement.  Figure 42. Average radial burnup discrepancy between the predicted P ave lin,BOC and the benchmark reference.

Conclusions
In this work, the developed CTF-PARCS core multi-physics computational framework for efficient Light Water Reactor (LWR) steady-state, depletion and transient uncertainty quantification is presented. Novel features of this framework include a versatile Python pre-processing script for streamlining the macroscopic cross-section generation process and the capability to estimate correlations between outputs of interest and the 56-group microscopic cross-sections, fission yields and all other multi-physics input parameters.
The developed computational framework is applied to three exercises of LWR-UAM Phase III that cover steady-state, depletion and transient calculations. The results show that the maximum fuel centerline temperature upper bound, defined as the 95% percentile, during the Beginning Of Cycle (BOC) Rod Ejection Accident (REA) is ∼2600 K, while during the End Of Cycle (EOC) REA, it is ∼2200 K. The results, in terms of correlations with the 56-group microscopic cross-sections, show that few isotope-reactions are found as statistically significant with strong correlations: • At BOC steady state, the 235 Uν p at low energies (<1 keV) and the 238 U inelastic and elastic scattering at high energies (>1 MeV, which corresponds to the 238 U inelastic scattering threshold). • At BOC REA, the 238 U inelastic and elastic scattering at high energies as well as the 238 Uν p andν d at high energies. • At EOC steady state, the 238 U inelastic and elastic scattering at high energies as well as 239 Pu fission and capture cross-section at low energies (<10 eV). • At EOC REA, the 238 U inelastic and elastic scattering at high energies, the 238 Uν p and ν d , and 239 Puν d at low energies (<100 eV).
The results highlight the efficiency and robustness of the computational framework and the ability to extract information, although subjective, up to the microscopic crosssections, which is something that can be very useful for identifying future experiments to reduce the nuclear data uncertainties.
Future steps based on this work will involve validating the computational framework and estimating the impact of different modeling approaches. Furthermore, parallel domain decomposition could be implemented in CTF to accelerate the calculations and allow even a pin-cell level of multi-physics coupling. The uncertainty quantification approach could also be investigated more thoroughly in order to derive a more objective way of discussing the obtained correlations for the microscopic cross-sections.

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