Uncertainty Propagation in SINBAD Fusion Benchmarks with Total Monte Carlo and Imprecise Probabilities

Abstract In this paper we perform nuclear data uncertain propagation with Total Monte Carlo, where the transport simulation is repeated for random evaluations of the data. The Oktavian Iron, Oktavian Nickel, and the Frascati Neutron Generator (FNG) neutron streaming SINBAD benchmarks were evaluated with OpenMC. Gaussian random deviates were drawn from the ENDF/B-VII.1 and TENDL-2017 libraries where the covariances were available. Uncertainty from multiple nuclides was propagated simultaneously assuming inter-nuclide independence. When the individual statistical uncertainty is negligible compared to the data uncertainty, then standard probability theory may be applied. If this is not the case and both need to be considered, we use Imprecise Probabilities (IP) to perform further analysis. We show how uncertain experimental data may be compared to uncertain simulation in the context of IP, and show how an uncertainty-based sensitivity analysis can be performed with IP.


I. INTRODUCTION
Neutronics has an important role in fusion research. Essential reactor subsystems like tritium breeder and heating blankets are neutronic systems. Various diagnostics, magnets, and the outside environment must also be shielded from radiation, either from the high-energy neutrons produced directly in the core or from the secondary particles produced from radiation interactions with matter. Therefore, reliable neutronics calculations are essential for the success of experimental fusion reactors and eventual power plants. The validation of Monte Carlo particle transport codes for fusion neutronics relies on legacy experimental data, the Shielding Integral Benchmark Archive Database (SINBAD) benchmarks, many of which were performed in the 1980s. SINBAD is a collection of reactor shielding (46), fusion blanket neutronics (31), and accelerator shielding (23) benchmarks. 1 The database contains experimental details, along with experimental data and some computer code models (MCNP). In previous decades there have been efforts to assess the quality of SINBAD, 2 particularly aimed at nuclear data validation purposes. The benchmarks greatly vary in quality and comprehensiveness, and it has been suggested that the database would greatly benefit from experimental reinterpretation and reevaluation with modern methods, which would require a large empirical effort.
In this paper we present three benchmarks evaluated with OpenMC (Ref. 3), an open-source and community-developed Monte Carlo code for neutron and photon transport. We evaluate the Oktavian Iron, 4 Oktavian Nickel, 5 and the Frascati Neutron Generator (FNG) neutron streaming 6 shielding benchmarks, with a particular focus on nuclear data uncertainty propagation. The uncertainty propagation was performed with the Total Monte Carlo method 7 (TMC), with a proposed slight modification that allows for the nuclear data uncertainty to be combined with the statistical uncertainty of the transport simulation in a consistent way. If the statistical uncertainty is negligible, then we propose that the individual simulated means may be used as data points in a standard probabilistic analysis. If this is not the case, as is relevant in a variant of the TMC called Fast-TMC (Ref. 8) where larger statistical uncertainty is traded for computational time, then we propose that the analysis is performed with Imprecise Probabilities 9 (IP). Imprecise Probabilities is a generalization of probability theory that allows one to work with sets of distributions, which we use when the distributional output of transport simulation cannot be ignored. It is shown how comparisons between two uncertain simulations can be performed in the context of IP, along with how uncertain experimental data and simulations can be compared.

II. NUCLEAR DATA AND UNCERTAINTY PROPAGATION
Nuclear data are a fundamental input to any nuclear application. They effectively define all of the complex particle-matter interaction rules for the transport simulation. Nuclear data libraries usually provide at least three quantities for every nuclide-interaction pair: 1. interaction cross sections: characterize the probability of a particular reaction 2. exit energy/angle distributions: if the reaction has an exit particle, its next energy/angle state is sampled 3. covariances: the variance, autocorrelation, and cross-correlation information of the interaction cross sections.
The creation of nuclear data inputs for nuclear applications is a research field by its own right, known as nuclear data evaluation. Here it is the nuclear data evaluator's task to statistically mix experimental reaction data with reaction models to produce his/her best estimate of the nuclear data quantities plus uncertainty. Due to the difficulty and cost of conducting nuclear reaction experiments, experimental data are sparse or often not present for the vast majority of nuclides and reactions. Excluding the main fission-related nuclides, which have been extensively studied and have low uncertainties, the uncertainty can be quite severe for most nuclides, with the worst case being that the covariance information is completely missing in some evaluations. This means that the uncertainty problem is quite severe in fusion compared to fission, since not only are the fusionrelevant nuclides less well studied, but the variety of nuclides that must be considered is far greater (i.e., the reactor material can be quite exotic).
The nuclear data are usually stored in the Evaluated Nuclear Data Format (ENDF) format, which is a low-rank approximation by polynomial expansion, the coefficients of which are stored. This allows for the data to be accurately approximated and sorted by a finite number of polynomials. (It is well known that any continuous function can be arbitrarily well approximated by a polynomial.) The covariance, which is effectively a two-dimensional function for continuous energy, is also stored in this way. Although accurate, for Monte Carlo codes it is expensive to reconstruct the cross sections from the expansion each time an event is sampled. The nuclear data are therefore converted to ACE format, which is effectively a lookup table with the nuclear data values defined on an energy grid that can be quickly interpolated by a Monte Carlo code.
Historically, perturbation methods have been used for nuclear data uncertainty propagation. However, these methods require the model to be linear in the range of the uncertainty and are too strong an approximation for this severity of uncertainty. Currently only Monte Carlo methods have been successfully applied to fusion neutronics problems. 10 Monte Carlo is a very general method for propagating a probability distribution. Here random evaluations of model inputs are sampled (from a distribution that must be known or assumed, including correlations), and samples of the output distribution are obtained by repeatedly evaluating the model for each random input. There are some drawbacks, for example, Monte Carlo simulation is generally expensive (although can be easily parallelized) and only an approximation of the output is obtained. The output moments, such as the mean and covariance, may be computed in relatively few samples. However, if the full output distribution is sought, a a very large number of samples is required, particularly for the tails of the output (rare events). Monte Carlo a As it should be in any serious uncertainty or risk analysis.
UQ SINBAD WITH IMPRECISE PROBABILITY · GRAY et al. 803 methods also require the input distribution to be completely known. If there is any uncertainty about this input distribution, perhaps if any of the marginals or dependencies are unknown, then multiple Monte Carlo loops are required.

II.A. Nuclear Data Sampling
There are two techniques for nuclear data sampling: 1. using the covariance and assuming normality 2. a Bayesian method with a nuclear reaction model code.
For the first method, one can use a well-known method for generating a correlated Gaussian random vector, e.g., like Cholesky factorization of the covariance matrix. This can also be used for continuous functions like a cross section and is the same approached used in Gaussian process regression and Gaussian random fields. This method is made available for nuclear data by SANDY (Ref. 11), which generates random pointwise-ENDF evaluations from an ENDF file where the covariances are available. It should be noted that although the input distribution is normally distributed, the output distribution will not be. Normality is only generally preserved with linear transformations. Using this method, it is possible to generate a correlated random vector of nuclear data quantities that is from a specific distribution other than Gaussian, for example, uniform or lognormal. This can be done by firstly generating a correlated Gaussian random vector and inverting the samples through a Gaussian cumulative distribution function (cdf). This produces a correlated random vector with standard uniform marginals distribution (with range ½0; 1�). Samples from a specific distribution can then be generated by transforming these samples through the inverse cdf of the selected distribution. This is a form of inverse transform sampling but for stochastic processes. It should be noted that the covariance of the transformed samples and the original Gaussian vector may not be identical since the marginal distributions play a role in shaping the distribution dependency structure.
The second method is used by Koning and Rochman in the TMC and can be used for nuclear data evaluation as well as uncertainty propagation. The method relies on a nuclear reaction model code, such as TALYS (Ref. 12), that can produce a complete evaluated data set. They use Bayesian updating, a now popular uncertainty characterization method, to construct distributions for the input parameters of TALYS with the available experimental data. This input distribution can then be propagated through TALYS with Monte Carlo to construct a distribution of evaluated nuclear data quantities. The distribution from a Bayesian updating is usually non-Gaussian. Figure 1 shows a visualization of a TALYS Evaluated Nuclear Data Library (TENDL) cross section. They then suggest that either this distribution is condensed into a mean evaluation plus a covariance, which can then be propagated with another method, or that this distribution is simply propagated further with Monte Carlo. TENDL (Ref. 13) was created using this method and includes many covariance files missing from other libraries.
For Monte Carlo particle transport applications, this uncertainty propagation scheme can be considered as second-order Monte Carlo; that is, since each execution of the Monte Carlo particle transport simulation produces a distribution (a mean value and deviation corresponding to the Monte Carlo error), then repeated execution will produce a set of distributions. Figure 2 is an illustration of this. The overall dispersion of the distributions is due to the uncertain nuclear data input, with the variance of each individual distribution being due to the Monte Carlo error of that particular simulation. If the individual simulations are well converged, with a low statistical error, then

804
GRAY et al. · UQ SINBAD WITH IMPRECISE PROBABILITY a standard probabilistic analysis can be performed with the output using the individual means as samples. If, however, this is not the case, then the statistical errors must be considered for an accurate representation of the uncertainty. This is particularly relevant in a variant of the TMC called Fast-Total Monte Carlo, 8 where computational time is exchanged for larger statistical uncertainty. This is simply done by reducing the number of particles simulated in each individual transport simulation. If N particles are used in a simulation without uncertainty propagation, then it is proposed that N n particles are used in the individual TMC calculations, where n is the number of nuclear data samples. This of course can lead to a nonnegligible statistical uncertainty. Rochman et al. suggest that the following quantities are calculated: and where μ Total is the overall mean, and σ 2 Total is the overall variance, which is the sum of the variance of the means and the average statistical variance. There are some arguments that can be made against representing the uncertainty this way. First, any distinction between the Monte Carlo error and the uncertainty from the nuclear data is lost. Equivalently, this representation mixes the Frequentest and Bayesian interpretations of probability, which represent aleatory and epistemic uncertainty, respectively. These are generally conflicting and should not be aggregated. This is because the cross section is a single but unknown function and does not have a frequency. Its evaluation comes from a Bayesian method and its probability distribution is interpreted as a belief. The Monte Carlo error alternatively is rooted in the variability of the transport model itself and is a frequency. Second, this representation only gives the mean and the variance, which gives little information about the overall probability distribution. If normality is assumed for further analysis or to propagate this uncertainty, this representation would lead to an underestimation of the uncertainty. Last, note that the variance in Eq. (1) is only approximated. They state that Eq. (1) is accurate if σ MC < , 0:5 σ Total , with the 0.5 not being a strict limit and application dependent. We believe that this quantification can be improved, where the statistical and nuclear data uncertainty can be represented in a consistent way that is a more exact characterization of the output uncertainty.

II.B. Imprecise Probabilities
Imprecise Probabilities is a generalization of probability theory where computations can be performed with sets of distributions. This allows different types of uncertainty to be combined and expressed by the same mathematical structure. Randomness, stochasticity, or aleatory uncertainty is captured by the individual probability distributions in the set; while imprecision and vagueness, or epistemic uncertainty, is captured by the differences in the distributions within the set. Epistemic uncertainty may also be modeled by probability theory (Bayesianism), however, it is thought by most in the uncertainty quantification community that epistemic and aleatory uncertainty should not be modeled in the same probability distribution since this mixes the two interpretations of probability (frequency and belief). This has therefore led to the theory of IP, a theory of structures that models both variability and imprecision. There are many ways to mathematically represent a set of probability distributions, including 1. intervals 14 : sets of real numbers or a set of distributions where only the range is known 2. random sets 15 : a probability distribution of intervals or other sets 3. fuzzy numbers 16 : sets with nonbinary indicator functions 4. probability boxes (p-boxes) 17 : bounds on cdfs.
These structures were discovered independently and are often synonymous and can be translated from one to another. Imprecise Probability links all these theories into one. For a comprehensive overview of IP, and for a formal description of uncertainty and information in terms of these structures, Reference 18 is recommended.
For this work we will use p-boxes. A p-box represents a set of distributions with the following three constraints: (1) two bounding cdfs, (2) interval bounds on moments, and (3) the shape (e.g., normal). That is, a distribution FðxÞ is a member of a p-box if: where FðxÞ and FðxÞ are the p-boxes lower and upper cdf bounds, respectively; μ x and v x are the mean and variance; and F is a distribution family (i.e., normal, uniform). Some of these constraints may be missing, for example, if the p-box shape is unknown then any distribution that meets criteria 1 and 2 is considered. Some of these constraints may be informed from the other. For example, bounds on the moments can be produced from the cdf bounds, and cdf bounds may be constructed from moment information. These constructions generally use the Chebyshev, Markov, and Cantelli inequalities from probability theory. 17 Figure 3 shows various examples of p-boxes. If all three criteria are missing but the range is known, the p-box is simply an interval.
A p-box can easily be constructed from the output of a TMC calculation by taking the upper and lower envelope over all the simulated cdfs. The p-box essentially bounds all the properties of the TMC output. Arithmetic operations can be performed with p-boxes, with known or unknown correlations, which we use as tally arithmetic to compare between simulation and experiments or two uncertain simulations. How arithmetic is performed with p-boxes is out of the scope of this paper, however, it is made available by the OpenSource software ProbabilityBoundsAnalysis.jl (Ref. 19) developed by the authors.

II.C. Uncertainty-Based Sensitivity Analysis
An uncertainty propagation calculation is often complemented by a sensitivity analysis, where the contribution of the individual uncertain inputs to the output uncertainty is calculated. Such calculations are relevant if, for example, one wants to reduce the

806
GRAY et al. · UQ SINBAD WITH IMPRECISE PROBABILITY uncertainty by empirical effort but has limited resources and must select the most relevant inputs. In this section we describe how sensitivity analysis can be performed with p-boxes. This calculation follows a similar example in a recent paper, 20 where a sensitivity analysis was performed on a random function (stochastic process) with p-box outputs. The goal of a sensitivity analysis is to rank the uncertain inputs of a computational model with respect to their contribution to the output uncertainty. The ranking is given by sensitivity indices, which are values in ½0; 1�, with 0 indicating no contribution to the output uncertainty from a particular parameter and a 1 indicating the output uncertainty comes entirely from that particular input. A simple formula for calculating the sensitivity indices is 21 where U T is the total uncertainty, with all nuclides being varied, and U i is the output uncertainty without nuclide i. If there is no influence from nuclide i, then U i � U T and S i � 0. Equation (3) requires some metric or measure of the output uncertainty U i , which must be selected by the analyst. For precise distributions, the variance or the Shannon entropy of the output distribution is a typical choice. In the case of an output which is p-box or other imprecise structure, then the Hartley function 22 is sometimes used. For this work, like in Ref. 20, the area metric will be used to measure the uncertainty contained in a p-box since it is easy to compute and has a nice interpretation. The area metric, which was originally proposed for model validation under uncertainty, 23 is a measure of distance between two probability distributions and is a true metric between distributions in the mathematical sense (nonnegative, symmetric, triangle inequality, and identity of indiscernibles). The area metric d between two cdfs F and G is given by and is the horizontal average, or area, between F and G. Unlike some other stochastic distance metrics, the area metric has some favorable properties such as being unbounded (d 2 ½0; 1�), returns the physical units of F and G (if F and G are in meters, their distance d will also be in meters), and is a stochastic generalization of Euclidean distance (when F and G are scalars, d returns the Euclidean distance). The area metric can be used to measure the amount of epistemic uncertainty contained within a p-box by measuring the distance between its two bounding cdfs, which we use to calculate U i and U T in Eq. (3).

III. BENCHMARK RESULTS
The Oktavian Iron, Oktavian Nickel, and FNG neutron streaming SINBAD benchmarks were evaluated with OpenMC. The SINBAD database's MCNP Constructive Solid Geometry (CSG) model and material input files were converted to an OpenMC input file with csg2csg (Ref. 24). The ENDF/B-VII.1 (Ref. 25) and TENDL 2017 (Ref. 13) nuclear data libraries were used. ProbabilityBoundsAnalysis. jl was used for the postsimulation uncertainty analysis with p-boxes. Five-hundred random samples were drawn from each library using SANDY, b with the individual nuclides in the libraries being sampled separately. With many of the covariance files missing from ENDF/B-VII.1, only the nuclides with available covariances were sampled. Generally no inter-nuclide covariances are provided (dependence between the data of different nuclides), we therefore propagate uncertainties from multiple nuclides simultaneously assuming inter-nuclide independence. After nuclear data sampling, this can be done by randomly selecting nuclides from the sampled library and creating a list of nuclides to be used in the transport simulation, for example, { 12 C -259, 54 Fe -461, 56 Fe -166, 57 Fe -375, . . .}, where the second number beside the nuclide is the index of the nuclide in the random library. This was done 500 times (the total number of random evaluations) and is done without replacement so that each random nuclide only appears once. Figure 4 shows the results for Oktavian Iron. Figure 4a shows the simulated leakage from the iron sphere, with the mean of ENDF/B-VII.1 and TENDL-2017 in green and blue, respectively, and the red and gray envelopes showing the 95% confidence interval drawn from the p-boxes over all the simulations. The leakages were normalized by dividing the output by the surface area of the iron sphere. Figure 4b shows a comparison of the simulation to experimental results, shown with 95% confidence intervals. The absolute relative difference was calculated by subtracting the experimental distribution from the p-box created from the OpenMC simulation, dividing by the mean, and taking the absolute value. The 95% interval was then taken from the resulting p-box. It can be seen that many of the data points cross the agreement line, but with very large errors. b There are some samples from the Bayesian method available from TENDL, however, the number of nuclides is currently quite limited.
UQ SINBAD WITH IMPRECISE PROBABILITY · GRAY et al. 807 Figure 5 shows the results for Oktavian Nickel, showing the range of the uncertainty as a 95% confidence interval again drawn from the p-boxes created from the OpenMC output and representing both the statistical and nuclear data errors. An evaluation from JEFF-3.3 (Ref. 26) is also shown but without uncertainty in black. A comparison between the two uncertain leakages is shown in Fig. 5b. This was again calculated by subtracting one simulated p-box from the other and normalizing by the mean for each point in energy. The comparison of the means is shown in red, with the 95% interval. There is a significant overlap between the two uncertain simulations, with the zero line of agreement being contained within the errors. Table I shows the results of an uncertainty-based sensitivity analysis for the principle nuclides in Oktavian Nickel, along with the weight fraction of each nuclide in the nickel sphere. Intuitively, most of the indices follow a similar ranking to their weight fraction, except perhaps 64 Ni which ranks as the third most influential nuclide but with a low weight. These indices were calculated by summing the individual neutron currents and constructing a p-box by enveloping. The indices can then be calculated by measuring the drop in uncertainty [Eq. (4)] when performing the simulation without a particular nuclide being varied and applying Eq. (3). Figure 6 shows a comparison between the leakages when all nuclides are varied and when the three most important nuclides ( 58 Ni, 60 Ni, and 64 Ni) are individually left out of the TMC calculation. The sensitivity indices may also be calculated for each energy bin; and Fig. 7 shows this for 58 Ni and 60 Ni.

808
GRAY et al. · UQ SINBAD WITH IMPRECISE PROBABILITY Figure 8 shows the spacial flux of the FNG benchmark using ENDF/B-VII.1. Figure 8b shows the mean, with Figs. 8a and 8c showing 5% and 95% quantiles of the uncertainty, respectively. P-boxes were used, and the uncertainty shows the joint contribution from the nuclear data and the statistical uncertainties. It can be seen that the range of the uncertainty is largest at the back of the shielding, ranging from ,0 to ,10 À 6 . We believe this is because of the large number of uncertain interactions the particles have taken in reaching the back of the shielding. An explanation for this is that the uncertainty is cumulative after each uncertain reaction, meaning the more interactions the neutron has in its history, the larger its uncertainty. This may be relevant for large experiments like ITER, which have many meters of  UQ SINBAD WITH IMPRECISE PROBABILITY · GRAY et al. 809 shielding. Figure 9 shows the 58 Ni (n,p), 93 Nb (n,a), and 27 Al (n,2n) reaction rates as a function of depth for FNG, showing the distribution of the means. The covariance files for 93 Nb and 27 Al were missing from ENDF/B-VII.1 and can be seen to have had a large impact on the calculated output uncertainty.

IV. CONCLUSIONS
Uncertainties in fusion neutronics are more than in fission simulations. This is due to the large number of nuclides that must be considered in fusion problems, many of which have large or missing covariances. In this paper we performed nuclear data uncertainty propagation with the TMC in three SINBAD benchmarks using OpenMC. The ENDF/B and TENDL nuclear data libraries were used, with Gaussian random samples being drawn from the libraries' covariance files. Uncertainty from multiple nuclides was propagated simultaneously assuming independence between the uncertainty between different nuclides. P-boxes were used to model both the statistical uncertainty of individual particle transport simulations and the propagated nuclear data uncertainty. Comparisons between uncertain simulations and experiments were performed with p-boxes. The nuclear data uncertainty was found to be generally large in all the

810
GRAY et al. · UQ SINBAD WITH IMPRECISE PROBABILITY simulated benchmarks, but particularly at the back of the shielding for the FNG benchmark. We believe this is due to the large number of uncertain interactions the particles must undergo to reach the back of shields. With uncertainty being cumulative, the larger the number of interactions, the larger the uncertainty. Such conclusions are relevant for ITER and other fusion experiments with many meters of shielding. The authors plan to pursue such calculations in the future. An example of a reaction rate tally with uncertainty was also shown, with and without covariances. It was shown that inclusion of the covariances in a reaction rate tally has a large impact on the output uncertainty.
It should be noted that this approach of using the nuclear data covariance matrix to sample from a Gaussian distribution provides only an approximation to the output uncertainty. Ideally, the data-correct distribution should be used in the form of random files provided by the libraries. At the time of writing, this is not usually provided by the libraries, with the exception of TENDL, which provides random ENDF and ACE files for some select few nuclides, although even here the variety is quite limited. Such approximations may significantly effect the quantiles and dependencies of the output uncertainty, and the authors would encourage more libraries to provide random nuclear data files. This would also simplify the calculation in this paper, as it would remove the need for Gaussian sampling.