High to Low pellet cladding gap heat transfer modeling methodology in an uncertainty quantification framework for a PWR Rod Ejection Accident with best estimate coupling

High to Low modeling approaches can alleviate the computationally expensive fuel modeling in nuclear reactor’s transient uncertainty quantification. This is especially the case for Rod Ejection Accident (REA) in Pressurized Water Reactors (PWR) were strong multi-physics interactions occur. In this work, we develop and propose a pellet cladding gap heat transfer (Hgap) High to Low modeling methodology for a PWR REA in an uncertainty quantification framework. The methodology involves the calibration of a simplified Hgap model based on high fidelity simulations with the fuel-thermomechanics code ALCYONE1. The calibrated model is then introduced into the CEA developed CORPUS Best Estimate (BE) multi-physics coupling between APOLLO3 R © and FLICA4. This creates an Improved Best Estimate (IBE) coupling that is then used for an uncertainty quantification study. The results indicate that with IBE the distance to boiling crisis uncertainty is decreased from 57% to 42%. This is reflected to the decrease of the sensitivity of Hgap. In the BE coupling Hgap was responsible for 50% of the output variance while in IBE it is close to 0. These results show the potential gain of High to Low approaches for Hgap modeling in REA uncertainty analyses.


Introduction
The improvement of nuclear reactors computational modeling together with the increasing requirements in safety analyses lead to the development of Best Estimate Plus Uncertainty (BEPU) approaches [1]. BEPU is a systematic approach where Best Estimate (BE) codes are used to calculate safety criteria with an estimation of their respective uncertainties. BEPU is a step forward from conservative approaches allowing for better quantification of the safety margins. BE codes are considered the codes that take into account the most important underlying physical phenomena, governing the different studied scenarios, with reasonable approximations. Different BE codes usually model different physics. The correct modeling of transient scenarios can require the coupling between different BE codes due to the multi-physics interactions. This creates many challenges for the uncertainty quantification (UQ). With the term UQ we mean the estimation of the statistical distributions of the outputs of interest (in particular, their means and standard deviations) and the global sensitivity analysis of the outputs (in particular, * e-mail: delipei.gregoryk@gmail.com their Shapley indices). Many code evaluations are required for the estimation of the statistical quantities and of the sensitivity indices rendering thus the computational expensive modeling the main limitation. A large scope of undergoing research, in both academia and industry, is the multifidelity High to Low modeling approaches that can reduce in an efficient way the computational cost while retaining a good predictive modeling.
The multifidelity approaches combine different sources of information coming from high fidelity models, low fidelity models and experiments. The fundamentals work of [2] and [3] created the theoretical framework that lead to papers such as [4] in structural dynamics, where a cheap surrogate model based on Polynomial Chaos was considered as a low fidelity model and corrected by a high fidelity code in order to accelerate the Bayesian calibration process using experimental data. In [5] a methodology for combining information from multifidelity codes and experiments in a Bayesian hierarchical model was developed using Gaussian processes. In [6] a novel machine learning based multifidelity approach was discussed in the context of Bayesian optimization. A comprehensive and detailed survey of multifidelity methods in uncertainty quantification, statistical inference and optimization can be found in [7]. In nuclear and more particularly in thermal-hydraulics the result of the PREMIUM project [8] indicated a strong user effect in the input uncertainty quantification. For this purpose in [9] a set of good practices is discussed regarding the input uncertainty quantification.
In this work we focus on the study of the Rod Ejection Accident (REA) in Pressurized Water Reactors (PWR) with a multi-physics coupling in an UQ framework. During the REA strong multi-physics coupling effects occur between neutronics, fuel-thermomechanics and thermalhydraulics. This necessitates a multi-physics modeling to take into account the interdisciplinary interactions increasing significantly the computational cost. The uncertaintly quantification for REA becomes thus very challenging. In this coupling the fuel-thermomechanics modeling is the most expensive since it solves at the fuel rod scale the thermal conduction equations coupled with the thermo-mechanical ones at a fine discretization. A High to Low modeling is necessary to introduce and use a simplified model for fuel-thermomechanics reducing the computational cost. The high fidelity modeling manages to model complex thermo-mechanical and physico-chemical phenomena occurring during a REA. Some examples are the high burn-up structure, fission gas releases in the pellet-cladding gap and the pellet cladding mechanical interaction (PCMI). All these phenomena impact the fuel rod behavior through the pellet cladding gap heat transfer (H gap ). For this reason the low fidelity fuel-thermomechanics modeling usually does not model explicitly these phenomena but directly use a predefined H gap value. There is, therefore, a strong interest in developing High to Low models for H gap in REA. In this work, we develop and propose a methodology to develop such High to Low models for H gap . The high fidelity fuel-thermomechanics code ALCYONE1 [10] is used to calibrate a simplified H gap model based on fuel thermal expansion. Compared to previous works such as [5], the low fidelity model is not a fitted surrogate but a physics based model that is calibrated in a preliminary way by minimizing the discrepancy with ALCYONE1 predictions. The model is then introduced into the CEA developed CORPUS [11] Best Estimate (BE) multi-physics tool, where core neutronics code APOLLO3 R [12] is coupled with core thermal-hydraulics code FLICA4 [13]. This creates an Improved Best Estimate (IBE) modeling of REA that is used for UQ. To the knowledge of the authors this is one of the first attempts to address H gap High to Low modeling in transients. Some considerations have been discussed in the context of the UAM (Uncertainty Analysis in Modelling) benchmark [14]. The benchmark consists of three phases starting from stand-alone neutronics uncertainty analysis for steady state and reaching up to full coupling at the system level for transient scenarios. In this benchmark the H gap is treated through lookup tables for the steady state and evolution calculations. However, there is not a dedicated treatment for the transient scenarios and this work and its prospects could be potentially integrated in the benchmark. Additionally in [15] the pellet cladding modeling was optimized including experimental data. The model's uncertainty was estimated as well.
In Section 2 we detail the two available coupling frameworks at CEA for the modeling of REA. The computationally expensive Best Effort coupling and the cheaper BE coupling. In Section 3 the specifications of the case study are presented. The geometry of the studied PWR core and the selected nodalization are provided together with the REA characteristics and the inputs and outputs of interest. The High to Low H gap modeling methodology is based on the results of an UQ analysis performed with the BE coupling (BE UQ). In Section 4 we present these results. In Section 5 we present the calibration of the High to Low Hgap model. The simplified H gap model is detailed and the calibration results are presented. As mentioned above, the obtained H gap model is introduced into the previous coupling creating the IBE. In Section 6 an UQ analysis is performed with the IBE coupling (IBE UQ) and the results are compared with the BE results of Section 4. Finally, we end this article with the main conclusions and perspectives of this work.

REA coupling framework
The REA is initiated by a control rod ejection due to mechanical malfunction inserting positive reactivity in the core. As a consequence power increases adiabatically until the fuel temperature starts increasing as well introducing a Doppler negative feedback that creates a power peak. The power then continues to decrease and when the heat flux reaches the moderator its density will start decreasing adding another negative feedback due to the negative moderator temperature effect. Depending on the core state at the moment of the ejection the power evolution and its damage on the first containment barrier (cladding) can vary significantly. During the transient strong multi-physics interaction effects occur between neutronics, fuel-thermomechanics and thermal-hydraulics necessitating a multi-physics coupling framework.
At CEA two different coupling schemes are established [16] for the modeling of REA based on CORPUS multiphysics tool as can be seen in Figure 1a. APOLLO3 R code is used for core neutronics, FLICA4 code for core thermal-hydraulics and ALCYONE1 code for pin fuelthermomechanics. The Best Effort scheme involves the full coupling between the codes. At each time step the codes exchange 3D fields. APOLLO3 R uses the fuel Doppler temperature T f from ALCYONE1 and moderator temperature T m from FLICA4 to compute the power generated in the fuel P f and the fluid P l . FLICA4 uses P l from APOLLO3 R and the cladding wall heat flux Φ w from ALCYONE1 to compute the T m and the external cladding wall temperature T ext c . ALCYONE1 uses the P f from APOLLO3 R and the T ext c from FLICA4 to compute the T f and Φ w . The Best Effort scheme models each discipline accurately but it is computationally very expensive for a full core REA and thus its use for the uncertainty propagation and global sensitivity analysis in the UQ framework is prohibited. For this reason the BE scheme is used where the coupling is reduced between APOLLO3 R and FLICA4. The complex fuel-thermomechanics phenomena

PWR geometry and modeling
The REA is studied in a large scale PWR core. The geometry has an 1/8 symmetry illustrated in Figure 2. It consists of 193 fuel assemblies with U O 2 and U O 2 − GdO 3 fuel compositions. Two different types of control rods are inserted at different depths. The black rods (B) with high neutrons absorption that are typically used for the shutdown of the reactor and the grey rods with less neutrons absorption that are used in the day to day reactivity control. The core is at Hot Zero Power (HZP) conditions at the end of the cycle. Around the fuel assemblies there is one ring of water reflector assemblies. The total height of the core is 468.72 cm with a bottom and top reflector of 21 cm leading to a fuel active height of 426.72 cm. Each assembly is a 17 × 17 lattice of fuel pins with pitch 21.504 cm. The control rod that will be ejected initiating the REA is located on the periphery as highlighted in Figure 2. It is inserted 97 cm from the top into the fuel active region. Due to the extraction of the control rod there is 1/2 symmetry for the REA.
The REA in the PWR core geometry is modeled using the BE and IBE coupling scheme discussed in Section 2. For core neutronics, APOLLO3 R code is used with a two group Diffusion approximation and void boundary conditions on the neutron current. The two group macroscopic cross-sections are parameterized in burn-up, boron concentration, moderator density and fuel temperature. The radial discretization is at the level of the quarter of assembly. For the axial discretization 34 meshes are used of which 30 are for the fuel active height and the rest for the top and bottom reflectors. For core thermalhydraulics, FLICA4 is used with a 4 equations porous modeling and a multi-1D axial flow approximation. The system of 4 equations consists of: mixture mass balance, vapor mass balance, mixture momentum balance and mixture energy balance. The boundary conditions are the inlet mass flow and enthalpy and the outlet pressure. For the radial discretization one thermal-hydraulic channel is used for each assembly. For the axial discretization  only the fuel active part is modeled using 30 meshes in accordance to APOLLO3 R modeling. One average fuel pin per thermal-hydraulic channel with 1D radial modeling is used in FLICA4 with a discretization of 25 radial regions for the fuel and 3 for the cladding. For the time discretization of the REA an adaptive end time is adopted based on the integral power evolution. For each transient in the UQ when the power surpasses half its nominal value then a SCRAM signal is sent. It is considered that from this time on 0.6 s are needed in order for the SCRAM to take place and end the modeling of the transient. The incremental time step is constant and equal to 0.001 s. The control rod is ejected in 0.1 s.

Initial state and reference transient
The PWR core is critical at the end of the cycle with HZP condition, meaning that the temperature is around 290 deg C and the power negligible (3.8W ). Since the core is at the end of the cycle the boron concentration is quite low at 95.5 ppm. Other characteristic conditions of the core are provided in Table 1. There is a burn-up distribution resulting from core evolution calculation performed in a previous internal CEA study. The obtained radial burn-up distribution at the end of cycle is illustrated in Figure 3. The burn-up radially averaged at the level of the assembly ranges from 10 to 52 GWd/t. The location of the ejected control rod is close to the periphery. In the core there is a radial and axial Xenon distribution that increases the control rod worth leading to an increased reactivity insertion and thus a violent prompt driven transient. The reference (without uncertainties) REA characteristics obtained with the BE coupling are presented in  Table 2 and Figure 4. The control rod worth is ρ worth = 1.2 $ indicating a prompt driven transient. In the Figure we can observe the created power pulse of width Γ = 38 ms with a maximum power of P max core = 2.54P nom at instant t max = 292 ms. The nominal power is P nom = 3800 MWth. Additionally, the F xyz deformation factor evolution in time is plotted. It starts from a value of 5 and reaches up to 25 when the control rod is fully ejected (0.1 s).

Input and outputs identification
Before presenting the BE UQ results we define the uncertain inputs that will be considered together with the outputs of interest. The inputs and outputs for the three disciplines are defined in Table 3. The statistical distributions of the inputs are presented in [17] and are included in Appendix A. In neutronics the two-group macroscopic cross-sections and the kinetic parameters are considered with a multivariate normal distribution. The rest of the inputs are considered independent of the neutronic inputs and between them. The thermal-hydraulic input distributions are mainly based on CEA experts opinions and are applied as random multiplicative factors with mean 1 on the different models. The UAM recommendations are used for the fuel-thermal inputs, i.e. for the thermal conductivities and specific heat capacities. The H gap is the quantity of focus in this work. In the BE coupling scheme H gap is an uncertain constant with a uniform distribution between the value for a complete open gap (2e 3 W m −2 K −1 ) and the one for a pellet-cladding contact (5e 4 W m −2 K −1 ). In Table 3. Inputs and outputs identification by discipline: neutronics, fuel-thermal and thermal-hydraulics.

Inputs (22)
Disappearance cross-section of group g NF g (2) νxfission cross-section of group g D g (2) Diffusion coefficient of group g S 1→2 Scattering cross-section of group 1 to 2 IV g (2) Inverse In neutronics we consider the maximum in time and space of the local linear power during the REA P max lin and the radial linear power distribution P 2D lin at the time instance and axial plane of P max lin . In fuel-thermomechanics we selected the maximum in time and space of the local stored enthalpy H max f while in thermal-hydraulics the minimum in time of the distance to Departure from Nucleate Boiling (DNB) DN B min . The latter output is defined as the difference between the Departure from Nucleate Boiling Ratio (DNBR) and the DNB threshold R crit . The UQ methodology is detailed and tested in [17] and [18]. The methodology consists of an initial screening of the inputs and the use of kriging models [19] to perform uncertainty propagation and global sensitivity analysis using Shapley indices [20,21]. More details about the methodology are discussed in Appendix B. The obtained results for the uncertainty propagation and global sensitivity analysis will be presented directly for both BE and IBE coupling schemes.

Best estimate uncertainty quantification
In this section we present the results of the BE UQ for the REA. The uncertainty propagation is performed by training kriging models between each scalar output and a reduced number of screened inputs. The screening is performed with a method based on HSIC statistical significance tests [22,23] presented in [17]. The kriging models are trained on a Latin Hypercube Sampling (LHS) of size 250. For the functional output Principal Component Analysis (PCA) [24] is used and kriging models are trained for two principal components that represent 95% of the outputs variances. For all the quantities of interest the obtained Q 2 of the kriging models, evaluated on a separate LHS of size 125, is close to 1 with smallest value 0.98 for DN B min and 0.945 for the second principal component of P 2D lin . More detailed results are provided in Appendix B.
The kriging models are used for brute force Monte Carlo uncertainty propagation with 10 5 samples and the results are presented in Figure 5. In each histogram we plot also the pdf of a normal distribution with the estimated mean and standard deviation (blue line) for a visual comparison of the normality. In Appendix B we present also the estimated relative standard deviations for the scalar quantities directly using the training LHS. It is important to notice that it is not recommended to use surrogate models for accurate estimation of the tails of the output pdf since they consist of rare events that the surrogate model has not been trained for. We observe a quite large relative standard deviation for P max lin of around 57%. The result for P 2D lin shows that the spatial distribution of the relative standard deviation does not vary significantly radially. For H max f a 20% relative standard deviation is obtained and the distribution is approximately normal. For DN B min a large mean value is obtained with large relative standard deviation of 57% resulting in a very small probability of reaching boiling crisis.
For the global sensitivity analysis, Shapley indices are estimated for the scalar outputs and aggregate Shapley indices for the functional ones as discussed in [17]. All the Shapley indices are estimated using the kriging models. For a more detailed focus on the Shapley indices estimation with surrogate models we refer to [25]. The R language [26] package "sensitivity" [27] (function "shapleyPermRand") was used for the estimation of the Shapley indices and their confidence intervals, that quantify the uncertainty of the estimation method. The results are presented in Figure 6. The inputs in this Figure are the selected ones by the screening process of the UQ methodology. For the outputs P max lin , P 2D lin and H max f the β ef f is responsible for 50% of the output variance with the other 50% attributed to the cross-sections T D 1 , D 1 , S 1→2 and IV 1 . All these cross-sections are highly correlated to each other and thus it is difficult to distinguish their separate contributions. For the output DN B min the H gap is the dominant input responsible for 50% of the output variance while the remaining 50% is attributed to the cross-sections and the β ef f . This is due to the large H gap uncertainty ranges and thus gives the incentive to improve its modeling.

Introduction
The gap heat transfer (H gap ) modeling in the BE coupling is performed through a constant value during the REA with a uniform distribution over a large interval. This is one of the most important modeling differences between Best Effort and BE coupling. There is a strong interest thus to improve the H gap modeling and to introduce it into a refined coupling that we will call Improved Best Estimate (IBE). In this subsection we address this challenge by developing a High to Low H gap modeling methodology. The methodology involves the calibration of a simplified H gap model that is based on fuel thermal expansion. We consider that this model is adequate for the REA and especially for the gap closing phase. The calibration is performed through decoupled ALCYONE V1.4 calculations with imposed power evolution. In this Section we first detail the model with its calibration parameters.
Afterwards, we present the methodology where we discuss issues such as how the power pulses are selected, how many H gap models will be created for the different fuel assemblies etc. Finally, we present the calibration results for the PWR core.

Pellet cladding gap heat transfer simplified model
The sharp power increase in the REA leads to a corresponding sharp fuel temperature increase. We assume that the H gap evolution is driven by the gap closing due to fuel thermal expansion and by the gas conductivity evolution in the gap. This is used to derive the simplified formulation for H gap defined by equations (1a)-(1c).
where: λ g is the gas conductivity in the gap and λ g,init = H init gap einit its initial value prior to the REA. The latter is calculated by the initial gap heat transfer H init gap and initial gap width e init . -T g is the gas temperature and T g,init the initial gas temperature prior to the REA. -E f is the energy stored in the fuel during the REA.
θ 1 and θ 2 are two calibration parameters that have to be estimated. e is the pellet-cladding gap width. It is assumed that only the fuel expansion is responsible for the gap evolution. r c,init is the initial internal cladding radius prior to the REA. r f = r f,init α f (T f ) is the fuel external radius. The fuel expansion is modeled using the fuel expansion coefficient α f (T f ) in ALCYONE1 which is a cubic function of the fuel temperature T f .
The H gap predicted by the proposed simplified model is based on fuel thermal expansion and depends on the evolution of the λ g and e. The λ g is considered a linear function of T g and E f . The latter allows to include a historical effect on the H gap . The two calibration parameters to be determined θ 1 and θ 2 are the coefficients of T g and E f respectively. The e evolution is assumed to depend only on the fuel thermal expansion while the cladding radius remains constant. For the modeling of the gas temperature T g the average between the external fuel temperature and the internal cladding temperature is used.

Methodology
Having defined the simplified H gap model the next step is to calibrate it. To this purpose we developed the High to Low methodology illustrated in Figure 7.
The starting point is the application of the UQ methodology in the BE coupling of Section 4. The temperature evolution in the fuel depends on the fuel assembly burn-up and the power seen by this assembly at its position in the core. Since we need to build a H gap model for every fuel spatial mesh (1080 meshes) a clustering is necessary. In step 1.1, we cluster the assemblies by similar radial average burn-up. This means that we will construct one H gap model for every identified cluster. While each cluster will have the same model, their different fuel initial conditions and the different power seen during the REA will lead to a 3D H gap . The next step 1.2 consists in selecting representative pulses from the BE UQ presented in Section 4. The pulses must cover most of the possible H gap variations inside the cluster due to both statistical and spatial aspects. For the statistical aspect three pulses are extracted: namely the ones which have produced the mean, the 97.5% upper quantile, and the 2.5% lower quantile of P max lin . For the spatial aspects, one selects the assembly that achieves the maximal power for the pulses corresponding to the mean and upper 97.5% quantile and one selects the assembly that achieves the minimal power for the pulse corresponding to the lower 2.5% quantile. For the spatial aspects, when the mean and the upper 97.5% quantile are imposed then the assembly seeing the maximum power is selected. Correspondingly, when the lower quantile is imposed the assembly seeing the lowest power is selected. This creates for each H gap model different representative axial and temporal profiles of linear power and external cladding temperature.
In step 2, the selected profiles are extracted and imposed in a decoupled ALCYONE1 REA transient calculation. One representative fuel pin is modeled. The resulting temperature and gap heat transfer profiles from ALCYONE1 are extracted and used for the H gap model calibration. In step 3.1 the calibration is carried out by minimizing an objective function using the BFGS optimization method. The mean square error on the H gap maximal and final values during the REA for each axial slice is considered as the objective function. Once the parameters that minimize this function are estimated the final step 3.2 is to quantify the calibrated model's uncertainty. The two main sources of uncertainties are the initial conditions and the calibrated parameters. The former one is quantified as two multiplicative factors on the initial gap width and the initial H gap with normal distribution with mean one and standard deviation 0.1. This is a result of a previous uncertainty propagation for fuel evolution calculations with ALCYONE1. The results also showed that the initial gap width and H gap are almost fully negatively correlated. This leads to make the assumption that the two coefficients are fully negatively correlated (ρ = −1) rendering thus one effective uncertain quantity for the initial conditions H g,i . The latter uncertainty source is the calibrated parameters. For simplifying their uncertainty quantification, they are considered as fully positively correlated (ρ = 1) with uniform distributions. The bounds of the distributions are calculated in order to account for the calibration error. The effective uncertain input representing the calibration uncertainty is H g,m . Finally, the model of equations (1a)-(1c), including its two effective uncertain parameters, is introduced into the multi-physics coupling creating the IBE modeling of Figure 1b. It can be seen as an intermediate modeling between the BE and Best Effort modelings. Since the H gap model is a simple analytic function the computational cost of IBE is similar to BE amounting to a significant gain compared to Best Effort coupling. The average Best Effort computational time for one REA is 3 h while the average BE and IBE is 6 min.

Results
The first step 1.1 consists in clustering assemblies with similar burn-ups. For each cluster one H gap model will be considered. In the PWR core there is a 3D burn-up distribution. This leads to a total of 193 × 4 × 30/2 = 11580 meshes (due to symmetry) with different H gap evolutions due to different burn-ups and power histories. In order to avoid constructing one H gap model for each mesh, the clustering of the assemblies of different burn-ups is carried out. At first we consider only radial burn-ups, by averaging the axial variations. Secondly, we observe that the burn-ups have radially small variations around three main values 15 GWd/t, 30 GWd/t and 45 GWd/t due to the PWR fuel loading pattern. We, therefore, choose to cluster the assemblies based on these three values and we add one cluster for the minimum 10 GWd/t and one cluster for the maximum burn-up 52 GWd/tt. This is done in order to have models covering all the burn-up variations. Additionally, it could be potentially used in the future for an application of a full 3D burn-up dependent H gap model by constructing models that interpolate the calibration parameters.
To summarize we consider a total of 5 fuel assembly clusters and for each cluster one H gap model will be constructed. We have to select for each cluster representative boundary conditions that vary both randomly and spatially since each cluster includes different assemblies. This is performed in step 1.2. The selected boundary conditions are presented in Figure 8. We know that the REA is a local phenomenon located in the upper part of the core around the ejected control rod position as seen in the radial cross-section of the Figure. We thus expect large variations of H gap on the upper part and low to negligible variations in the lower part. For the random aspects we use the results from Section 4. More specifically, we consider that P max lin gives a good indicator of H gap variations in REA. Based on this we select samples corresponding to the mean, the upper 97.5% and lower 2.5% quantiles of P max lin . From these samples we extract the linear power and cladding wall temperature axial and temporal evolutions. We combine both random and spatial aspects by selecting representative assemblies in the upper part for the mean and upper quantiles while we select their mirror assemblies from the lower part. The selected assemblies are presented in Figure 8 where each assembly has the burn-up value of its cluster. The green circles correspond to the selection for the mean and upper quantile while the yellow circles for the lower quantile. For the 10 GWd/t and 52 GWd/t clusters there is only one possible assembly for each cluster. For the other clusters, from the many available options we prefer the assemblies close to the ejected control rod location. For 15 GWd/t and 30 GWd/t clusters we select two different assemblies while for the 45 GWd/t we select three, since we consider that this cluster will have the largest H gap variations due to its high burn-up.
In step 2 the extracted boundary conditions are imposed in ALCYONE1 and the REA stand-alone calculations are performed for each cluster. A representative fuel pin for each selected assembly in each sampling is modeled with the same axial discretization as the BE    Figure 9 and the estimated values of the calibration parameters in Table 4. We observe that the calibration errors increase in general with burn-up. The maximum error is of the order of 8% for the 45 GWd/t cluster.
For the uncertainty quantification of the calibration parameters θ 1 , θ 2 we make the assumption that they are positively fully correlated, which makes it possible to simplify significantly their uncertainty quantification. Additionally, as for all the other inputs they are also fully correlated spatially. This means that the calibration parameters of each cluster vary homogeneously. We assign uniform distribution to each calibration parameter with ranges that cover the calibration errors (the uniform distribution maximizes the entropy amongst distributions with prescribed ranges). The results for the estimated ranges for each parameter are shown in Table 5.
The bounds of the uniform distributions are used for the prediction of the H gap evolution during REA by the calibrated models. The results are compared to the ALCYONE1 calculations and are illustrated in Figures 10,11 and 12. For the comparison the axial slice with the maximum H gap value for the three main clusters is presented. The plotted H gap predictions are also the ones with the largest errors and we can see that for all the predictions the ALCYONE1 calculation is within the uncertainty bounds created by the calibration parameters.
Finally, there are two different sources of uncertainties for the calibrated models. The first one is due to the calibration error and is quantified by the distributions of the calibration parameters. The parameters are assumed fully positively correlated. The second one as discussed previously is due to the initial conditions. At the end, two new uncertain parameters replace the constant H gap of the BE modeling. The first one is related to the calibration error of the models H g,m and the other is related to the initial conditions H g,i .

Improved best estimate uncertainty quantification
The calibrated H gap models from Section 5 are introduced into the BE coupling to create an intermediate IBE multi-physics coupling. The complete UQ methodology detailed in [17] is applied on this improved modeling. It is important to mention that the computational cost does not increase. As in the BE UQ, kriging models were trained on a LHS of size 250 and used in brute force Monte Carlo for uncertainty propagation with 10 5 samples for each output. The results are presented in Figure 13. The validation of the kriging model's training show a Q 2 close to 1 with the smallest being 0.987 for DN B min and 0.945 for the second principal component of P 2D lin . More detailed results are provided in Appendix B.
The obtained histogram for H max f is approximately normal as for the BE UQ. We observe small increase of 2% in the mean value and a reduced relative standard deviation from 20% to 16%. For P max lin the mean value significantly decreases by 10% with similar relative standard deviation. The relative standard deviation of P 2D lin is not affected by the improved H gap model. The output quantity that is the most impacted is the DN B min with an increase of 14% of the mean value and a decrease of the relative standard deviation from 57% to 42%. This means that there is significantly smaller probability to reach boiling crisis. In Appendix B we present also the estimated relative standard deviations of the scalar quantities directly using the training LHS. Although there are some discrepancies due to the limited number of samples, the same trends are observed as well.
The impact on the different mean values is attributed to the more realistic modeling of the H gap evolution during the REA. In the BE modeling the mean constant value of H gap is 2.4e 4 W/m 2K , much larger than the one predicted by the calibrated models, and it is applied during the whole duration of the transient. This leads to more heat transferred from the fuel to the coolant. The fuel temperatures are lower with a corresponding weaker Doppler feedback and thus a higher maximum linear power. The increased heat extracted from the coolant in the BE UQ explains also the smaller minimum distance to boiling crisis compared to the IBE UQ. The lower fuel temperatures induce also the observed lower stored enthalpy.
The global sensitivity analysis is performed by estimating Shapley indices for all the outputs. The results are presented in Figure 14. For P max lin , P 2D lin and H max f as in the BE modeling the β ef f is responsible for around 50% of the outputs variances and the T D 1 and D 1 are responsible for the remaining 50%. A significant difference in observed for the DN B min sensitivities. The gap heat transfer is not any more the dominant input, instead as for the other outputs the β ef f and the T D 1 and D 1 account for most of the DN B min variations. This explains also the significant reduction of the DN B min relative standard deviation.
Finally, we studied the 3D H gap (H 3D gap ) as an output of interest and included it in the uncertainty quantification methodology. For H 3D gap two principal components are needed to represent 95% of its variance. The estimated means and standard deviations for the radial and axial cross-sections at the location and instant of the local maximum are presented in Figures 15,16. As expected by the calibrated H gap models, the assemblies with higher Table 5. Estimated distributions of Hgap model calibration parameters.

Burn-up group:
15 GWd/t 30 GWd/t 45 GWd/t    burn-up have also higher H gap . The maximum value is obtained at the assemblies on the right and left of the assembly with the ejected control rod. This is due to the important power seen by these assemblies in combination with their high burn-up. The relative standard deviation distribution exhibits strong variations from 10% up to 32%. The 10% lower bound on the H gap is related to the initial conditions uncertainties. The aggregate Shapley indices for H 3D gap are estimated and presented in Figure 17. In this case the H g,i is the dominant input responsible for 80% of the outputs variances. The remaining 20% is mainly explained by β ef f . This result is not surprising since the initial conditions determine the H gap evolution.

Conclusion and perspectives
In this work we have presented a High to Low H gap model methodology for REA uncertainty quantification in a PWR core. A simplified H gap model based on fuel thermal expansion has been considered. The CORPUS multi-physics coupling framework is used for the transient modeling. The fuel-thermomechanics code ALCYONE1 is used as the high fidelity code. The available Best Estimate (BE) coupling scheme involves a coupling between core neutronics code APOLLO3 R and core thermal-hydraulics code FLICA4. For APOLLO3 R a two-group diffusion modeling is used while for FLICA4 an axial multi-channel 1D modeling. The internal fuel thermal modeling of  FLICA4 is updated to include the H gap predicted by the calibrated model. This creates an Improved Best Estimate (IBE) coupling scheme. The calibration methodology involves different steps. The starting point is a previous BE REA uncertainty propagation that creates the set of available results from which representative ones will be selected for the calibration. Then the first step is to cluster the assemblies and use one H gap model for each cluster. The assemblies are clustered based on their burn-up, which creates 3 main clusters: 15, 30 and 45 GWd/t. The next step is to select the fuel-thermomechanics boundary conditions on which ALCYONE1 is run for each cluster. REA power pulses are selected from the available BE results in such a way to cover the extreme cases and in locations where the higher and lower power variations are observed. The ALCYONE1 calculations are performed for each cluster power pulses and the H gap and temperatures evolutions are extracted. The H gap model calibration is performed for all the axial slices at two time instances: the instance of the maximum H gap during the REA and the last value of H gap . The calibration parameters are estimated for these datasets and the calibration error for all the time steps is computed. The assumption has been made that the calibration parameters are fully positively correlated. Their uncertainty bounds are selected in order to create H gap intervals that cover the ALCYONE1 results. The results indicate that the uncertainty of the parameters increases with burn-up, which shows the limitation of such a simplified model for high burn-up fuel rods. Additionally, an uncertainty of 10% on the initial rod conditions is considered as a result of ALCYONE1 evolution calculations. Finally, the last step is to include the calibrated H gap model into the BE coupling and to create the IBE.
An uncertainty quantification study is performed using the IBE and compared to the BE results. The most significant impact is on the distance to boiling crisis with an increase of 14% in its mean value reducing the probability of boiling crisis. The relative standard deviation of DN B min is reduced from 57% to 42%. The reason is the decreased sensitivity on the H gap due to its better uncertainty quantification and modeling. This is observed in the estimated Shapley indices, where the H gap sensitivity reduces from 50% to almost 0. These results, although in an academic proof of concept study, show the potential gain of High to Low approaches for H gap modeling in REA uncertainty quantification.
Many future prospects open following this work. On the one hand a more elaborated H gap model could be sought with incorporation of the burn-up parameter to avoid creating one model for each burn-up. Besides that, a more direct High to Low modeling could be investigated in order to inform specific low level code models, such as simplified dynamical H gap models, from the high level code. On the other hand, different transient scenarios and different reactor types could be studied in order to cover a larger spectrum of applications since they all share the common need for High to Low H gap modeling. APOLLO3 R is registered trademark of CEA. We gratefully acknowledge CEA, Framatome, and EDF for their long-term partnership and support. We would like also to thank the APOLLO3 R and FLICA4 development teams for their support and successful exchanges. This work would not be possible with-  have performed scientific advisory roles and provided support through expert judgement, critical verification, and proofreading of the manuscript.

Appendix A: Input uncertainty quantification
The input uncertainty quantification is performed by using the UAM recommendations when they are available and by using expert opinions for the rest. The resulting pdf for the inputs are presented in Figure A.

1.
A multivariate normal distribution is used for the neutronic inputs. The mean vector of the pdf is the reference cross-sections produced at CEA Σ CEA and the covariance matrix C UAM is empirically estimated by the results of UAM. The UAM provides a dataset of two-group macroscopic cross-sections obtained by neutronic lattice uncertainty propagation simulations. This dataset of 100 realizations was adapted to the cross-sections used by APOLLO3 R by assuming: 1. Negligible uncertainties on the up-scattering crosssection S 2→1 . 2. Negligible n − 2n, n − 3n . . . cross-sections uncertainties.
The adapted macroscopic cross-sections dataset is used to estimate the correlation matrix and the relative standard deviations of the cross-sections. Finally, C UAM is calculated using the correlation matrix, the relative standard deviations and the reference CEA cross-sections. The UAM correlation matrix is illustrated in Figure A.1, where we can observe large positive and negative correlations. The T D 1 is strongly positively correlated with S 1→2 and IV 1 and negatively correlated with D 1 . The β ef f is also strongly positively correlated with λ ef f .
The rest of the inputs are considered independent of the neutronic inputs and between them. The thermalhydraulic inputs are mainly based on CEA experts opinions and are applied as multiplicative coefficients on the different parameters. The R crit usually is penalized to 1.3 and in this work we consider this value as the 95% upper quantile of a normal distribution with mean value 1 and standard deviation 0.15. For H DN B a uniform distribution is used reflecting the limited current knowledge about this phenomenon. In fuel-thermomechanics inputs and more specifically for the thermal conductivities and specific heat capacities the UAM recommendations are used. For the Rowlands temperature a uniform distribution is considered on the weight fraction of the fuel centerline temperature. In the reference situation the Rowlands temperature has a 4/9 weight on the fuel centerline temperature and 5/9 on the fuel external surface temperature. By using an uncertain multiplicative factor with U(0, 1) distribution on the centerline temperature weight we consider that it can only decrease uniformly between 4/9 and zero with a corresponding increase in the external surface temperature weight. The H gap is a particular input quantity since in the BE coupling it is an uncertain constant in FLICA4 with uniform pdf bounded from below by its value for a complete open gap (2000 Wm −2 K −1 ) and from above by its value for a pelletcladding contact (50000 Wm −2 K −1 ). In the IBE coupling H gap is a result and thus is not considered, instead the calibration parameter's uncertainty and the initial conditions of the model are used. Finally, the power radial profile uncertainty is modeled by a multiplicative factor on the fuel external surface power with a pdf resulting from the convolution of a normal and a uniform distribution. The power radial profile is peaked towards the periphery for the high burn-up fuel pins. The deformation is increasing with burn-up and to model it an explicit function of burn-up is used. The uncertainty is a convolution of two effects: the uncertainty of the function used N (0, 0.0175) and the uncertainty due to the presence or not of a guide tube near the fuel pin U(0, 0.04). The latter effect is a result of a previous CEA study. independent LHS of size 125. For all the outputs the Q 2 is close to 1 with smallest values for DN B min (0.987) and the second principal component of P 2D lin (0.945). Finally, we present in Table B.5 the estimated relative standard deviations using directly the training LHS code evaluations. We observe some differences compared to the kriging model's brute force Monte Carlo due to the limited number of samples. However, similar trends are observed T D1, N F2, D1,S1→2, β ef f , Hgap P 2D lin,pc2 T D1, N F2, D1,S1→2, β ef f , Hgap H max f T D1, N F2, D1, β ef f , Cp f , Hgap,TR DNB min T D1, N F2, D1,S1→2,IV1, β ef f , Hgap, Hc