Algorithm for Calculating Uncertainty in the Computational Simulation for Radiochronometry of Nuclear Materials

핵물질 연대추정을 위한 전산모사 불확도 계산 알고리즘

Abstract

Nuclear forensics is an essential part of nuclear material control and nuclear non-proliferation verification. Radiochronometry for nuclear forensics is used to estimate the timing of refinement and production of nuclear materials based on decay chain characteristics and the Bateman equation. The results of radiochronometry have uncertainties because the decay constant and number of nuclides are statistics derived from analyses or repeated experiments and involve uncertainties. The aim of this study was to develop an uncertainty calculation algorithm by performing computational simulation to overcome the limitations of the existing uncertainty calculation method for radiochronometry based on the Bateman equation. The results of the proposed uncertainty calculation algorithm were comparable to those of the existing method. The algorithm allowed for more than two generations of uncertainty calculations and mitigated the underestimation of the decay constant during the uncertainty calculation.

핵감식은 국제사회에서 핵물질 통제 및 핵비확산 검증에 필수적인 부분으로 인식되고 있다. 핵감식을 위한 연대추정은 붕괴계열 특성 및 베이트만 방정식을 기반으로 핵물질의 정제 및 생산시기를 추정한다. 연대추정을 위한 요소 중 붕괴상수와 핵종 수는 분석이나 반복 실험을 통해 도출된 통계로 불확도를 가지기 때문에 연대추정의 결과도 불확도를 가진다. 본 연구는 기존의 베이트만 방정식 기반의 연대추정 불확도 계산의 한계를 극복하기 위해 전산모사를 통한 불확도 계산 알고리즘을 연구하였다. 불확도 계산 결과 기존 불확도 계산 방법과 동등한 수준의 결과가 나타났으며, 기존의 한계를 극복하여 2세대 이상의 불확도 계산이 가능하였고 불확도 계산 중 붕괴상수의 과소평가도 개선되었다.

Ⅰ. INTRODUCTION

Radiochronometry is the measurement of radioactive materials and their decay products to determine the age of the material. Radiochronometry of nuclear materials involves uncertainties and affects the reliability of the results. Therefore, the uncertainty in nuclear activity radiochronometry requires an appropriate calculation that is neither overestimated nor underestimated as a representative value of the accuracy and reliability of the radiochronometry results. However, the calculation method for uncertainty estimation proposed by the LLNL(Lawrence Livermore National Laboratory) and others is dominated by the uncertainty in the activity ratios of parent and daughter nuclides. Therefore, it underestimates the uncertainties in decay constants and initial nuclides. By contrast, the calculation method for estimating the uncertainty of nuclear activity radiochronometry by applying the Bateman equation can only estimate the uncertainty in the one-step relationship between the parent and the daughter nuclides owing to the complexity of the equation. Therefore, problems with the reliability of the radiochronometry results, in terms of the uncertainty and the inability to select different signature nuclides, are encountered. The LLNL proposed that the development of an uncertainty estimation algorithm utilizing Monte Carlo methods was the only approach for uncertainty calculations that included decay constants and other factors to mitigate the underestimation of uncertainty factors, which was a limitation of the calculations for uncertainty estimation based on the Bateman equation, and to estimate uncertainty for signature nuclide relationships in decay chain of more than two steps, such as parent and daughter nuclides^[1,2].

The aim of this study is to prevent the underestimation of uncertainty factors in radiochronometry calculations by utilizing computational simulation methods. Additionally, the aim is to develop an algorithm that allows for uncertainty estimation even when selecting signature nuclides relationships beyond the two step within the decay chain.

The calculation and computational simulation methods for uncertainty estimation based on the Bateman equation are investigated to develop an algorithm that can estimate uncertainty even when selecting signature nuclides in more than two steps of relationships within a decay chain. The results of the developed computational simulation method are compared with the calculation results of the conventional radiochronometry method based on the Bateman equation. The developed algorithm for estimating uncertainty in radiochronometry can be utilized as a nuclear-forensics-based technology to enhance nuclear non-proliferation capabilities globally and strengthen nuclear material control capabilities.

Ⅱ. MATERIAL AND METHODS

1. Development of an algorithm for estimating uncertainty in nuclear material radiochronometry

Uncertainty is a non-negative parameter that characterizes the variance of a measurement over a quantity based on the information used^[3]. Uncertainty is categorized into two types based on how it is evaluated: Types A and B. Type A uncertainty is defined as the uncertainty calculated by statistical analysis of the observations, and Type B uncertainty is defined as the uncertainty calculated by methods other than statistical analysis of the observations ^[3]. The conventional estimation of uncertainty in radiochronometry based on the Bateman equation is Type B uncertainty, which is calculated by arithmetic combined standard uncertainty. However, methods based on the Bateman equation can only derive uncertainties up to the first generation of the relationship and are limited by the underestimation of the decay factor. This has led to research on methods for calculating Type A uncertainty, which is derived from statistical analysis. This study analyzed the mathematical method based on the Bateman equation for Type B uncertainty and investigated the computational simulation method for Type A uncertainty estimation.

2. Estimating uncertainty with mathematical models

In this study, the uncertainty propagation model and influencing factors were analyzed to investigate the uncertainty estimation method using a mathematical model based on the Bateman equation, and the uncertainty estimation algorithm using the mathematical model and influencing factors was derived.

2.1. Uncertainty propagation model and influencing factors

Nuclear activity radiochronometry utilizes the Bateman equation^[4], which is a system of coupled differential equations for variables, such as the initial abundance of each nuclide, branching ratio, and decay constant, to represent the abundance of nuclides over time when multiple nuclides decay simultaneously and is expressed as Equation (1).

\(\begin{aligned}N_{n}(t)=N_{1}(0) \times\left(\prod_{i=1}^{n-1} \lambda_{i}\right) \times \sum_{i=1}^{n} \frac{e^{-\lambda_{i} t}}{\prod_{j=1, j \neq i}^{n}\left(\lambda_{j}-\lambda_{i}\right)}\end{aligned}\)       (2)

N : Number of nuclide

λ : Decay constant [y^-1]

t : Time [y]

For nuclear activity radiochronometry, if the initial nuclear material (parent nuclide) is assumed to be in a pure state unmixed with daughter nuclides at the time of production (or purification), the following equation holds between the parent and daughter nuclides^[5].

\(\begin{aligned}N_{2}(t)=\frac{\lambda_{1}}{\left(\lambda_{2}-\lambda_{1}\right)} N_{1}(t)\left(1-e^{\left(\lambda_{1}-\lambda_{2}\right) t}\right)\end{aligned}\)       (2)

Solving Equation (2) for time t leads to Equation (3), which can be used to estimate the timing of the production (or purification) of the nuclear material under investigation^[1].

\(\begin{aligned}t=\frac{1}{\left(\lambda_{1}-\lambda_{2}\right)} \ln \left[1+\frac{N_{2}(t)}{N_{1}(t)} \frac{\left(\lambda_{1}-\lambda_{2}\right)}{\lambda_{1}}\right]\end{aligned}\)       (3)

Equation (3) is used for the calculation of the elapsed time, t, from the main factors of the Bateman equation, t, N, and λ. The main influencing factors for calculating t in Equation (3) are N, which is the uncertainty in the analysis of the nuclear material, and λ, which is the decay constant and has the uncertainty of statistics derived from repeated experiments. The uncertainty of N varies depending on the material and method, but the analytical uncertainties of ²³⁴U and ²³⁰Th at each of the laboratories utilized in this study were within 1%^[6]. The uncertainty of λ is smaller than the uncertainty of N, with an uncertainty of approximately 0.1% for most nuclides. However, according to the half-life information provided by the National Nuclear Data Center, it has a high uncertainty of approximately 15% for nuclides with short half-lives, such as ²¹⁸At, which can be an important uncertainty factor in age estimation using the characteristics of the decay chain. Thus, as N and λ are uncertain values, the radiometric value t derived from those factors is uncertain^[7].

The N derived from gamma spectroscopy, alpha spectroscopy, and mass spectrometry can be expressed as R, which is the nuclide ratio of the parent to daughter nuclides. The uncertainties (u_c,γ, u_c,α, and u_c,ICP-MS) from gamma, alpha, and ICP-MS analyses propagate into the uncertainty (u_R) of R, which propagates into the uncertainty in radiochronometry (u_t), similar to the uncertainties (u_λ1, u_λ2) due to decay constants. Fig. 1 shows the schematic of the propagation model.

Fig. 1. Radiochronometry uncertainty propagation model.

2.2. Method of estimating uncertainty with the mathematical model

The calculation method for estimating uncertainty in radiochronometry using the Bateman equation is calculated using Equation (3), which is derived from Equations (1) and (2) in Section 2.1.1, and is shown in Equation (4) with λ and R as the influence factors for calculating the uncertainty estimate.

\(\begin{aligned}t=\frac{1}{\left(\lambda_{1}-\lambda_{2}\right)} \ln \left[1+R \frac{\left(\lambda_{1}-\lambda_{2}\right)}{\lambda_{1}}\right]\end{aligned}\)       (4)

R : Activity ratio

The uncertainty in radiochronometry is calculated in Equation (5) by using the sensitivity coefficients derived by the partial differentiation of Equation (4) for each factor and the factor-specific uncertainties.

\(\begin{aligned}u(t)=\sqrt{\left(u\left(\lambda_{1}\right) \frac{\partial t}{\partial \lambda_{1}}\right)^{2}+\left(u\left(\lambda_{2}\right) \frac{\partial t}{\partial \lambda_{2}}\right)^{2}+\left(u(R) \frac{\partial t}{\partial R}\right)^{2}}\end{aligned}\)       (5)

To derive Equation (5) for calculating the uncertainty in radiochronometry, Equation (4) is partially differentiated for each factor, and the calculated sensitivity coefficients for each factor are shown in Equations (6) - (8). The contribution rates for each factor in Equation (5) are calculated in Equations (9) - (11)^[8].

\(\begin{aligned}\frac{\partial t}{\partial R}=\frac{1}{\lambda_{1}+R\left(\lambda_{1}-\lambda_{2}\right)}\end{aligned}\)       (6)

\(\begin{aligned}\begin{array}{c}\frac{\partial t}{\partial \lambda_{1}}=\left(\frac{1}{\lambda_{1}-\lambda_{2}}\right)\left(\frac{1}{1+\frac{R\left(\lambda_{1}-\lambda_{2}\right)}{\lambda_{1}}}\right)\left(\frac{\lambda_{2} R}{\lambda_{1}^{2}}\right) \\ -\frac{1}{\lambda_{1}-\lambda_{2}} \ln \left(1+\frac{R\left(\lambda_{1}-\lambda_{2}\right)}{\lambda_{1}}\right)\end{array}\end{aligned}\)       (7)

\(\begin{aligned}\begin{array}{l} \frac{\partial t}{\partial \lambda_{2}}= \frac{1}{\lambda_{1}-\lambda_{2}} \ln \left(1+\frac{R\left(\lambda_{1}-\lambda_{2}\right)}{\lambda_{1}}\right) \\ -\frac{R}{\lambda_{1}\left(\lambda_{1}-\lambda_{2}\right)\left(1+\frac{R\left(\lambda_{1}-\lambda_{2}\right)}{\lambda_{1}}\right)}\end{array}\end{aligned}\)       (8)

Equations (6) - (8) are the sensitivity coefficients for each factor (R, λ₁, and λ₂) to the uncertainty in radiochronometry and are used in conjunction with the uncertainty of each factor to derive the mathematical modeling of Equation (5) for the uncertainty in radiochronometry.

The contribution of each factor (R, λ₁, and λ₂) to the uncertainty in radiochronometry is given by Equations (9) - (11).

\(\begin{aligned} \begin{array}{l} Contribution \; of \;R=\\\frac{\left(u(R) \frac{\partial t}{\partial R}\right)^{2}}{\left(u\left(\lambda_{1}\right) \frac{\partial t}{\partial \lambda_{1}}\right)^{2}+\left(u\left(\lambda_{2}\right) \frac{\partial t}{\partial \lambda_{2}}\right)^{2}+\left(u(R) \frac{\partial t}{\partial R}\right)^{2}}\end{array}\end{aligned}\)       (9)

\(\begin{aligned} \begin{array}{l} Contribution \; of \; \lambda_{1}=\\\frac{\left(u\left(\lambda_{1}\right) \frac{\partial t}{\partial \lambda_{1}}\right)^{2}}{\left(u\left(\lambda_{1}\right) \frac{\partial t}{\partial \lambda_{1}}\right)^{2}+\left(u\left(\lambda_{2}\right) \frac{\partial t}{\partial \lambda_{2}}\right)^{2}+\left(u(R) \frac{\partial t}{\partial R}\right)^{2}} \end{array}\end{aligned}\)       (10)

\(\begin{aligned}\begin{array}{l} Contribution \; of \; \lambda_{2}=\\\frac{\left(u\left(\lambda_{2}\right) \frac{\partial t}{\partial \lambda_{2}}\right)^{2}}{\left(u\left(\lambda_{1}\right) \frac{\partial t}{\partial \lambda_{1}}\right)^{2}+\left(u\left(\lambda_{2}\right) \frac{\partial t}{\partial \lambda_{2}}\right)^{2}+\left(u(R) \frac{\partial t}{\partial R}\right)^{2}}\end{array}\end{aligned}\)       (11)

3. Estimating uncertainty with the Monte Carlo method

The calculations for estimating uncertainties in nuclear activity radiochronometry using the Bateman equation are described in Sections 2.1. and 2.2. have the disadvantage that they are valid only for certain conditions, such as the parent and daughter nuclide relationships, and cannot be applied to the selection of daughter nuclides and further relationships due to the complexity of the derivative calculations.

In this study, the Monte Carlo method for uncertainty estimation was developed based on the computational simulation to overcome the limitations of the uncertainty estimation method of the mathematical modeling of the Bateman equation. Uncertainty estimation methods and influencing factors were analyzed, and an uncertainty estimation algorithm using the LHS(Latin hypercube sampling) method was derived for reliable and efficient computational simulation.

3.1. Monte Carlo uncertainty estimation method and influencing factors

The Monte Carlo method utilizes the property in which the uncertainty factors, λ and R, derived from the Bateman equation, are values with uncertainty based on the uncertainty distribution. As the results are derived by calculating the activity ratio of the signature nuclide and estimating the activity ratio appropriate to the analyzed value of R, the Monte Carlo method enables uncertainty estimation by calculating the uncertainty factor λ used to calculate the activity ratio and the uncertainty factor R for radiochronometry to randomly vary based on the uncertainty distribution, deriving values through an iterative process and analyzing the distribution of the calculated results. However, the Monte Carlo method requires many samples to ensure high reliability because statistical errors different from the existing probability distribution occur in a small number of samples. Problems are associated with the computing time and efficiency related to sample collection and analysis because the complex calculation process of radiochronometry must be carried out in iterative operations involving calculations for each sample.

3.2. Latin hypercube sampling (LHS) algorithm for estimating uncertainty

An efficient method for extracting the variable is needed to prevent statistical errors when the statistical probability distribution analyzed for an unknown variable with a mean of \(\begin{aligned}\bar x\end{aligned}\) is known and the number n of samples extracted from this probability distribution is small. The LHS method, proposed by Mckay, Beckman, and Conover (1979), is a method for obtaining effective statistical results with fewer samples than that through random sampling for computational simulations using the Monte Carlo method^[9]. LHS is a method of dividing the range of the probability distribution into n to prevent biased samples from being extracted from the entire probability distribution and then extracting n from each of the following intervals one by one to avoid overlapping extraction. For a uniform distribution, as shown in Fig. 2, the probability density function equally divides each bin, and samples are drawn from the equally divided bins, as shown in the cumulative distribution function^[10].

Fig. 2. Probability density function (left) and Cumulative distribution function (right) in uniform distribution.

When expanded to two dimensions (k = 2) with two variables, the grid is divided into equal probability intervals equal to the number n^k, and random sampling is performed such that each row and column is selected only once, allowing the entire range of each variable to be evenly sampled. Fig. 3 displays that when two variables follow a uniform distribution between 0 and 1, if a sample of size 5 is drawn to the LHS, each variable is drawn between 1 and 5 to obtain [5, 3, 2, 1, 4] and [1, 3, 2, 5, 4], and the position of each grid is determined when the two variables are paired. Only one point was drawn from the corresponding interval for each row and column, indicating that the values of the variables were sampled at the same rate across the entire range of the probability distribution.

Fig. 3. Sampling example from LHS method.

To utilize LHS in the case of a normal distribution, we define the range of the normal distribution as [F, G], divide it into n intervals between F and G, and extract n samples, one from each interval, to avoid overlapping samples. For a normal distribution (Fig. 4), if the area of each interval, representing the sum of the probabilities of being extracted from the probability density function (extraction probability), is equally divided, the equal area of the probability density function means that the probability distribution range is divided by the same extraction probability. This is equivalent to evenly splitting the y-axis because the y-axis of the cumulative distribution function represents the cumulative sum of the probabilities. An LHS with a normal distribution can be sampled by extracting the y-axis from these evenly divided probability intervals and finding the x-axis corresponding to the extracted y-axis from the cumulative distribution function.

Fig. 4. Probability density function (left) and Cumulative distribution function (right) in normal distribution.

Random number generators used by most computer programs first generate random numbers that follow a uniform distribution because it is relatively simple and efficient to generate uniformly distributed random numbers. Probability distributions are transformed using the aforementioned inverse transformation process to extract random numbers from probability distributions other than the uniform distribution.

Fig. 5. Form of probability density function (left), Cumulative distribution function (middle), and Inverse cumulative distribution function (right) in normal distribution

This is the inverse transformation method, which uses the inverse cumulative distribution function to convert random numbers generated from a uniform distribution into a different probability distribution and generate random numbers with the desired probability distribution. Therefore, the horizontal axis of the inverse cumulative distribution function is divided into n equal widths between 0 and 1, and random numbers are extracted for each interval to change the probability distribution using LHS in the program, as shown in Fig. 6. Then, through the inverse cumulative distribution function, extracted sample values, which are converted to the probability distribution to apply, can be obtained by pairing each with the y-axis^[11].

Fig. 6. Transform using Inverse cumulative distribution function.

In this study, an algorithm using the LHS method was implemented in Matlab computer software (MathWorks, USA, R2023a Update2)^[12] for computational simulation for random sampling owing to the characteristics of the Monte Carlo method. The radiochronometry algorithm in the program proceeds as shown in Fig. 7 and is divided into two main steps: radioactivity calculation and radiochronometry calculations.

Fig. 7. Computer simulation uncertainty estimation process.

First, the value with uncertainty in Equation (1), which is derived on the Bateman equation, is the constant λ, and each time a set of samples is calculated, the value of λ within the uncertainty range is randomly sampled by the LHS method to replace the previous constant value of λ to calculate the activity ratio of the two nuclides selected as signature nuclides over time. Then, time t was estimated using the "vpasolve()" function^[13], which finds an approximate solution to a given equation using input values within Matlab, given an activity ratio, R, analyzed for the selected radionuclide. The activity ratio, R, is also a value with uncertainty, and each time one set of samples is calculated. The R value within the uncertainty range is randomly sampled through the LHS method, and the radiochronometry results corresponding to the set number of sample sets are collected and statistically analyzed to calculate the new standard deviation and variance to derive the uncertainty.

Ⅲ. RESULT

1. Nuclear material radiochronometry uncertainty model analysis

The mathematical model presented in Section 2 was evaluated. The uncertainty estimation model was examined by using the Monte Carlo and LHS methods, and the uncertainty estimation results were analyzed by using the algorithm.

1.1. Mathematical method uncertainty model analysis

The sensitivity of the calculations for uncertainty estimation based on the Bateman equation was analyzed. Analysis of the sensitivity equations (6) - (8) revealed that each sensitivity coefficient shows a unique change depending on the sensitivity factor, as shown in Fig. 8. The common factor is that the sensitivity coefficient infinitely increases in the radial equilibrium interval.

Fig. 8. Sensitivity coefficient and radial equilibrium graph.

As a result of the change in the contribution rate using Equations (9) - (11) assuming an uncertainty of 1% for each factor (R, λ₁, and λ₂), the uncertainty contribution of λ₁ is initially as high as 50%. However, over time, the uncertainty contribution of the nuclide ratio, R, becomes dominant, as presented in Table 1. The decay constants of the parent nuclide have a greater impact than the daughter nuclide, and the contribution of R increases with time.

The reason for the temporary increase in λ₂ at 100,000 years in Table 1 is that the sensitivity coefficient of each factor increases as it approaches the radial equilibrium interval, as depicted in Fig. 8. As λ₁ shows a sharp change after approaching the equilibrium activity ratio and λ₂ slowly increases, the temporary rise may have been caused by the preceding effect of increasing sensitivity. The reason for the sharp change in R and λ₁ when changing from 100,000 to 1,000,000 years in Table 1 is that the sensitivity coefficient infinitely increases toward the radial equilibrium interval, as shown in Fig. 8. As the rate of increase in R is relatively higher, the decay constant is underestimated by the synthetic uncertainty characteristic because the error synthesis converges to a large factor in the uncertainty estimation calculation based on the Bateman equation. As the uncertainty in the decay constant is 1% or less and is a fixed analytical value, the uncertainty in the nuclide ratio, R, of the signature nuclide from the actual analysis is determined to be an important factor in determining the uncertainty in radiochronometry.

Table 1. Contribution rate by Radiochronometry uncertainty factor(²³⁰Th/²³⁴U)

The uncertainty in radiochronometry based on Equation (5) depends on the radial equilibrium properties of the two nuclides chosen for signature nuclides, with nuclides that take a long time to equilibrate maintaining a low uncertainty for a long period of time, while nuclides that rapidly equilibrate are characterized by a spike in the uncertainty in a short period of time. The left graph in Fig. 9 shows the uncertainty evolution with ²²⁸Th/²³²U as the signature nuclide, with uncertainties exceeding 50% at times of 20 years or less because radial equilibrium appears quickly, whereas the uncertainty is less than 10% over approximately 100,000 years with ²³⁰Th/²³⁴U as the signature nuclide.

Fig. 9. Uncertainty change graph according to signature nuclides.

As the half-life uncertainty factor characteristics and radial equilibrium characteristics of nuclides differ depending on the signature nuclide, the sensitivity, contribution rate, and radiochronometry range of each uncertainty estimation contributor should be considered in the radiochronometry uncertainty estimation.

1.2. Monte Carlo and LHS method uncertainty model analysis

To compare the reliability of the Monte Carlo method and the LHS, the convergence of each method to the original distribution when sampling is compared when the probability distribution is normal, as shown in Fig. 10^[14-16]. The left side of Fig. 10 shows a comparison of the Monte Carlo and LHS methods for 100 samples, with the Monte Carlo method showing less convergence to the original distribution to the point where the original distribution is unrecognizable when the number of samples is small. On the right is a comparison of the Monte Carlo and LHS methods for 10,000 samples, and while both methods are normally distributed, the LHS method is relatively more convergent to the normal distribution than the Monte Carlo method. The convergence rate of the Monte Carlo method and the LHS method theoretically approaches the classical distribution as the sample size increases. However, studies have shown that the number of samples (n) squared of the number of samples (n²) in the LHS method is required to achieve accuracy in the Monte Carlo method using the LHS method^[17].

Fig. 10. Comparison of convergence between MC and LHS methods in normal distribution. (100 samples left and 10,000 samples right)

The standard deviation of the uncertainty estimation results as a function of the number of samples (iterations) of the LHS method is analyzed to examine the convergence of the distribution as the number of samples (iterations) of the computational simulation method increases. Table 2 lists the radiochronometry results and the standard deviation of the results as a function of the number of iterations (number of samples) using the LHS method. For radiochronometry, the signature nuclides ²³⁰Th/²³⁴U, an activity ratio of 0.001, a nuclide analysis uncertainty of 10% for ²³⁰Th and ²³⁴U, and an age of 109.208 years were used.

Table 2. Standard deviation of uncertainty results according to the number of repetitions

As the number of samples increases, as shown in Table 2, the results are closer to a uniform value with a small standard deviation, which can be plotted in Fig. 11, showing that the standard deviation logarithmically decreases as the number of iterations increases. Increasing the number of iterations has a large impact on accuracy in intervals where the standard deviation rapidly decreases, but the impact decreases in later intervals.

Fig. 11. Variation of standard deviation according to the number of repetitions.

The computation time increases as the number of samples (iterations) increases, such the appropriate number of iterations should be set for efficient uncertainty estimation. As described in Section 3.1.1, whether the same tendency of uncertainty evolution with time was observed in the method using computational simulated random sampling was analyzed, where the uncertainty infinitely increased as the activity ratio of the signature nuclide approached radiometric equilibrium in the mathematical model^[1].

Fig. 12 shows a graph of the uncertainty curve of the radiochronometry results derived from the LHS method and the mathematical modeling method as a function of the change in activity ratio. For radiochronometry, the analysis uncertainty of signature nuclides ²³⁰Th/²³⁴U, ²³⁰Th, and ²³⁴U was set at 1% each, and the number of samples was set at 200. As shown in Fig. 12, both the LHS and mathematical model methods showed the same tendency as the equations and algorithms used in the radiochronometry process. As the radial equilibrium ratio increased and approached radial equilibrium, the uncertainty of the LHS method became relatively higher compared to the mathematical modeling method. This suggested that the impact of the uncertainty in the decay constant, which was underestimated in the previous mathematical modeling method, was appropriately reflected in the uncertainty. However, the maximum activity ratio of the radial equilibrium of the signature nuclides ²³⁰Th/²³⁴U was more than 1.44, allowing the mathematical model to estimate the uncertainty of the equilibrium up to the approximation of the equilibrium. However, for the LHS method, more than 1.42 caused an algorithm error. As for the cause of the error, when the activity ratio R is a value with uncertainty and a value of 1.42 or more is entered in the process of sampling according to the uncertainty distribution, if a value of 1.44 or higher is sampled, the problem of not being able to find a solution occurs as it is outside the maximum activity ratio of the selected signature nuclide, requiring additional research to solve future problems.

Fig. 12. Uncertainty changes according to changes in radioactivity ratio of LHS method (left) and mathematical model method (right).

2. Validation of nuclear material radiochronometry uncertainty estimation results

In this study, the uncertainty estimation results of the LHS method were compared with the uncertainty estimation results of the mathematical model and radiochronometry proposed by international institutions to verify the uncertainty estimation results of the computer-simulated random sampling method. The uncertainty estimation results for the selection of signature nuclides in the two-generation and aforementioned relationships were analyzed.

2.1. Comparative analysis of uncertainty estimation results

The results of four institutions, LLNL (USA), Los Alamos National Laboratory (LANL, USA), Japan Atomic Energy Agency (JAEA, Japan), and Commissariat a l'Energie Atomique et aux Energies Alternatives (CEA, France), with the information of ²³⁴U and ²³⁰Th, and the results of the estimation of the radiochronometry uncertainty in radiochronometry calculated by employing the mathematical model and LHS method were compared and shown in Table 3 to evaluate the results of the calculation of the uncertainty in radiochronometry^[6]. The measurement uncertainties of ²³⁴U and ²³⁰Th were calculated using the measurement uncertainties of ²³⁴U and ²³⁰Th analyzed by each institution and the analytical results. The comparison shows that the results of the mathematical modeling and uncertainty estimation implemented by the LHS method are close to each other, as shown in Table 3. The expected differences are due to differences in decimal place handling and differences in the decay constants and uncertainties of the decay constants^[18], which are the eigenvalues of the nuclides applied in the calculations by each institution.

Table 3. Comparison of the calculation results of Radiochronometry uncertainty using ²³⁰Th and ²³⁴U

1) Half-life using mathematical model and LHS methods : λ²³⁴U = 2.8234×10^-6/year, λ²³⁰_Th = 9.1954×10^-6/year ^[7]

2) Half-life of use by institution : λ²³⁴_U = 2.8263×10^-6/year, λ²³⁰_Th = 9.1580 ×10^-6/year

3) CEA's Radiochronometry estimation uncertainty is expressed up to the first decimal place

4) JAEA's Radiochronometry estimation uncertainty is expressed as an average value

5) Number of repeated LHS methods (number of samples): 200

6) Comparative data: Table 1 of Reference ^[6] [Appendix 1]

2.2. Analysis of uncertainty estimation results for two generations or more

In this study, the LHS method was used to analyze the change in uncertainty when performing radiochronometry by selecting signature nuclides that differed by more than two generations in the decay chain, which were difficult to calculate using mathematical model equations. Table 4 lists the change in uncertainty for each generation by selecting signature nuclides (²²⁶Ra, ²²²Rn, and ²¹⁸Po) within the ²³⁴U decay chain that differ by more than two generations and entering an activity ratio that is approximately 100 years old. For radiochronometry, the analysis uncertainty of the ²³⁴U decay chain, ²³⁴U, and selected nuclides was set at 1% each, and the number of samples (n) was set at 500. As shown in Table 4, the half-life uncertainty of the calculation process is reflected as the related generation increases, confirming that the uncertainty result of each generation increases as the generation increases for the same radiochronometry t. However, uncertainty estimation beyond the second generation is currently lacking comparative validation data, and further research is needed for validation.

Table 4. Uncertainty of Radiochronometry of 2 generations or more

Ⅳ. DISCUSSION

In this study, the uncertainty estimation methods for nuclear material and nuclear activity radiochronometry were investigated by using mathematical modeling and LHS methods to estimate the uncertainties that should be considered during the radiochronometry process and improve the reliability and accuracy of nuclear activity radiochronometry. A mathematical modeling method based on the Bateman equation for estimating the uncertainty in radiochronometry using the decay chain of nuclear materials such as uranium or spent nuclear fuel was analyzed. An LHS method based on the Monte Carlo method of computational simulation with random sampling was developed to overcome the limitation of not being able to date more than one generation.

The characteristics of the calculation equation and sensitivity coefficient for uncertainty estimation of the Bateman equation-based mathematical model were determined. Each uncertainty component in the mathematical model increased in coefficient with time and diverged to infinity as it approached the radial equilibrium interval.

As the time to reach the radial equilibrium is different for each nuclide, the appropriate signature nuclide must be selected based on the nuclear material conditions for radiochronometry. Furthermore, the uncertainty in the decay constant of the parent nuclide and the nuclide ratio of the signature nuclide are the most influential factors in radiochronometry uncertainty, and the contribution rate of the uncertainty in the nuclide ratio increased as time increased. For the computational simulation sampling method, the Monte Carlo method and the LHS method were compared to derive the utility and characteristics of the LHS method. The LHS method is highly reliable with a small number of samples compared to the Monte Carlo method, and the standard deviation logarithmically decreased as the number of samples increased since the LHS method is also based on the Monte Carlo method. The LHS method also showed the same spike in uncertainty in the radial equilibrium interval as the mathematical model.

To compare the radiochronometry calculations, the results of the mathematical modeling and LHS methods were compared with those of the LLNL(Lawrence Livermore National Laboratory), LANL(Los Alamos National Laboratory), JAEA(Japan Atomic Energy Agency), and CEA(Commissariat a l'Energie Atomique et aux Energies Alternatives), and they showed similar performance in estimating the uncertainty in radiochronometry, even for two or more generations of signature nuclide relationships. However, It appears that there are factors beyond the differences in the accuracy of input values and differences in decimalization of output values are present. The computational simulation methods, when compared to the conventional Bateman equation-based approach, have reduced the underestimation of uncertainty factors such as decay constants. Furthermore, while the conventional Bateman equation-based approach faced challenges in estimating radiochronometry uncertainties in relationships beyond two steps due to the complexity of uncertainty formulas, the computational simulation methods enable the estimation of uncertainties in radiochronometry for relationships beyond two steps. Further research is needed to compare and verify the calculation results that account for complex decay chains and non-homogeneous conditions, reflecting the characteristics of realistic sample types and nuclear forensics environments.

The results of this study can be used as a research and development tool for identifying signature nuclides and developing related technologies in the field of nuclear forensics, as well as a tool for comparing and verifying the calculation results of operators during independent verification for spent nuclear fuel management in the future. Finally, it can be used as a basic technology in nuclear forensics for strengthening nuclear material control capabilities at domestic nuclear facilities and for independent and active participation in nuclear non-proliferation verification in neighboring countries.

Ⅴ. CONCLUSION

This study conducted research on a computational simulation-based uncertainty calculation algorithm for radiochronometry of Nuclear Materials. Radiochronometry uncertainty increased over time and diverged to infinity as it approached the radial equilibrium interval. The Latin Hypercube Sampling method exhibited higher reliability with fewer samples compared to the Monte Carlo method. The computational simulation methods, when compared to the conventional Bateman equation-based approach, have reduced the underestimation of uncertainty factors such as decay constants. Furthermore, while the conventional Bateman equation-based approach faced challenges in estimating radiochronometry uncertainties in relationships beyond two steps due to the complexity of uncertainty formulas, the computational simulation methods enable the estimation of uncertainties in radiochronometry for relationships beyond two steps. In the future, this study can be applied as a foundational tool in the field of nuclear forensics, addressing complex decay chains and non-homogeneous conditions while reflecting the characteristics of realistic sample types and nuclear forensic environments.