Lognormal-Based Sampling for Fission Product Yields Uncertainty Propagation in Pebble-Bed HTGR

Uncertainty analyses of fission product yields are indispensable in evaluating reactor burnup and decay heat calculation credibility. Compared with neutron cross section, there are fewer uncertainty analyses conducted and it has been a controversial topic by lack of properly estimated covariance matrix as well as adequate sampling method. Specifically, the conventional normal-based sampling method in sampling large uncertainty independent fission yields could inevitably generate nonphysical negative samples. Cutting off these samples would introduce bias into uncertainty results. Here, we evaluate thermal neutron-induced U-235 independent fission yields covariance matrix by the Bayesian updating method, and then we use lognormal-based sampling method to overcome the negative fission yields samples issue. Fission yields uncertainty contribution to effective multiplication factor and several fission products’ atomic densities at equilibrium core of pebble-bed HTGR are quantified and investigated. )e results show that the lognormal-based sampling method could prevent generating negative yields samples and maintain the skewness of fission yields distribution. Compared with the zero cut-off normal-based sampling method, the lognormal-based sampling method evaluates the uncertainty of effective multiplication factor and atomic densities are larger. )is implies that zero cut-off effect in the normal-based sampling method would underestimate the fission yields uncertainty contribution. )erefore, adopting the lognormalbased sampling method is crucial for providing reliable uncertainty analysis results in fission product yields uncertainty analysis.


Introduction
Reactor design and safety analysis rely on accurate calculations of system responses with properly evaluated uncertainties. ere has been an increasing need for evaluating the credibility of reactor safety. Pebble-bed high temperature gas-cooled reactor (pebble-bed HTGR) is a multiphysics nonlinear coupled system, including neutron transport and complex heat transfer hydraulics behaviour [1]. In order to systematically and thoroughly investigate the uncertainties propagation in pebble-bed HTGR, an IAEA Coordinated Research Plan (CRP) [2,3] has been initiated after the start of OECD/NEA UAM-LWR [4]. Recent advances regarding uncertainty propagation analysis in pebble-bed HTGR mainly concern the nuclear cross section uncertainties propagation in reactor neutronic calculations [5][6][7][8][9]. As pebble-bed HTGR allows fuels recirculation during fuel cycles and adopts higher fuel enrichment (8.5 wt.%), fuels usually could achieve larger burnup values, and then fission product yields could be nonnegligible uncertainty sources in reactor burnup and decay heat calculations. eir uncertainties contributions to important reactor burnup responses need to be considered properly for evaluating the credibility of rector safety-related quantities of interest (QoI), e.g., maximum fuel pebble temperature.
Fission product yield describes the fraction of a certain fission product produced per fission. During the measurements of fission product yields, correlated errors or covariances may exist when using the same equipment or methods [10]. However, they are ignored in evaluated nuclear data library. Also, self-consistent fission yields data set should follow several physical constraints such as binary fission, mass conservation, and charge conservation [11]. ese constraints could introduce covariances between fission yields data. As the fission yields' covariances in current releases of evaluated nuclear data libraries are still absent, e.g., ENDF/B-VII.1, many methods are developed to estimate these covariances information based on the imposed physical constraints. e Bayesian updating method is widely used in data assimilation, data adjustment, and model fitting problems. It refines parameters by taking both the prior information about those parameters and the likelihood which refers to new data into consideration [10,12]. It allows estimating the covariance matrix of fission yields by sequentially introducing the above physical constraints. In the domain of fission yields adjustment, it is introduced by Kawano and Chadwick [13] to update Pu-239 fission yields with chain yields to reduce the independent fission yields discrepancy in ENDF/B-VII.1. Pigni et al. [14] expand it to involve cumulative fission yields into covariances estimation. e difference between chain yields-based updating and cumulative fission yields-based updating is further investigated by Fiorito et al. [15,16]. Based on the provided independent and cumulative fission yields uncertainties information in ENDF/B-VII.1, this work adopts the Bayesian updating method to estimate the independent fission yields covariances.
Sampling-based methods for uncertainty analysis or stochastic UQ methods [17] require properly perturbed samples to provide reliable uncertainty analysis results of QoI. As it is observed from the evaluated nuclear data library, independent fission yields generally have larger uncertainties. Random sampling on these yields under normal distribution could generate nonphysical negative samples. Cutting off these negative yield samples by setting them to zero could artificially omit some information from the original fission yields distribution, resulting in biased uncertainty analysis results. is zero cut-off effect on quantified uncertainty has not been well studied.
is raises question whether normal distribution is sufficient for describing inherently positive random variables with large uncertainties only given their mean values and covariance matrix. Smith et al. [18] propose to replace normal distribution with lognormal distribution by the principle of maximum entropy [10], andŽerovnik et al. [19,20] investigate this method in the sampling resonance parameters where negative samples problem was encountered as in fission yields. is work proposes an implementation of the lognormal-based sampling method in fission product yields sampling.
e present work is organized as follows: Section 2 describes the nomenclature of fission product yields and the burnup calculation of pebble-bed HTGR. An implemented stochastic UQ method for fission yields uncertainty propagation is described in Section 2.3. e Bayesian updating method and the lognormal-based sampling method are detailed in Section 3. Finally, results of fission yields uncertainty contributions to effective multiplication factor and several important fission products atomic densities are provided and discussed in Section 4.

ENDF/B-VII.1 Fission Product Yields Sublibrary.
Fission product yield characterizes the fraction of a particular fission product nuclide produced per fission. A compound nucleus is formed when a fissile nucleus is bombarded by an incident neutron. As its energy overcomes the fission barrier, this compound nucleus could undergo fission. A brief description of the fission process is illustrated (see Figure 1). After the scission of compound nucleus, primary fission fragments are produced and they would undergo deexcitation by releasing prompt neutrons due to their high neutron to proton ratios. After the emission of prompt neutrons, the remaining fission fragments referred to as fission products would undergo β decay, isomeric transition, or particle emission along their corresponding decay chain and finally reach stable nuclides. Each fission product is identified by its mass number A, charge number Z, and isomeric state I and is denoted as the triplet (A, Z, I) A detailed description about the nomenclature of fission product yields could be found in [11] and they are briefly summarized as follows. IFYs and CFYs determine the fraction of a fission product at different stages in the fission process. IFYs denoted as y(A, Z, I) are the fraction of a fission product produced directly from one fission after the emission of prompt neutrons but prior to any radioactive decays. Because IFYs are produced before any radioactive decay in the fission system, they should be subject to the physical constraints of fission system, e.g., binary fission, conservation of mass, and charge number. CFYs denoted as c(A, Z, I) determine the total fraction of a fission product produced over all time after one fission. It involves not only the direct production from fission but also the contributions from the decay of other products. e current releases of ENDF/B-VII.1 fission yield sublibrary provide fission yields data for 31 fission actinides from -227 to Fm-255. ough energy-dependence issues within fission spectrum are highlighted in current releases of evaluated nuclear data library and neutron induced Pu-239 fission yields at 2.0 MeV are supplemented to allow users to linearly interpolate yields between 0.5 MeV and 2.0 MeV for high accuracy purpose [22], other fission actinide fission yields data are taken directly from ENDF/B-VI evaluated by England and Rider [23] in 1993. ree fission systems for U-235 are evaluated with respect to incident neutron energy, namely, 0.0253 eV thermal energy, 0.5 MeV fission spectrum energy, and 14.0 MeV high energy. IFYs and CFYs are evaluated for 1,247 fission products in thermal neutron induced U-235 fission yield (see Figure 2). e relationship between IFYs and CFYs [11] is referred to as (1), where b(A′, Z ′ , I ′ ⟶ A, Z, I) is the branching ratio: (1) It could be found that most IFYs appear in the upper region of β-stability line and they are most likely to undergo β − decay to reach a stable state. As CFYs involve the production of a certain fission product from the decay of other fission products as shown in (1), the peaks of CFYs distribution in neutron-charge number figure tend to be closer to the β-stability line (see Figure 2). e evaluation of fission yields data requires a combined work of experimental measurements and theoretical model predictions. It is natural for the evaluated fission yields possessing uncertainties originated from measurement errors and theoretical model parameters uncertainties. Although England and Rider provide the uncertainties (standard deviation) of each fission product yield in their original work, covariances information between fission yields has not been provided since then. ose covariances information is crucial for representing the physical constraints imposed on IFYs, and they should be estimated properly in order to generate self-consistent IFYs  Figure 1: Neutron induced fission process [21]. e fission products refer to the fission fragments after the emission of prompt neutrons. Independent fission yields (IFYs) characterize the fraction of a fission product produced before any radioactive decay, whereas cumulative fission yields (CFYs) describe the fraction of that produced product over all time after a fission. perturbations.
is work focuses on the propagation of thermal neutron induced U-235 fission yields uncertainties to burnup simulation of pebble-bed HTGR based on ENDF/ B-VII.1. e estimation of covariances information will be detailed in Section 3.1.

Pebble-Bed HTGR Burnup Model and Built-In Fission
Yields Analysis. Pebble-bed HTGR core (see Figure 3(a)) consists of spherical fuel elements or fuel pebbles. Each of these pebbles is composed of a spherical graphite matrix in the centre where thousands of small coated particles (known as TRISO particles) are embedded. ese particles contain UO2 kernel in the centre with four structural coating layers surrounding it (see Figure 3(b)). During reactor operation, these fuel pebbles are consistently flowing downward from the top of the core to the bottom and are irradiated at different core spectrum regions randomly. Fuel recirculation is a characterized fuel cycling procedure adopted in pebblebed HTGR which is different with that applied in Light Water Reactor (LWR). Such recirculation allows fresh fuel pebbles being loaded into the core and spent fuel pebbles being discharged online without shutting reactor down. More importantly, this recirculation permits fuel pebbles running through core multiple times before they are finally being discharged. Because of the fuels recirculation, there exist running-in phase and equilibrium core states. e equilibrium core state refers to the nuclei compositions inside the core kept unchanged with time and therefore effective multiplication factor being stable at a certain value. is could give a more flattened power distribution across the core and higher average discharge burnup value. e V.S.O.P. computer code system [26] is developed to perform burnup calculation of pebble-bed HTGR by simulating the fuels recirculation process stepwise and conduct spectrum calculation online at each spectrum region inside the core. A detailed description of this simulation process could be found in these articles [27]. e built-in fission product chain in V.S.O.P. code involves 44 fission products and among these 44 fission products' data, 14 are taken as IFYs, while 30 are taken as CFYs. ese data are taken from ENDF/B-IV and ENDF/B-V. An additional "nonsaturating" fission product is evaluated to account for the sum of many lumped fission yields which are not explicitly included in the chain [26]. e comparison between these built-in fission yields and replacing them with the current releases in ENDF/B-VII.1 is conducted to examine the availability of V.S.O.P. code for fission yields uncertainty propagation. is investigation is conducted on HTR-PM [28] with 8.7% fuel enrichment (while 8.5% enrichment is applied in actual design). e impact of each built-in fission yield on k eff at equilibrium core state is investigated individually by replacing them with ENDF/B-VII.1. It should be noted that built-in fission yields library in V.S.O.P. includes a combination of IFYs and CFYs, and they are presented separately (see Tables 1 and 2). IFYs are evaluated by subtracting the total contributions of its precursors from experimental measured CFYs. With the improvement of CFYs measurements, the evaluated IFYs become more precise. It could be seen from the table that IFYs in ENDF/B-VII.1 are lower than the built-in fission yields used in V.S.O.P. Except the large discrepancy in the fission yield of Mo-95, all the impacts from replacing fission yields are lower than 20 pcm. e overall impact is 67 pcm (see Table 3) when all the yields are replaced without FP-44. e difference is acceptable in effective multiplication factor calculations when substituting built-in V.S.O.P. fission yields with ENDF/B-VII.1 fission yields. e V.S.O.P. burnup model is further used to conduct fission product yields uncertainty propagation as described in Section 2.3.

Uncertainty Quantification Scheme.
e HTR-PM [28] reactor core is modelled in V.S.O.P. computer code system to analyse the uncertainty propagation of fission yields in this work. 15 times recirculation of fuel is adopted and the average discharge burnup value is around 90, 210 MW · d/tU with fresh fuel having 8.5 wt.% enrichment. As fuel recirculation tightly couples the neutronics and burnup calculation spatially inside the core, it is difficult to separate the uncertainty propagation step by step. Stochastic UQ method is used to investigate the uncertainty propagation in equilibrium core state. An uncertainty propagation scheme is proposed in this work (see Figure 4).
Two sampling methods are implemented in this work, namely, normal-based sampling and lognormal-based sampling. Different from normal-based sampling, lognormal-based sampling requires a lognormal transformation of the original mean vector and covariance matrix. When the IFYs samples are generated, their corresponding CFYs are calculated and combine them to form self-consistent fission yield samples. ese combined IFYs and CFYs samples are propagated to V.S.O.P. HTR-PM model for further uncertainty analysis. Detailed Bayesian updating method description and lognormal-based sampling procedures will be introduced in Section 3.

U-235
ermal Neutron-Induced IFYs Covariances Estimation. Bayesian updating method or the generalized least square method (GLSM) is a data adjustment method, which allows the prior data being updated with combination of new knowledge about these data. Such knowledge could be measured data or physical constraints imposed on these prior data. e present work applies Bayesian updating method to estimate the covariance matrix of IFYs based on ENDF/B-VII.1 thermal neutron induced U-235 fission yields sublibrary. e specification of this method is briefly recalled as follows.
Consider a multivariate linear regression model shown in where c and y ∈ R n×1 are observables and parameters to be updated or estimated, respectively. X ∈ R n×n is the design  Figure 3: Pebble-bed HTGR core. (a) Core geometrical of PBR250 design [24]. (b) Fuel pebbles [25].
Science and Technology of Nuclear Installations matrix that represents linear mapping between estimating parameters and observables. ε ∈ R n×1 are the measurement errors of observables with expectation E[ε] � 0 ∈ R n×1 and variance Var[ε] � V ∈ R n×n . By the principle of maximum information entropy, it is objective and plausible to assign Gaussian distribution on this error. Similarly, estimating parameters y could also be assigned Gaussian distribution given their expectation E[y] � y 0 and variance Var[y] � Z 0 . e generalized least square problem [29] is formulated by the following minimization in the domain of estimating    parameters to find the best least square estimated parameters as e above minimization process could also be interpreted in the perspective of Bayesian updating. Consider the estimated parameters have a prior of Gaussian distribution with density p(y) in And likelihood function determines the probability of any candidate estimated parameters appearing in the observables distribution. en, likelihood function p(c | y | ) is given as e posterior distribution of estimated parameters y is therefore calculated by Bayesian theorem and it gives Considering the conjugacy between Gaussian prior and likelihood, the posterior estimated parameters follows Gaussian distribution as well. Under quadratic loss, the optimal estimates of true values and their uncertainty are the mean vector and covariance matrix of the posterior distribution. It is worthwhile to mention that the estimated mean vector could maximize the exponential term in (3) and this could also lead to the solution of GLSM in (3). e posterior estimated parameters are obtained as where Z 1 is the posterior covariance matrix of estimated parameters and it is shown in (8), and after applying Woodbury matrix identity, it is reformed as (9): Here, regarding IFYs as estimated parameters y with prior covariance matrix Z 0 (diagonal matrix with only consideration of each fission yields uncertainty in ENDF/B-VII.1), observables c represent the evaluated CFYs in ENDF/B-VII.1, total independent yields, fission system total mass number, and charge number, respectively. e corresponding design matrix could be formulated as follows: (1) Consistency with CFYs: c � My, where M is the Qmatrix proposed in [11]. It could be formulated from the linear mapping in (1) with the provided branching ratios data in ENDF/B-VII.1 decay sublibrary. is updating process follows Luca Fiorito's updating procedures [15] on CFYs consistency in JEFF-3.1.2. Different than in previous work [14], this work explicitly constructs this design matrix with branching ratios rather than obtaining each element via direct perturbations using a burnup code. Such procedures allow direct examination of consistency between IFYs and CFYs in the current releases of ENDF/B-VII.1. Total IFYs, total mass number, and total charge number conservations are implemented following the procedures proposed in Pigni et al.'s work [14]. e updating results of IFYs' covariance matrix are in (2) Conservation of binary fission: T y � U T y, where U ∈ R n×1 is a unity vector. e sum of total yield T y is 2.0 with summation precision of σ 2 sum � 1.0 × 10 −5 . e updated covariance matrix subsequent to (10) is listed in (11). It should be noticed that ternary fissions may occur; however they are not considered in ENDF/B-VII.1 and these ternary fissions are not included in this updating process: (3) Conservation of fission system mass number: where N ∈ R n×1 whose element corresponds to the mass number of each fission product. e total mass number of fission system is conserved to 233.57915 (considering the average prompt neutrons released at 0.0253 eV is 2.42085 recorded in ENDF/B-VII.1 and mass defect of U-235 is not considered). e assumed variance of total mass number is 1.0 × 10 −5 . e updated covariance matrix subsequent to (11) is shown in (4) Conservation of fission system charge number: with each element being the charge number of each fission product considered. e total charge number of fission system is conserved as 92.05318. is total charge number is calculated from the charge numbers of each fission product weighted by their corresponding IFYs provided in ENDF/B-VII.1. It is observed in this work that if we take the total charge number as exactly 92.0, the calculated CFYs calculated from updated IFYs will have large discrepancy with CFYs provided in the library. And this discrepancy will be narrowed when we take the decimal digits into consideration. e updated covariance subsequent to (12) is shown in Correlation matrix of updated IFYs is plotted (see Figure 5). ese correlations are introduced sequentially to cooperate the consistency with CFYs, conservation of binary fission, mass number, and charge number of fission system. Figure 5(a) shows that there is a significantly two-humped tendency in the correlation distribution. is tendency is similar with the two-humped distribution of IFYs, where many correlations are introduced from the conservation constraints in fission system, while fewer correlations are introduced between humped part and valley part. And Figure 5(b) presents a close look of the correlations among fission product index range from 65 to 245. It could be noticed that the diagonal of this correlation matrix is divided into several small groups regarding different decay chains. IFYs within each decay chain have negative correlation with each other introduced from the consistency of CFYs. e updated IFYs are compared with the prior fission yields recorded in ENDF/B-VII.1 (see Figure 6). It could be seen that small adjustment is introduced to fission product yields in the two-humped part, while larger adjustment is introduced in the valley and two tail parts. is is mainly because IFYs in those parts have smaller prior fission yields and they are not as accurately evaluated as those larger ones in the two-humped part. erefore, more adjustments are expected in those regions. e updated and prior standard deviations are presented and compared (see Figure 7). It could be seen that the adopted updating procedures could reduce the uncertainty of updated IFYs. is is mainly due to the introduced constraints that further constrain the uncertainty of these fission yields and introduce covariances among them. e final updated covariance matrix of IFYs Z 4 and the posterior IFYs mean vector y 4 are applied to generate the perturbation samples of IFYs. e detailed sampling procedures are further discussed in the following section.

Lognormal-Based Sampling Procedures. Considering
IFYs are inherently positive, random sampling under normal distribution could draw unphysical negative samples. ese negative samples would appear significantly when the sampled parameters have large uncertainty (e.g., relative difference σ/μ > 30%). Smith et al. concluded that when the relative uncertainty of a random variable exceeds 30%, the probability distribution of this parameter chosen to represent its physical uncertain information tends to be skewed noticeably [18] and the drawn negative samples fraction tends to grow. It could therefore be concluded that normal distribution is not adequate to describe inherently positive random variables whose uncertainties are large, because it could not capture the skewness of random variable distribution. By the principle of maximum information entropy, lognormal distribution is suggested to be the optimal choice for inherently positive parameter when only expectation and variance are known about this parameter [10,29]. Larger relative uncertainty would result in a more skewed distribution (shown in Figure 8). Lognormal distribution is shifting to a normal-like distribution as its relative uncertainty becomes lower than 30%, where skewness of the distribution is not significant. e updated posterior IFYs relative uncertainties are compared with prior relative uncertainties (see Figure 9) in our previous work [30]. Except for a few fission products which have their relative uncertainties increased, most fission products have their corresponding relative uncertainties decreased to around 42%. e increased relative uncertainty fission products are Ag130m0, Cd129m0, Sn127m1, Cd126m0, In126m0, Sb124m1, Zn123m1, Ag115m0, Y93m1, Y93m0, Se85m1, and Ge77m0. eir relative uncertainties increased due to their updated smaller posterior mean values. From Figure 10, it could be seen that most fission yields standard deviations have been reduced because of the updating process. However, the above fission products have their mean value updated even smaller and that makes their relative uncertainties increased. Compared with the listed monitor fission products for fission of U-235 in Fiorito et al.'s work [15], they are not included and we may think they are less relevant to the reactor burnup and critical calculation. When applying simple random sampling procedures under normal distribution, drawing samples in R n×S from the N(y 4 , Z 4 ), where n is the number of fission yields and S is the sample size, it is almost impossible to draw a sample set with all positive yields as the yields domain is too large (e.g., n > 900).
In this work, lognormal random sampling procedures are applied to generate IFYs perturbation samples. e sampling follows the development ofŽerovnik et al. [19] and applies it into the generation of IFYs samples. Multivariate lognormal distribution is defined as where y is the posterior IFYs with expectation y 4 and covariance matrix Z 4 estimated by Bayesian updating method discussed in Section 3.2 and L ∈ R n×1 is the natural logarithmic value of independent yields. μ l and Z l are the corresponding mean and covariance matrix of IFYs in the natural logarithmic domain. e detailed derivation of their relation with parameters in original domain (y 4 and Z 4 ) could be found in [20]. e basic idea is recapped in the following. Consider the preservation of probability; the relation between random variables in original domain and logarithmic domain is formulated in e lognormal distribution density is therefore derived as in   Before update A er update    Figure 8: Lognormal distribution of random variable X in terms of its relative uncertainty. Relative uncertainty R � (σ X /μ X ) is ranged from 10% to 80% and μ X � 2.0. Dashed line shows the distribution with relative uncertainty lower than or equal to 40%, whereas solid line indicates the distribution with relative uncertainty larger than 40%.   [30]. e increased relative uncertainty fission products are Ag130m0, Cd129m0, Sn127m1, Cd126m0, In126m0, Sb124m1, Zn123m1, Ag115m0, Y93m1, Y93m0, Se85m1, and Ge77m0.
With the logarithmic density function, each element in μ l and Z l is derived as cov ln x i , ln x j � ln cov y i , y j where cov(y i , y j ) and μ[y i ] are retrieved from the posterior updated IFYs covariance matrix Z 4 and updated IFYs mean vector y 4 . With the calculated distribution parameters μ l and Z l , the lognormal-based IFYs sampling procedures could be conducted as follows: (1) Obtain prior IFYs information including IFYs value y 0 as well as its covariance matrix Z 0 from ENDF/B-VII.1 fission yield sublibrary. Implement Bayesian updating procedures detailed in Section 3.2 on the prior information and obtaining the updated IFYs mean vector y 4 and the estimated covariance matrix Z 4 . (2) Consider IFYs follow lognormal distribution, and transform y 4 and Z 4 into natural logarithmic domain with (17) and (18). e normal distribution parameters of natural logarithmic yields are obtained as mean vector μ l and covariance Z l .
(3) e transformed logarithmic domain covariance could not remain symmetric positive definite (SPD) due to the numerical error in the transformation procedure. A nearest-SPD searching algorithm [31] is therefore applied to search for the nearest SPD approximation of the calculated covariance matrix in the sense of least Frobenius norm difference. e approximated SPD logarithmic domain covariance matrix is thus obtained as Z l ′ . (4) Implement the simple random sampling procedures in the logarithmic yield domain with distribution parameters mean μ l and approximated SPD covariance matrix Z l ′ . And the generated logarithmic fission yields sample matrix P n×S is obtained where n denoted the number of fission products considered and S is sample size. (5) Take the exponential transformation of each element in sample matrix P n×s and the sampled negative-free samples are generated and denoted as Y n×S .
e nearest-SPD searching algorithm approximates non-SPD covariance matrix Z l by an approximated matrix Z l ′ with relative difference in Frobenius norm (‖Z l − Z l ′ ‖ F /‖Z l ‖) � 9.74% and their corresponding eigenvalues distributions are presented in Figure 11. e nearest-SPD searching algorithm could approximate a non-SPD covariance matrix while most of its eigenvalue unchanged. e approximation that resides in the above sampling procedures is the SPD approximation of calculated covariance matrix.
is approximation could affect consistency of each drawn IFYs sample with the physical constraints imposed on it.
ere are 1,000 IFYs samples drawn with the lognormal sampling procedures. And the sample mean and standard deviation (STD) for each fission product yield and Pearson's correlation coefficient between these fission yields are calculated and justified by comparison with its corresponding population values in updated y 4 and Z 4 (see Table 4). Table 4 indicates that the proposed lognormal sampling procedures could obtain an overall representation of IFYs population distribution considering the lower RMSE. However, there still exist a few fission products listed in Figure 12 having large biases compared with their corresponding population values considering the maximum of absolute relative difference. After comparing these fission products with the monitor fission products for thermal neutron induced U-235 fission listed in Fiorito et al.'s work [15], they are not included and could be considered less relevant to reactor burnup and criticality calculations. ese outliers' appearance could result from the nearest-SPD procedures and a further investigation regarding this will be conducted in future work. Figure 13 presents the sampled Pearson's correlation coefficients relative difference to their corresponding population values. It could be seen that simple random sampling procedure is not an efficient sampler for sampling low correlation fission yields (|ρ| < 0.1) as shown in the neighbour around 0.00 in this figure. However, these low correlations could have little impact on the uncertainty quantification of fission yields compared with large correlations. As for the larger correlations (|ρ| > 0.25), 1,000 samples are sufficient for maintaining the Bayesian updated correlations and this discrepancy could be further reduced when increasing the sample size. A more efficient sampler, like Latin Hypercubic Sampler (LHS), could be adopted to guarantee more precise results when using 1,000 samples and this will be adopted in future work. e consistency of IFYs samples with these imposed physical constraints is justified in Table 5. e conservation parameters (e.g., total fission yields, total mass number, and total charge number) are calculated for each yield samples and the mean and standard deviation are summarized to compare with the target conservation value. It is found that although the consistency is not strictly restored as the standard deviation of the total yield is larger than the imposed 10 − 5 , their mean values are close enough to the target value indicating the constraints are maintained. e large standard deviation is originated from the approximation mentioned above. In order to examine the difference between normalbased sampling and lognormal-based sampling, 1,000 samples are drawn from the IFYs distribution of Zr95m0, Mo95m0, and Cs134m0. Notation m0 indicates these fission products are at ground state. e IFYs of these three fission products are explicitly involved in V.S.O.P. burnup calculation and are important for reactor decay heat release calculations. Especially for Cs134m0, it is one of the main decay heat contributors of UOX fuels in long-term after reactor shutdown [32]. e updated relative uncertainty of Zr95m0 IFY is 16.1% while Cs134m0 and Mo95m0 have their relative uncertainties of 38.4% and 65.7%, respectively. From the sampled histogram of these fission products IFYs samples (see Figures 14-16), lognormalbased sampling procedures (blue bars) could effectively capture the skewness of these fission yields and permit "negative-free" samples. It is also observed that the skewness of these fission products would become larger as their relative uncertainties become larger (e.g., Mo95m0 and Cs134m0).

Uncertainty Analysis of the Effective Multiplication Factor at Equilibrium Core.
e unperturbed burnup calculation is conducted with V.S.O.P. built-in fission yields library and ENDF/B-VII.1 posterior fission yields. Figure 17 shows that reactor achieved the equilibrium state after operating longer than 2500 days. Effective multiplication factor calculated from ENDF/B-VII.1 posterior fission yields is compared with that calculated from V.S.O.P. built-in fission yields and the total discrepancy at equilibrium core state (which is at the end point of fuel cycle time in Figure 17) is lower than 50 pcm which is small enough for the following fission product yields uncertainty propagation analysis. e comparison between ENDF/B-VII.1 posterior fission yields predicated k eff (black dashed line) and builtin yields predicted k eff (orange dashed line) are shown in Figure 18 is discrepancy is within the sampling distribution of k eff .
1,000 fission yields samples are generated with normalbased sampling procedures and lognormal-based sampling procedures and they are propagated to V.S.O.P. burnup calculation to obtain k eff samples under equilibrium core state (3049 days). e sample distributions from these two sampling procedures are drawn and compared (see Figure 18). It is obvious from the comparison that normalbased samples contain fewer distribution information compared with lognormal samples as its distribution range is smaller than that in lognormal samples. is is due to the zero cut-off procedure of the uncontrolled negative samples. Such procedure artificially omits certain information in the original fission yields distributions and could not provide a Relative difference (%) Original log-domain relative covariance matrix Searched log-domain relative covariance matrix  correspondingly reasonable and satisfied sampling distribution of k eff . In this sense, lognormal sampling procedures overcome this problem by imposing a more plausible distribution on fission yields and allow the generation of smaller perturbed samples. erefore, it leads to a negative skewness (long tail in left) of effective multiplication factor distribution and permits a more rational and persuasive sampling distribution. e uncertainty analysis results are presented (see Table 6). e propagated sampled distribution of k eff from normal-based sampling method passes the normality test with p value 0.3737 and the quantified relative uncertainty is around 1.09 × 10 − 4 . Lognormal samples provide a skewed k eff distribution and fails the normality test with p value smaller than 0.05. e quantified relative uncertainty from  e k eff quantified from lognormal-based sampling method is larger than that from normal-based sampling method, and this shows that the zero cut-off effect in normal-based sampling method could cause underestimation of fission product yields uncertainty contribution to QoIs.

Uncertainty Analysis of Certain Fission Products Atomic
Densities. In this section, fission products Zr95m0, Mo95mo, and Cs134m0 atomic densities uncertainties contributed from fission products yields are quantified. Specifically, their uncertainties differences from different sampling methods are compared and discussed. From the  e horizontal axis indicates the average burnup values of these fuels. As HTR-PM allows recirculation of fuels, 15 times recirculation is adopted in this analysis, which indicates these fresh fuels will be reloaded into the core 15 times before they are finally discharged. e discharged burnup value or the end point of the horizontal axis is 90210.44 MW·d/tU. roughout the burnup process, the thermal power of reactor core is kept at 250 MW. e atomic densities of Zr-95m0 fluctuate along with the increases of fuels burnup value. is fluctuation is due to the fuel recirculation procedures adopted in V.S.O.P. burnup calculations. ere are total 14 lower valleys that appeared in dashed line of Figure 19 which corresponds to the 14 times reloading of the fuels from the bottom of the core to the top. For each reloading, the fuels will be irradiated again during their passes through the core. As it could be seen from   solid line). e blue and orange shadings in these figures are the 95% confidence interval of relative uncertainty computed by bootstrap method. From these figures, it is worth to mention that lognormal-based sampling quantified atomic density relative uncertainties are larger than that quantified from normal-based sampling for all of these three fission products. is is reasonable as zero cut-off adopted in normal-based sampling method would artificially omit some information provided by fission yields distributions, and this would result in an underestimated atomic density relative uncertainty quantification result. After closely comparing the atomic density relative uncertainties underestimation for Zr95m0 and Cs34m0, it could be seen that this underestimation effect will be enlarged when the fission products IFYs have larger relative uncertainties (Zr95m0: 16.1% and Cs134m0: 38.4%).
is is because lognormal distribution would resemble normal distribution when the random variate has smaller relative uncertainty, as discussed in Section 3.2. And in this case, lognormal-based sampling results would be in agreement with those calculated from normal-based sampling. erefore, this underestimation would be narrowed.
Besides, another interesting phenomenon is observed here. is underestimation seems not positively correlated with the relative uncertainty of random variates, as it is seen from comparison between Mo95m0 and Cs134m0. Although Mo95m0 has its relative uncertainty (65.7%) larger than Cs134m0 (38.4%), the underestimation effect observed from Figures 20 and 21 shows that the underestimation effect of Mo95m0 is smaller than that of Cs134m0. One possible reason could be the decay of these fission products. As Mo95m0 is the direct descendant of Zr95m0 whose half-life is around 64 days, its atomic density relative uncertainty is contributed both from its own fission yields uncertainty and the atomic density uncertainty of Zr95m0. As Zr95m0 atomic density uncertainty is less underestimated, the atomic density relative uncertainty underestimation in Mo95m0 is therefore counterbalanced. While Cs134m0 is treated as stable fission products in V.S.O.P. burnup fission product chains, its atomic density relative uncertainty is directly related to its fission yields uncertainty and large  atomic density relative uncertainty underestimation could be seen. e atomic density relative uncertainties of all these three fission products quantified at 90210.44 MW·d/tU are summarized in Table 7.

Conclusions
e present work proposed a stochastic UQ method for propagation fission products yields uncertainties. V.S.O.P. code [26] is used to conduct the burnup calculation of HTR-PM reactor core with allowing 15 times recirculation of fuel pebbles [30]. Uncertainties of thermal neutron induced U-235 IFYs are investigated in this work based on ENDF/B-VII.1. Bayesian updating method is applied to estimate the covariance matrix of IFYs. Lognormal-based sampling method is implemented to generate perturbations of yields samples. e differences of quantified uncertainties between conventional normal-based sampling method and lognormal-based method are addressed and investigated. Specifically, the effect of zero cut-off procedures used in normalbased sampling method is studied and discussed. From the above investigation, conclusions are summarized as follows: (1) Lognormal-based sampling method could effectively overcome the negative samples generation caused by the large relative uncertainties in fission yields data. Compared with normal-based sampling method, it could provide reasonable and negative-free fission yields samples to permit a more plausible and reasonable QoI sampling distribution for further uncertainty analysis.
(2) e contribution of thermal neutron induced U-235 fission yields uncertainties in ENDF/B-VII.1 to k eff of pebble-bed HTGR at equilibrium core is 0.0258%. is contribution is smaller than that from neutron cross section 0.48% at equilibrium core [33].
(3) e zero cut-off procedures used in conventional normal-based sampling method to overcome the negative fission yields samples appearance would underestimate the uncertainty analysis results. For relative uncertainty of effective multiplication factor, it would underestimate the results by 0.0149% which is around 42% of results obtained from lognormal-based sampling method. For atomic density relative uncertainty, the underestimations are also observed, and especially for Cs134m0, this zero cut-off effect would underestimate the atomic density relative uncertainty by 22% compared with lognormal-based quantified results.
It is worth to mention that there are several approximations and simplifications made during the Bayesian updating process and implementing of lognormal-based sampling methods in this work. e considered constraints for Bayesian updating independent yields covariance matrix are preliminary in this work and a more complete and comprehensive study regarding this will be conducted in future work. Also, the effect of using nearest SPD algorithm in implementing lognormal-based sampling method will be investigated in the future. For the following work, additional fission systems will be investigated with the proposed uncertainty propagation scheme. And a sensitivity analysis of effective multiplication factor to fission yields should be conducted to determine the reason behind the formation of effective multiplication factor skewed distribution.

Nomenclature
IFYs or y(A, Z, I): Independent fission yields CFYs or c(A, Z, I): Cumulative fission yields A: Nuclide mass number Z: Nuclide charge number I: Nuclide isomeric state b(A ′ , Z ′ , I ′ ⟶ A, Z, I): Branching ratio k eff : Effective multiplication factor y 4 : Bayesian updated IFYs mean vector Z 4 : Bayesian updated IFYs covariance matrix μ l : Natural logarithmic value of IFYs mean vector Z l : Natural logarithmic value of IFYs covariance matrix Z l ′ : Nearest-SPD approximated Z l SPD: Symmetric positive definite μ: Mean σ: Standard deviation ρ: Pearson's correlation coefficient R n×1 : n-dimension real vector R n×n : n-dimension real matrix.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.