Importance of parametric uncertainty in predicting probability distributions for burst wait-times in fissile systems

Abstract A method of uncertainty quantification in the calculation of wait-time probability distributions in delayed supercritical systems is presented. The method is based on Monte Carlo uncertainty quantification and makes use of the computationally efficient gamma distribution method for prediction of the wait-time probability distribution. The range of accuracy of the gamma distribution method is examined and parameterised based on the rate and magnitude of the reactivity insertion, the strength of the intrinsic neutron source and the prompt neutron lifetime. The saddlepoint method for inverting the generating function and a Monte Carlo simulation are used as benchmarks against which the accuracy of the gamma distribution method is determined. Finally, uncertainty quantification is applied to models of the Y-12 accident and experiments of Authier et al. (2014) on the Caliban reactor.


Introduction
In a fissile system, the fluctuations in neutron population over time are driven by the branching of neutrons chains, an inherently random process. When the neutron population is large, the Law of Large Numbers determines that the outcome of the many branching processes taking place in the system, will tend towards the average (or ''expected") behaviour; the likelihood of significant deviations from the expected behaviour is small and can often be neglected. It is well known, however, that the behaviour of a fissile system when the neutron population is small, such as a reactor start-up in the presence of a weak neutron source, cannot be accurately modelled without considering the stochastic nature of the growth in these neutron chains. In these cases, significant deviations from the average behaviour are to be expected.

Relative strength of an intrinsic neutron source
A useful qualitative indication of the relative strength of a given neutron source was derived by Hansen (1960) who noted that a source should be considered weak if, 2Ss mC 2 where S is the neutron source strength in n/s, s is the prompt neutron lifetime, m is the average number of neutrons released per fission, and C 2 ¼ mðm À 1Þ= m 2 . This expression is approximately equivalent to the more simple inequality, 1.2. The wait-time The implications of a low neutron population are well documented and have been demonstrated experimentally, by Hansen (1960) and Authier et al. (2014), amongst others. In both examples, a fast burst reactor, GODIVA (Hansen) or Caliban (Authier et al.), was brought multiple times from a subcritical state to a delayed supercritical state, and the time taken to reach a pre-defined fission rate threshold was measured. The time taken to reach the threshold is known as the wait-time and it was shown to vary significantly between each realisation of the experiment, despite identical experimental conditions.
In a delayed supercritical system, the neutron population at any given moment consists of prompt neutrons emitted from fission and delayed neutrons emitted from those fission fragments which are delayed neutron precursors. In a delayed supercritical system, the population of delayed neutron precursors increases over time. This leads to an increase in the number of prompt neutron chains initiated and a corresponding increase in the fission rate and the https://doi.org/10.1016/j.anucene.2018.04.026 0306-4549/Ó 2018 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). prompt neutron population. The wait-time is the time taken for this build up in the neutron population to occur. It varies because the growth rate in the neutron population depends on the different events which can happen to each neutron emitted in the system (absorption, scattering, leakage, etc.) and these events are inherently random. The sequence of events can never be identical for two realisations of the same experiment and this can be observed on the macroscopic scale as variation in the wait-time.
The wait-time is an important concept that can have significant implications for the severity of accidental power excursions and for the safe start-up of nuclear reactors. During reactor startup, for example, control rods may be withdrawn to increase the reactivity of the system. Assuming the reactor behaves in a deterministic manner, the reactor power will begin to increase as soon as the system reactivity increases, and once the system is critical, an exponential increase in power should be observed. This will cause a rise in temperature, leading to negative reactivity feedback through material expansion and Doppler broadening, limiting the overall reactivity of the system and preventing an excessive increase in power. However, if the neutron source is too weak, or the withdrawal of the control rods too fast, there is likely to be a delay between reaching positive reactivity and any significant rise in reactor power. Meanwhile the reactivity of the system continues to rise in the absence of any negative temperature feedback. When the power output does finally begin to rise, the system reactivity may already be large enough to produce a dangerous power excursion. One potential objective in seeking to characterise the waittime is to know at what rate the control rod can be withdrawn so the probability of a dangerous power excursion remains below a specified safety probability, e.g. 10 À8 .

Methods for predicting the wait-time probability distribution
Methods for determining the wait-time probability distribution include Hansen's method (Hansen, 1960), the Fourier series method (see Abate and Whitt, 1992) and the saddlepoint method of Hurwitz et al. (1963). Hansen's method is approximative and based on neutron survival probabilities. Hansen's method considers that the wait-time consists of two parts: the time taken before a persistent neutron chain is initiated and the time for the neutron population due to the first persistent chain to build up to the waittime threshold. Persistent chains sponsored after the first are not considered to influence the wait-time significantly, because prompt neutron chains grow very rapidly; and unless the delay between the initiation of the first and second persistent chains were on the order of the generation time, the neutron population due to the second chain would be insignificant compared to that due to the first. Hansen notes that the initiation of a second persistent chain on this timescale is unlikely in a weak source scenario.
The Fourier series and saddlepoint methods rely on inversion of probability generating functions to obtain the probability distribution of the neutron population. This approach is more rigorous but also far more computationally expensive. These methods do not rely on the concept of the first persistent chain and are therefore able to account for the possibility of overlapping chains (initiated by different source neutrons) contributing to the neutron population at the moment the wait-time threshold is exceeded. This is important in delayed supercritical systems, where persistent fission chains consist of finite prompt chains linked by delayed neutron precursors. Since the delayed neutron precursors decay on a long timescale, compared to the generation time, overlapping fission chains become a significant possibility. The generating function methods can be applied to point models or space-dependent models. They can also be used with single or multiple energy groups.
A less expensive alternative to the generating function approach is to approximate the wait-time probability distribution using the gamma distribution method. This method was first proposed by Harris (1964) and relies on the fact that the neutron population in a multiplying system will tend towards a gamma distribution. There are cases when the neutron population does not conform to the gamma distribution so this method does not work for all scenarios, however it can be highly accurate and fast when applied to certain problems. This method has been applied by Williams (2016) to the Caliban experiments, with close agreement between the predicted wait-time probabilities and the experimental results of Authier et al. (2014) -see Section 1.5.

Types of uncertainty
For the purposes of the discussion that follows, it will be useful to distinguish between aleatoric uncertainty and parametric uncertainty. Aleatoric uncertainty will hereafter refer to uncertainty resulting from the random, stochastic nature of the buildup of neutron chains in fissile systems and parametric uncertainty will refer to that resulting from epistemic uncertainty in the input parameters. When calculating the wait-time probability distribution using any method, there will inevitably be some epistemic uncertainty in the input parameters, particularly when simulating accidental excursions where the exact chain of events may be unknown. The uncertainty in input parameters adds to the uncertainty already present due to aleatoric uncertainty.
The objective of this paper is to demonstrate a method for quantifying the impact of epistemic uncertainty in the input parameters. Uncertainty quantification (UQ) will be carried out using the Monte Carlo approach, which requires a fast method for determining the wait-time probability distribution, so that many calculations can be performed for a range of randomised sets of input parameters. The gamma distribution method will be used for this purpose, with verification of the predicted probability distribution using the saddlepoint method.
The ability to quantify the impact of epistemic uncertainty in the input parameters on the resulting wait-time probability distribution has important implications for criticality safety and safe reactor start-up.

Purpose of the current work
The method of uncertainty quantification presented in this paper makes use of the gamma distribution using the method presented in Williams (2016). In his paper, Williams shows that the wait-time probability distributions observed during the experiments of Authier et al. (2014) on the Caliban reactor, can be accurately predicted using the gamma distribution method. Starting from the forward form of the generating function equation, Williams derives equations for the mean neutron population, precursor group populations, detector counts and their corresponding covariances. These same equations will be used in this paper as the basis for the gamma distribution method.
The purpose of the current work is to incorporate the gamma distribution method demonstrated by Williams into a Monte Carlo algorithm, for the purpose of uncertainty quantification. The gamma distribution method for determining the wait-time probability distribution is particularly amenable to Monte Carlo uncertainty quantification due to its excellent computational efficiency compared to alternative, more rigorous methods.
The gamma distribution method has been shown to produce accurate results for certain scenarios, however it will be shown in Section 2.3 of this paper that the neutron population does not always conform to a gamma distribution. The present work will examine the parameters influencing the accuracy of the gamma distribution method in order to establish the range of transients to which this method of uncertainty quantification can be applied.

Gamma distribution method
The neutron population, n, is modelled by the following probability density function (the Gamma Distribution) in which n is approximated as a continuous random variable, with the time, t, as a parameter.
Pðn; tÞ ¼ gðtÞ nðtÞCðgðtÞÞ The function CðzÞ is the complete gamma function, defined as follows CðzÞ ¼ Two incomplete gamma functions also exist: the lower and upper incomplete gamma functions, where the upper or lower limit of integration, respectively, is replaced with a finite limit. The lower incomplete gamma function, for example, is defined as follows The time-dependent parameter gðtÞ is the ratio of the squares of the mean and standard deviation, The methodology used in this paper for calculating the wait-time probability distribution from the gamma distribution is based on that described in Williams (2016) with some minor modifications. The probability distribution, P w ðtÞ, of the wait-time may be determined by differentiating the cumulative probability density function Q ðn Ã ; tÞ, where n Ã is the instantaneous neutron population corresponding to the wait-time threshold fission rate, with respect to time.
The values of Q ðn Ã ; tÞ are obtained directly from the ratio of incomplete and complete gamma functions, as follows, Q ðn Ã ; tÞ ¼ cðgðtÞ; gðtÞn Ã = nðtÞÞ where cðx; yÞ and CðxÞ are the lower incomplete and complete gamma functions discussed above. The wait-time may also be calculated with respect to the cumulative detector count (Z Ã ), as in Williams (2016). In this case, the wait-time probability is obtained by differentiating the cumulative probability density function, Q ðZ Ã ; tÞ, for the detector count, Z Ã , corresponding to the wait-time threshold fission rate.
The parameter gðtÞ is replaced by g z ðtÞ, The value of n Ã or Z Ã is determined from the pre-defined wait-time fission rate threshold by calculating the expected fission rate, FðtÞ, using the following equation from Williams (2016), where kðtÞ is the time-dependent neutron multiplication factor, m is the mean number of neutrons released per fission, s is the prompt neutron lifetime and NðtÞ is the time-dependent neutron density. The value taken for n Ã is the value of NðtÞ at the moment FðtÞ crosses the wait-time threshold. Likewise, if the wait-time is calculated based on cumulative detector count, then the value of Z Ã is the value of ZðtÞ when FðtÞ crosses the fission rate threshold. The mean neutron population, cumulative detector count and their standard deviations were calculated by solving the system of ordinary differential equations found in Williams (2016) using the ODE solver of Shampine and Gordon (1975).
The model is zero-dimensional, one-speed and without thermal feedback. It is derived by differentiating the one-speed, point model forward equation (Bell, 1963), for the probability generating function.

Method overview
The wait-time probability distributions calculated by the gamma distribution method were compared against calculations using the saddlepoint method to confirm their accuracy. A detailed description of the saddlepoint method can be found in Williams and Eaton (2017) and only a brief summary will be presented here. The cumulative probability density function Q ðn Ã ; tjsÞ is obtained from Eq. (13). It represents the probability that the neutron population has exceeded n Ã by time t if a source was present since time s. It is equivalent to the Q ðn Ã ; tÞ of Eq. (8) if s is equal to the start time.
where G S ðz 0 ; tjsÞ is the backward form of generating function for the probability distribution of the neutron population at time t, given a source present since time s, and r 0 is given by, where G 0 S and G 00 S are the first and second derivatives of the generating function G S with respect to z. Finally, z 0 is the value of z which satisfies the equation, where n Ã is the neutron population corresponding to the wait-time fission rate threshold. The backward generating functions for a point model were obtained from Williams and Eaton (2017) and solved using the ODE solver of Shampine and Gordon (1975). The saddlepoint method is actually an approximation, albeit a very accurate one, because it relies on the method of steepest descents to approximate a line integral before arriving at the expression in Eq. (13). The error in the method does, however, become significant for values of Q ðn Ã ; tjsÞ close to unity, and Hurwitz notes that for values of n Ã ) n; Q ðn Ã ; tjsÞ tends to e ffiffiffiffi 2p p ¼ 1:08444 . . .. This results in a wait-time probability distribution whose integral is $ 1:08 and it is therefore necessary to renormalise the distribution to have an integral of 1:0. Notwithstanding these limitations, which result in a small, quantifiable error in the saddlepoint method's predictions, the accuracy of the saddlepoint method has been evaluated by Williams and Eaton (2017) who compared it to the exact Fourier series method, details of which can be found in Abate and Whitt (1992). Willams and Eaton found that the values of Q ðn Ã ; tjsÞ calculated using the saddlepoint method were within 0.5% of the exact values calculated using the Fourier series method, in all the cases they tested. It is therefore assumed throughout this paper that the results of the saddlepoint method are close to the true waittime probability distribution, and any significant discrepancies between the probability distributions predicted by the two methods are attributable to error in the gamma distribution method.
The reason that the saddlepoint method is valid in all cases, whereas the gamma distribution is not, is that the saddlepoint method inverts the generating function to obtain the neutron population probability distribution directly from known quantities (such as neutron cross-sections). The gamma distribution method, on the other hand, assumes a probability distribution for the neutron population based on its mean and standard deviation. While this method is accurate in many cases, it is known that the probability distribution for the neutron population can deviate significantly from the assumed distribution.
As an additional check to confirm the accuracy of the saddlepoint method, one set of data, for a situation where the gamma distribution method fails to give accurate results, was also compared against a Monte Carlo simulation. Details of this comparison can be found in Section 2.3.1.

Finding the value of z 0
The value of z 0 required by Eq. (14) is determined by finding the roots of the Eq. (15). Because the generating functions G S and G 0 S are functions of z 0 whose values can only be determined by solving a system of ODEs, the roots cannot be found analytically and an iterative solution is necessary. Determining the value of z 0 for all time steps can therefore be an expensive task. Hurwitz et al. (1963) provided the following approximation for z 0 which can be used as a starting guess for iterations, The method of bisection is a reliable technique for finding the roots of Eq. (15), and thereby determining z 0 . Nonetheless, this method requires a relatively large number of iterations to achieve the required accuracy in z 0 , leading to long computing times, particularly when modelling fast systems where the system of ODEs for the generating functions tends to be particularly stiff. A more efficient solution is proposed by Williams and Eaton (2017) who suggest that the roots of the equation could be found using Newton-Raphson iterations. This was found to reduce the number of required iterations dramatically; a typical simulation requiring approximately 20-30 iterations with bisection was found to need only 4-8 iterations using Newton-Raphson.
The Newton-Raphson method allows the roots of an equation f ðxÞ ¼ 0 to be determined through a series of iterations. In each iteration, an estimate for the value of the root x n is calculated based on the value of the previous estimate x nÀ1 . In most cases, the estimate converges rapidly on the true value of the root. Each new estimate of the root x n is determined from the following relation, where f 0 ðxÞ is the first derivative of f ðxÞ with respect to x.
In order to implement the Newton-Raphson method to find z 0 it is necessary to express it in the form f ðz 0 Þ ¼ 0 and this can be done in many ways. The form given in Eq. (18) was chosen, from which the derivative in Eq. (19) immediately follows.
The second derivative of the generating function G 00 S is already required in Eq. (14) to determine r 0 so no extension is required to the system of ODEs.

Range of applicability of the gamma distribution method
The gamma distribution method works well for simulating fissile systems in which the neutron population conforms to a gamma distribution. If this is not the case, however, then the predicted wait-time probability distribution will not be correct. It is therefore necessary to determine the range of applicability of the gamma distribution method if it is to be used with confidence. This was achieved by verifying the results of the gamma distribution method against the saddlepoint method for a range of different scenarios, varying the magnitude and rate of the reactivity insertion, as well as the prompt neutron lifetime and the neutron source strength. All calculations in this section include 6 groups of delayed neutrons. Fig. 1 shows mean wait-times and standard deviations predicted by the gamma and saddlepoint methods for a 0.7$ step insertion over a range of different intrinsic neutron source strengths. The means and standard deviations predicted by the two methods converge as the neutron source strength is increased, showing that the accuracy of the gamma distribution method improves with increasing neutron source strength. The relative offset in the mean predicted wait-time converges to a value around 1.2%, indicating that the mean value predicted by the gamma distribution method is 1.2% larger than that predicted by the saddlepoint method. This is likely due to the previously mentioned weakness in the saddlepoint method, where the cumulative probability distribution for n ( n Ã tends to a value greater than unity. This skews the distribution towards lower values of n and leads to a slight underestimate of the overall mean of the wait-time probability distribution. Fig. 2 shows the wait-time probability distributions calculated using both the saddlepoint and gamma distribution methods for a range of reactivity step insertions in the presence of a neutron source emitting 30 n/s. The agreement between the two methods is very close for the 0.1$ reactivity step insertion, however the gamma distribution method becomes progressively less accurate as the magnitude of the reactivity insertion increases. For the larger reactivity insertions, the gamma distribution method overestimates the likelihood of longer wait-times, resulting in higher means and standard deviations than those predicted by the saddlepoint method (see Table 1).
The accuracy of the gamma distribution method improves when the neutron source strength is increased to 90 n/s (see Fig. 3). Significant inaccuracies are still observed, however, in the wait-times predicted for larger reactivity insertions. Fig. 4 shows the minimum source strengths required for the gamma distribution method to achieve acceptable accuracy over a range of reactivity step and ramp insertions. The system modelled in this figure has a prompt neutron lifetime of 45 ls. Acceptable accuracy was defined as a mean wait-time predicted by the gamma distribution method that was within 2% of the mean wait-time predicted by the saddlepoint method. No criterion was imposed on the standard deviation because it was found that the standard deviations converged with the predicted means (see Fig. 1). Fig. 4 shows that the minimum required source strength for the gamma distribution method to be close to the saddlepoint method increases with the magnitude of the reactivity insertion for step insertions and 1 s ramp insertions. The required minimum source strength also increases for a 10 s ramp insertion but the effect is smaller. For a 30 s ramp insertion, there is no correlation between the magnitude of the reactivity insertion and the minimum source strength required for accurate wait-time prediction by the gamma distribution method.
The minimum required source strength was also found to be a function of the prompt neutron lifetime. Fig. 5 shows how the required minimum source strength varied with the magnitude of the reactivity insertion and the prompt neutron lifetime. The source strength corresponding to 2% accuracy initially increased with decreasing prompt neutron lifetime, however for prompt neutron lifetimes shorter than $10 ls, no significant further increase in the minimum required source strength was observed. Fig. 5 includes markers to show the source strength and reactivity insertion size corresponding to two sets of burst wait-time experiments: the Caliban experiments of Authier et al. (2014) and the GODIVA experiment of Hansen (1960). From the figure it would be expected that the gamma distribution method would prove accurate for the Caliban experiments but not for Hansen's experiment on GODIVA. Indeed this is the case: it has already been shown by Williams (2016) that the Caliban experiments can be modelled accurately using the gamma distribution method. Hansen's experiment, on the other hand, requires a more rigorous      method due to the combination of a relatively large and rapid reactivity insertion and a weak neutron source.

Verification with Monte Carlo simulation
A 0:7$ step insertion case was simulated using a Monte Carlo code, the results of which are shown in Fig. 6 and summarised in Table 2. The Monte Carlo simulation was performed using the code described in Cooling et al. (2016). The code performs a large number of realisations (5000 in this case), each with the same input parameters, and within each realisation, randomly determines the fate of each neutron produced at each time step according to the physical probabilities of different events.
The results match well with those of the saddlepoint method but not the gamma distribution method, which supports the hypothesis that the saddlepoint method works well even for scenarios where the gamma distribution method is not valid. The probability distribution generated by the Monte Carlo code has a jagged appearance, whereas the distributions generated by the gamma distribution method and saddlepoint methods are smooth. This is because the Monte Carlo code works by calculating the wait-times for a large number of system histories. In order to construct a continuous probability distribution, it is necessary to group similar wait-times together into 'bins'. The probability corresponding to each wait-time bin is then inferred from the proportion of system histories whose wait-time falls within each bin.
It was also attempted to use the code of Cooling et al. to simulate the 0.1$ and 0.5$ step insertion cases, however it was found that the value of N thresh (see Cooling et al., 2016 for details) required to converge the resultant probability distribution was too large to allow the calculations to be run in a reasonable time. Essentially, the low reactivity in these cases allowed a large number of delayed neutron precursors to build up before the power increased to a level where its behaviour could be considered deterministic. Simulating the decay of such a large number of delayed neutron precursors proved too expensive in terms of computing time.

Types of uncertainty
Monte Carlo uncertainty quantification was carried out in two case studies to examine the impact of parametric uncertainty on the predicted wait-times. Four different wait-time probability distributions are presented: Deterministic This probability distribution is a delta function, representing the expected wait-time when all uncertainty (both aleatoric and parametric) is neglected. The deterministic waittime is the time at which the mean (or expected) fission rate crosses the fission rate threshold. Aleatoric This is the distribution of wait-time probabilities resulting from the stochastic nature of the processes leading to the growth and branching of neutron chains. No uncertainty in input parameters is taken into account. Parametric This is the distribution of deterministic wait-time probabilities when uncertainty is applied to one or more of the input parameters. No aleatoric uncertainty is taken into account. Aleatoric-Parametric This is the wait-time probability distribution when aleatoric and parametric uncertainty are both taken into account. It is expected that the distribution incorporating the combined aleatoric and parametric uncertainties will be broader and less peaked than either the aleatoric or parametric distributions.

The Monte Carlo code
The Monte Carlo uncertainty quantification was implemented using the Message Passing Interface (MPI) to run multiple instances of the gamma wait-time code in parallel. In each instance, the Monte Carlo code randomised the variable to be made uncertain, before passing this to the gamma wait-time code, which returned either the deterministic wait-time or the calculated waittime probability density function.
The parametric only probability distribution was then constructed by establishing wait-time bins of equal width across the range of expected wait-times. The MPI code was used to calculate 100,000 values of the deterministic wait-time by running 100,000 realisations of the wait-time code. These values were then sorted into the wait-time bins, and the probability distribution inferred from the proportion of values returned in each bin.
For the aleatoric-parametric probability distribution, the MPI code was used to calculate 10,000 aleatoric wait-time probability density functions by running 10,000 realisations of the gamma wait-time code. The combined probability density function was where P w ðtÞ is the probability of exceeding the wait-time threshold at time t; N rn is the number of Monte Carlo realisations and P w;i ðtÞ is the probability of exceeding the wait-time threshold at time t according to distribution i. Note that P w ðtÞ is the wait-time probability across all realisations representing the full range of input parameters, whereas P w;i ðtÞ is the wait-time probability for a single realisation with a single set of input parameters. A larger number of realisations is required to obtain the parametric-only distribution than the parametric-aleatoric distribution. This is due to the need to sort the calculated deterministic wait-times into bins so that the continuous probability distribution can be approximated. The total number of realisations must be sufficiently large that the number of wait-time values in each bin is large compared to the standard deviation. The total number of realisations required depends on the number of bins used. On the other hand, each realisation of the parametric-aleatoric simulation produces a continuous probability distribution, therefore no binning is required and the number of realisations can be smaller.
Execution times depended strongly on the scenario simulated. The simulations presented in the next Section take approximately 60 min each, running on a 20-core desktop computer. This corresponds to 0.4 s per realisation for the aleatoric-parametric calculation and 0.04 s for the parametric only calculation. Individual realisations of the parametric only calculation are faster than the aleatoric-parametric calculation because only a single wait-time value is returned, whereas in the aleatoric-parametric calculation, gamma functions must be evaluated in order to construct the wait-time probability distribution for each realisation. Since the Monte Carlo method is easily parallelised, execution time decreases in inverse proportion to the number of cores available.
The prompt neutron lifetime has a particularly significant impact on the required computation time, as shorter prompt neutron lifetimes increase the stiffness of the ODEs that must be computed in order to calculate the gamma distributions. Williams (2016) observed that the wait-time probability distribution is sometimes insensitive to the prompt neutron lifetime; in cases where the relative standard deviation in the neutron population (or detector count) reaches a constant value before the probability of reaching the threshold value has risen above a negligible value. In these cases, the execution time can be reduced by selecting a longer prompt neutron lifetime, thereby reducing the stiffness of the system of ODEs. This technique was used to accelerate the execution of the simulations presented in Section 4.2.

Mean and standard deviation of the distributions
The means and standard deviations were calculated for each of the wait-time probability density functions obtained from the Monte Carlo uncertainty quantification. The method used to calculate these values depended on the type of distribution.
For the aleatoric and aleatoric-parametric distributions, which are continuous distributions, the mean wait-times and standard deviations were calculated by numerically integrating the following expressions using Simpson's rule, where l is the mean wait-time, r is the standard deviation in the wait-time and PðtÞ is the probability of exceeding the wait-time threshold at time t.
For the parametric only distributions, which are discreet, the mean wait-times and standard deviations were calculated from the following expressions, where N rn is the number of Monte Carlo realisations and t i is the wait-time calculated in the ith realisation.

Case study I: The Y-12 accident
The gamma distribution method was applied to predict the timing of the first power excursion in the criticality accident that occurred in 1958 at the Y-12 National Security Complex in Oak Ridge, Tennessee. As with most criticality accidents, there is some uncertainty around the exact conditions (flow rates, concentrations, etc.) prior to the accident. The impact of uncertainty in the flow rate of the uranium-rich solution will be examined by means of Monte Carlo uncertainty quantification.

Description of the accident
A 208 litre drum was being used for the collection of water drained from pipework during a leak test. According to Patton et al. (1958), the water used for leak testing was to be transferred to the drum via a system of pipes, however a leaking valve meant that these pipes had become filled with a solution of 90% highly enriched uranium (HEU) in the form of uranyl nitrate. The solution of HEU in the pipework was deposited in the drum and was followed by water from the leak testing. Initially, the HEU deposited in the drum remained subcritical, however as water from the leak test flowed into the drum, the degree of moderation was increased, causing the solution to become critical. Water from the leak test continued to flow, making the solution more dilute, and eventually the system became subcritical once again. The total yield was estimated at 1:3 Â 10 18 fissions -a particularly large yield relative to other criticality accidents.
According to McLaughlin et al. (2000), it is possible that the drum had already reached a prompt critical configuration by the time of the first power peak. This feature makes an analysis using the method presented in this paper very relevant.

Model parameters
The aim of the model is to obtain the probability distribution of the wait-time between the system reaching criticality and the fission rate exceeding some pre-defined threshold. It is estimated in Zamacinski et al. (2014) that the maximum fission rate during the first power peak was between 10 16 and 10 17 fissions per second. Preliminary calculations indicate that this fission rate is likely to have occured after prompt criticality was reached. The gamma distribution method has not been validated in the prompt critical region so the threshold value will be set several orders of magnitude lower than the peak power level so that the delayed supercritical phase of the excursion may be examined.
The equations for the mean and covariances, derived in Williams (2016) from the forward form of the probability generating function, will be used to determine the time-dependent mean and standard deviation of the neutron population which characterise the gamma function. These equations do not include thermal-hydraulic feedback so they are only valid for the period before the fission power output rises to a level at which feedback becomes significant. This is another reason to set the wait-time fission rate threshold at a relatively low level so that the power output is low enough that thermal-hydraulic feedback may be neglected.
Delayed neutron precursors are represented in 6 groups using the same parameters as those used by Zamacinski et al. (2014) for their point kinetic model of Y-12. The reactivity profile as a function of time is also taken from Zamacinski et al. (2014), as is the mass and concentration of 235 U in the fuel solution. Zamacinski et al. determined that the prompt neutron lifetime varied approximately linearly from 4:1 Â 10 À5 s at the start of the transient, to 5:5 Â 10 À5 s at the end. Eq. (54) from Zamacinski et al. (2014) was used to model the prompt neutron lifetime in this analysis.
The probability distribution for the number of neutrons released per fission was taken from Table 1 of Zucker and Holden (1986).
An analysis of the various phenomena contributing to the intrinsic neutron source in uranium-fueled reactors with nonirradiated fuel can be found in Harris (1960). The report identifies three main sources of neutrons, 1. spontaneous fission of 235 U and 238 U, 2. (a,n) interactions where a particles emitted in radioactive decay of 238 U induce neutron emission in light nuclei, 3. cosmic ray induced emission of neutrons from within the fuel solution.
Experimental measurements of the intrinsic neutron source of solutions of uranyl nitrate can be found in Hankins (1966) and Seale and Anderson (1991) gives a total intrinsic source of 315 n/s. It should be noted that there is considerable uncertainty in this figure as it relates to the Y-12 accident. Firstly, it is not clear whether or not the uranyl nitrate solution present at the Y-12 accident had been previously irradiated. If it had been, there would be a higher concentration of radioactive isotopes, whose decay would contribute to an increased intrinsic neutron source through (a,n) reactions. Secondly, as noted above, the actual uranium concentration was slightly lower than the value in Hankins' experiment, which would be expected to reduce the intrinsic source due to reduced alpha emission from 238 U as well as reduced spontaneous fission.
Due to the uncertainty in the intrinsic neutron source, results will be presented for two intrinsic source strengths: 30 n/s and 315 n/s, so that the degree of sensitivity to the intrinsic source strength can be examined.

Parametric uncertainty
In order to examine the impact of uncertainty in the fuel solution flow rate on the resulting wait-time probability distribution, uncertainty was applied to the reactivity insertion rate. The reactivity insertion rate was modelled using Eq. (25), based on Eq. (53) from Zamacinski et al. (2014), who calculated the reactivity as a function of time from the available hydraulic and geometric data using MCNP. The randomly varying parameter a was included to make the reactivity insertion rate uncertain. Its value was randomly varied between 0:9 and 1:1 according to a uniform probability distribution. The time axis of the equation has also been modified so that the system reaches criticality at t ¼ 0.
R ex ðtÞ ¼ 1:8636 Â 10 À2 at À 5:5338 Â 10 À5 ðatÞ 2 þ 6:8570 Â 10 À8 ðatÞ 3 ð25Þ Fig. 7 shows the system reactivity as a function of time, including the non-randomised reactivity profile corresponding to a ¼ 1, as well as minimum and maximum values corresponding to a ¼ 0:9 and a ¼ 1:1. Considering, for example, the effect of a value of a > 1, two effects are notable. Not only does it result in the system reactivity rising faster, making shorter wait-times more likely; it also means the system is more deeply subcritical during the time before criticality is achieved. The latter effect offsets the first to some degree, because at the moment of reaching criticality, the population of delayed neutron precursors will be smaller in a system with large a than in a system with small a.

Deterministic wait-time and aleatoric uncertainty
The deterministic wait-time can be obtained by solving the system of ODEs referred to in Section 2.1 for the moment at which the expected fission rate FðtÞ crosses the wait-time threshold of 2 Â 10 9 fissions per second. The value obtained is 45.7 s for an intrinsic source strength of 30 n/s, and 43.0 s for an intrinsic source strength of 315 n/s. Unsurprisingly, the stronger intrinsic source results in a shorter wait-time. Fig. 8 shows the aleatoric uncertainty associated with these wait-times; that is uncertainty due to the stochastic nature of fission chain growth, assuming certain input parameters. Wait-time probability distributions, calculated using the gamma distribution method and saddlepoint method, are shown for comparison. The predicted distributions are similar: for the 30 n/s case, the gamma distribution method gives a mean (and standard deviation) of 47.0 s (2.1 s) and the saddlepoint method predicts 47.1 s (2.6 s). For the 315 n/s case, the gamma distribution method gives a mean (and standard deviation) of 43.1 s (0.8 s) and the saddlepoint method predicts 43.0 (0.8 s). Fig. 9 shows the four probability distributions (deterministic, aleatoric, parametric and aleatoric-parametric) for the wait-time to reach 2 Â 10 9 fissions per second in the Y-12 accident. Fig. 9a shows the probability distributions for the case where the intrinsic neutron source strength is 30 n/s and Fig. 9b shows the distributions for the 315 n/s case.

Parametric uncertainty quantification
Despite the uniform probability distribution adopted for the randomisation of the reactivity insertion rate, the parametric wait-time probability distributions are slightly biased in favour of shorter wait-times. The is because the impact of varying alpha on the wait-time predicted by the deterministic model is not linear.
The parametric uncertainty makes the overall combined (aleatoric-parametric) probability distributions broader, compared to the aleatoric probability distributions, however, in both cases, there is no significant change in the mean. For 30 n/s case, the mean (and standard deviation) of the aleatoric-parametric waittime probability distribution is 47.1 s (2.9 s) compared to 47.0 s (2.1 s) for the aleatoric distribution. For the 315 n/s, the aleatoric-parametric distribution has a mean (and standard deviation) of 43.3 s (2.0 s) compared to 43.1 s (0.8 s) for the aleatoric distribution.
The parametric uncertainty has a greater impact on the standard deviation of the probability distribution in the 315 n/s case. This is because there is less aleatoric uncertainty, so the parametric uncertainty is more significant by comparison. This effect is clearly visible in Fig. 9b, where the combined aleatoric-parametric uncertainty is dominated by parametric uncertainty, whereas in Fig. 9a, the combined aleatoric-parametric distribution is influenced in approximately equal parts by the aleatoric and parametric uncertainties.
These results show that a small degree of uncertainty in the reactivity insertion rate results in a significant increase in the overall uncertainty in the wait-time.

Summary of results
The mean wait-times and standard deviations discussed in this section are summarised in Table 3.

Case study II: the Caliban experiments
The Caliban reactor is a fast burst reactor constructed of solid metal disks of uranium/molybdenum with a 235 U enrichment of 93.5%. Authier et al. (2014) describe a series of experiments in which varying amounts of reactivity were inserted into the reactor so that the wait-time could be measured. It has already been demonstrated in Williams (2016) that the gamma distribution method to determine the wait-time probability distribution is accurate for these experiments and it is clear from Fig. 4 that these experiments fall within the range of applicability of the gamma distribution method.

Description of the experiments
Full details of the experimental set-up can be found in Authier et al. (2014). Reactivity control in the Caliban reactor is achieved in part by means of calibration rods and a safety block. The calibration rods are set so that the desired degree of supercriticality is achieved when the safety block is lifted to its maximum position. The reactor is also equipped with a burst rod but this was not used for the delayed supercritical experiments modelled in this paper.
Each delayed supercritical excursion was initiated by raising the safety block progressively closer to the reactor until it was in its maximum position. According to Fig. 3 of Authier et al. (2014) the safety block is moved in 17 small movements over a period of 62 s until the block is in the desired delayed supercritical position. In the model described here it was chosen to lump together some of these movements and the 17 movements of the safety block were approximated by 6 small ramp reactivity insertions (see Eq. (26)).
Authier et al. count the wait-time from the moment the safety block reaches its maximum position to the time at which the neutron detector indicates a fission rate of 2 Â 10 9 fissions per second. According to Table I of Authier et al. (2014) the reactor reaches a supercritical configuration approximately 4 s before the safety block arrives at its maximum position.

Model parameters
The experiment selected for uncertainty quantification was the 0.272$ reactivity insertion. The wait-time threshold was set to 2 Â 10 9 fissions per seconds, the same value used by Authier et al. (2014). Delayed neutron precursors are represented in 6 groups using the parameters for fast fission of 235 U from Keepin (1965).
R ex ðtÞ ¼ À16:7$ for t 6 À62:1 s À16:7$ þ tþ62:1s 16:7 s Â 1:7$ for À 62:1 s < t 6 À45:4 s À15:0$ þ tþ45:4 s 16:7 s Â 3:6$ for À 45:4 s < t 6 À28:8 s À11:5$ þ tþ28:8s 16:7 s Â 6:6$ for À 28:8 s < t 6 À12:1 s À4:93$ þ tþ12:1s 8:1 s Â 4:9$ for À 12:1 s < t 6 À4:0 s 0:066$ þ tþ4:0 s 4:0 s Â 0:206$ for À 4:0 s < t 6 0:0 s 0:272$ for t > 0:0 s The reactivity profile is shown in Eq. (26). It was constructed using data from Table 1 of Authier et al. (2014). The final reactivity increase in the table was omitted, as this corresponds to the insertion of the burst rod (prompt critical experiments only), and the reactivity values were adjusted by þ0:077$ so that the final reactivity corresponded to 0:272$. The beta effective used in the model was 0:00659. For the intrinsic neutron source of the Caliban reactor, Authier et al. (2014) used a value of 200 n/s in their model, so the same value will be used for this analysis. No external source was present during the experiments. The prompt neutron lifetime in the Caliban reactor has been evaluated as 1:2 Â 10 À8 s À1 (see Casoli et al., 2009). Attempting to model the reactor with such a short prompt neutron lifetime results in a stiff set of ODEs, which increases the computational expense of the Monte Carlo analysis and limits the maximum number of realisations achievable. It is noted however, in Williams (2016), that the relative standard deviation in the cumulative detector count reaches a constant value, becoming insensitive to prompt neutron lifetime within 100 s of reaching criticality. This implies that the wait-time probability distribution is also insensitive to prompt neutron lifetime after this point. Since the expected wait-time is much longer than 100 s, it should be possible to increase the prompt neutron lifetime to improve program run-time, without impact on the predicted wait-time probability distribution. This hypothesis was tested by comparing the wait-time probability distribution predicted with a prompt neutron lifetime of 12 ns to a simulation with a prompt neutron lifetime of 65 ls.
The results, shown in Fig. 10, confirm the insensitivity of the wait-time probability distribution predicted for the Caliban experiment to the prompt neutron lifetime.

Parametric uncertainty
Two areas of parametric uncertainty were examined for the Caliban experiments: time taken to move the safety block and the yields of the delayed neutron precursor (DNP) groups.
Eq. (26) shows the reactivity profile representing the movement of the safety block in the Caliban experiments. The final ramp insertion lasts 4:0 s and takes the reactor from a slightly supercritical state (0:066$) to the target reactivity of 0:272$, for a reactivity insertion rate of 5:15 Â 10 À2 $/s. The duration of this step was varied between 0 and 20 s according to a uniform probability distribution so that the reactivity insertion rate varied from a step to a ramp of 1:03 Â 10 À2 $/s. The yields of the delayed neutron precursor groups are subject to some experimental uncertainty in their measurement. Table I of Keepin et al. (1957) states the standard deviation associated with each yield. From this value, an upper and lower bound (see Table 4) was calculated by assuming a uniform probability distribution. The values of each yield were then varied randomly within these ranges.
The total delayed neutron fraction, b, was calculated for each realisation as the sum of all precursor group yields. Its value therefore also varied, with a minimum of 0.00611 and a maximum of Table 3 Summary of calculated wait-times to reach a fission rate of 2 Â 10 9 s À1 for the simulation of the Y-12 accident. Mean values shown with standard deviation in brackets.

Deterministic
Aleatoric Parametric Aleatoric-Parametric 0.00707. The reactivity profile was fixed before randomisation of the delayed neutron precursor yields so that the time-dependent k eff was the same for each realisation. The reactivity of the reactor measured in dollars therefore varied with the changing value of b.
For the safety block in its final position, the reactivity in dollars varied between a minimum of 0:254$ and a maximum of 0:293$.

Deterministic wait-time and aleatoric uncertainty
The deterministic wait-time for the system to reach 2 Â 10 9 fissions per seconds is 207 s. This value ignores both aleatoric, and any parametric uncertainty. When aleatoric uncertainty was modelled using the gamma distribution method, a probability distribution was obtained, with mean and standard deviation 211 s (12.8 s). The distribution obtained is shown in Fig. 11.

Safety block manipulation time
Adding parametric uncertainty to the timing of the safety block movement in the deterministic wait-time model produces the parametric probability distribution shown in Fig. 12. The parametric probability distribution has mean and standard deviation 202 s (1.8 s).
The results confirm that the uncertainty applied to the timing of the final movement of the safety block has a significant impact on the resulting wait-time. This is because the system reaches criticality before the safety block reaches its final position. During this time the population of delayed neutron precursors in the system is increasing exponentially. The time taken to carry out this movement therefore will determine the delayed neutron precursor population at the start of the experiment with a resulting impact on the wait-time observed.
The aleatoric-parametric distribution shown in Fig. 12 shows the impact of the parametric uncertainty when it is combined with the aleatoric uncertainty inherent in the neutron population growth rate. The aleatoric-parametric probability distribution has mean and standard deviation 205s (13.4s) so the effect of the parametric uncertainty is to make the probability distribution broader and shift it slightly in favour of shorter wait-times.

Delayed neutron precursor yields
The uncertainty applied to the yields of the delayed neutron precursor groups results in the parametric probability distribution shown in Fig. 13. This parametric uncertainty applied to the deterministic wait-time (without aleatoric uncertainty) results in a mean and standard deviation of 207 s (5.2 s).
The aleatoric-parametric distribution shown in Fig. 13 shows the impact of the parametric uncertainty when it is combined with the aleatoric uncertainty inherent in the neutron population 2.23 2:68 Â 10 À3 0.007 2:60 Â 10 À3 2:76 Â 10 À3 5 0.496 8:44 Â 10 À4 0.008 7:52 Â 10 À4 9:35 Â 10 À4 6 0.179 1:71 Â 10 À4 0.003 1:37 Â 10 À4 2:06 Â 10 À4 Fig. 11. Deterministic wait-time and aleatoric wait-time probability distributions for a cumulative detector count corresponding to a fission rate of 2 Â 10 9 fissions per second in the Caliban experiment.  growth rate. The aleatoric-parametric probability distribution has the same mean value of 211 s as the aleatoric-only distribution but the standard deviation is increased slightly from 12.8 s to 13.8 s. The uncertainty in the measurement of the delayed neutron precursor yields therefore has a small but non-negligible impact on the overall wait-time probability distribution.

Summary of results
The mean wait-times and standard deviations discussed in this section are summarised in Table 5.

Conclusions
A non-intrusive method has been demonstrated for quantifying the impact of parametric uncertainty on the probability distribution for the wait-time in a delayed supercritical system. The method makes use of the computational efficiency of the gamma distribution method to determine the wait-time probability distribution which makes uncertainty quantification feasible via the Monte Carlo method.
This method of uncertainty quantification has been applied to the criticality accident that occurred in 1958 at the Y-12 National Security Complex in Oak Ridge, Tennessee. Uncertainty in the flow rate of liquid into the drum was modelled by varying the reactivity insertion rate, and the resulting uncertainty in the predicted waittime probability distribution was evaluated. Simulations were run with two different values of the intrinsic neutron source strength in order to demonstrate the sensitivity of the wait-time probability distributions (aleatoric, parametric and combined aleatoricparametric) to this parameter.
The method has also been applied to a model of an experiment on the Caliban reactor where two areas of uncertainty were examined. In the first case, the method was used to quantify the impact of uncertainty in the timing of the movement of the reactor safety block on the predicted wait-time probability distribution. The results show significant sensitivity of the wait-time probability distribution to the timing of the safety block movement. In the second Caliban case, the model was used to quantify the impact of epistemic uncertainty in the yields of delayed neutron precursor groups on the predicted wait-time probability distribution. The results show that reported uncertainty in the measurement of the delayed neutron precursor group yields is sufficient to have a small but significant impact on the wait-time probability distribution.
The range of applicability within which the gamma distribution method can be relied upon to produce an accurate prediction of the wait-time probability distribution has been examined. A range of applicability for a reactivity insertion in a hypothetical system is presented in terms of the magnitude and rate of the reactivity insertion, the strength of the neutron source and the lifetime of prompt neutrons in the system.
Wait-time probability distributions predicted by the gamma distribution method have been verified against predictions made using the more rigorous saddlepoint method. A method for rapid solution of the saddlepoint method using Newton-Raphson iterations to accelerate convergence has been demonstrated.

Data statement
In accordance with EPSRC funding requirements all supporting data used to create figures in this paper may be accessed at the following URL: http://doi.org/10.5281/zenodo.1215970.