Simulating Three-Dimensional Plume Migration of a Radionuclide Decay Chain through Groundwater

In this study, we present a semi-analytical model for simulating three-dimensional radioactivity transport of a radionuclide decay chain and assessing the radiological dose impact on the general public. The mechanisms and processes considered in the model include the one-dimensional advection, hydrodynamic dispersion in longitudinal and two lateral directions, linear equilibrium sorption, and first-order radioactive decay reactions. The semi-analytical model is derived for a semi-infinite domain, and the solutions in the Laplace domain for members of the decay chain are first generalized in a compact format. The concentrations in the original domain of each nuclide are independently evaluated with the help of the efficient and robust Laplace numerical inverse algorithms. The accuracy of the derived semi-analytical model is demonstrated by comparison of our developed model with an existing analytical model described in the literature. The results of the verification exercise indicate that the derived semi-analytical model is accurate and robust. The developed semi-analytical model is applied to an illustrative example that simulates the three-dimensional plume migration of a radionuclide decay chain on both the temporal and spatial scales. Moreover, the time histories of the radiological doses at different distances from the inflow source boundary are presented to understand the potential radiological impact on the general public. The developed model facilitates rapid assessment of the radiological impact posed by the presence of radionuclides in the environment because of leakage from a nuclear waste repository or accidental discharges from nuclear facilities.


Introduction
Nuclear technology has been extensively employed for the development of nuclear weapons, production of electricity, medical and industrial applications, as well as a variety of other civilian uses for approximately seventy years. Although society has greatly benefited from the use of nuclear technology, there exist potential risks posed by the presence of radionuclides in the environment, whether released from nuclear waste repositories or through accidental discharges from various nuclear facilities including nuclear power plants, fuel reprocessing plants, and uranium mining and milling operations. It is important to build a reliable emergency and long-term risk management system for assessing the radiation dose received by the people. Groundwater is one of the most important environmental media for the transport of radionuclides after their release from waste disposal areas or from certain classes of accidental releases from nuclear facilities. Actual observation of long-term transport behavior of radionuclides released from the waste disposal areas cannot be performed however. Groundwater transport models along with experiments, measurements, and observations performed in the field provide efficient means to calculate the expected radioactivity concentration of radionuclides following their release into the groundwater. The groundwater transport model for simulating the long-term radioactivity transport can be developed by solving the transport governing equations using analytical solutions or numerical approaches. Despite their inherent limitation to homogeneous media with relatively simple and regular geometries, analytical models are more rapid and efficient tools for predicting the long-term concentration of radioactivity at specific points as compared with time-marching numerical models which generally require an extremely large number of time steps and place severe demands on computational time. The analytical groundwater transport model can be used for certain types of analyses where the currently available data do not warrant a more complicated study. Such models may frequently be adequate for regulatory needs provided the model parameters are chosen conservatively.
Numerous analytical groundwater transport models have been reported in the literature for the simplified case of unidirectional groundwater flow with one-, two-, and three-dimensional hydrodynamic dispersion in saturated and homogeneous geological media . However, such analytical groundwater transport models have mostly been derived for a single dissolved substance. Radionuclides decay into other radioactive products or stable species, called daughter species or progeny. The chain decay processes of radionuclides are particularly important for modeling the transport of actinides and transuranics. The parent-daughter interactions are neglected in the analytical groundwater transport models that account for the transport of only a single dissolved substance, thus leading to significant error in the predicted transport behaviors of all members of a radionuclide decay chain. Nair et al. [24] developed a numerical model for simulating three-dimensional radioactivity transport of a long radionuclide decay chain with consideration of distinct retardation factors for each individual nuclide and showed that the model with decay chain transport gives total effective doses about 1000 times higher than those calculated using a model without decay chain transport. This showed the importance of the inclusion of decay chain transport in modeling the radionuclide in the subsurface system.
In order to incorporate the chain decay processes of radionuclides, the advection-dispersion transport equation for each progeny must include the term for ingrowth from the preceding parent radionuclides. Analytical groundwater transport models that use a set of advection-dispersion equations sequentially coupled with radioactive decay reactions can serve as rapid tools for simultaneously simulating the long-term transport behaviors of the original species and its progeny of a radionuclide decay chain. However, only a few analytical solutions solved for coupled multiple radionuclide transport equations are currently available in the literature for evaluating the radionuclide decay chain migration. Developing analytical models for predicting the radionuclide decay chain transport generally involves cumbersome and complicated mathematical manipulations because of the need to solve a system of coupled advection-dispersion equations simultaneously.
Most of the analytical models for the decay chain transport currently reported in the literature are limited to one-dimensional advective-dispersive transport [17,[25][26][27][28][29][30][31][32][33][34][35]. Multidimensional models for decay chain transport for radionuclide decay chains would be much more suited for real-world applications. Bauer et al. [36] presented a pioneering analytical model for one-, two-, and three-dimensional decay chain transport. Montas [37] presented an analytical model for three-dimensional advection-dispersion equations coupled with first-order decay reactions for a three-species decay chain. Quezada et al. [38] extended Clement's [31] solution technique to derive analytical solutions in the Laplace domain for a decay chain with an arbitrary number of members. Sudicky et al. [39] presented a semi-analytical model for investigating the three-dimensional plume migration of a decay chain subject to a first-type inflow boundary condition, but their model was limited to four-level decay chain and did not obtain generalized expressions for arbitrary target species. It should also be noted that the first-type inflow boundary condition generally causes mass conservation errors, especially for a subsurface system with a large longitudinal dispersion coefficient [40][41][42][43]. Chen et al. [44] developed an exact analytical model for the two-dimensional plume migration of a decay chain in a finite-length system subject to a third-type inflow source boundary condition. A concise mathematical expression for an arbitrary member of a decay chain was obtained, but the solution involves the summations of three infinite series expansions. Numerical evaluations of the analytical solutions derived for a finite-length domain with a summation of infinite series expansion are always time-consuming and computationally inefficient [45][46][47].
Based on the literature review, this study develops a semi-analytical model for rapidly simulating the three-dimensional radioactivity transport of a radionuclide decay chain. The derived semi-analytical model includes a general mathematical expression for an arbitrary member of a radionuclide decay chain, making it easy to write a computer program code for implementing the computation. The novelty of the developed semi-analytical model is three-fold. First, the solutions are computationally efficient. Second, the solution can evaluate concentration for larger Peclet values. Third, the solutions are derived for three-dimensional radionuclide transport, which has more practical applications. The derived analytical model should have more practical applications for a rapid evaluation of the concentration of radioactivity on a spatial and a temporal scale as well as the radiological dose for a real-world radionuclide decay chain.

Mathematical Model
This study considers the movements of the original radionuclide and its serial progeny, sequentially coupled by the first-order radioactive decay reactions, producing a radionuclide decay chain. Assuming a homogeneous aquifer with a steady and unidirectional flow field moving along the x direction with three dispersion regions along the x, y, and z directions, as schematically described in Figure 1, the governing equations for the ith member of a radionuclide decay chain can be expressed as where C i (x, y, z, t) represents the radioactivity concentration of nuclide i [Bq/m 3 ]; x, y, z denote the three spatial coordinates [m], respectively; t is time [y]; D x , D y , and D z are the hydrodynamic dispersion coefficient in the x, y and z directions; v is the linear average seepage velocity; λ i is the radioactive decay constant of nuclide i [y −1 ]; R i is the retardation factor of nuclide i [-]; N is the total number of the nuclides involved in the radionuclide decay chain. This study considers the movements of the original radionuclide and its serial progeny, sequentially coupled by the first-order radioactive decay reactions, producing a radionuclide decay chain. Assuming a homogeneous aquifer with a steady and unidirectional flow field moving along the x direction with three dispersion regions along the x , y , and z directions, as schematically described in Figure 1, the governing equations for the ith member of a radionuclide decay chain can be expressed as  It is assumed that the aquifer is free of radionuclides before their release from the source area and can be stated as When simulating the transport of a radionuclide decay chain, the radionuclide sources are important because they govern the magnitude of the radioactivity concentrations and the shapes and sizes of the nuclide plumes. Herein, we consider a source region that contains specified amounts of the nuclide sources per cross-sectional area perpendicular to the groundwater flow direction and assume the release rate of each nuclide is proportional to the amounts of nuclides staying in the repository. A set of coupled equations is used to compute the total amount of each nuclide member as a function of time. The coupled equations and their required initial conditions can be mathematically expressed as follows [24,48]: where M i (t) is the amount of nuclide i per unit cross-sectional area perpendicular to the groundwater flow direction at time t [Bq/m 2 ]; γ i is the proportionality constant for nuclide i that quantifies the release rate from the repository, and M 0 i is the amount of nuclide i at the initial time. In this study, the release of the radionuclide is considered as an inflow source boundary condition. Based on the continuity of the radioactivity flux for each nuclide across the source boundary, a third-type inflow boundary condition can be specified as where ϕ is the effective porosity. Boundary conditions for all species are set at infinity as follows: The other boundary conditions required for obtaining unique solutions to Equations (1) and (2) are where The solution strategy for solving Equations (11)- (20) involves conversion of the coupled partial differential equations (PDEs) into a set of simultaneous ordinary differential equations (ODEs) via a series of integral transform, after which the solutions can be easily obtained. The solution for radioactivity concentration of each individual nuclide in the transformed domain can then be sequentially obtained by applying the general approach for solving a set of second-order ODE with constant coefficients. Taking the Laplace transforms on Equations (11), (12), (14), and (15) with respect to T and substituting the initial conditions as defined in Equations (13) and (16) for C i (X, Y, Z, s) and M i (T), respectively, yields where where s represents the Laplace transform parameter. The transformed boundary conditions from Equations (17)-(20) turn into Solving for M i (s) in sequence using Equations (23) and (24), one has Considering the second-order derivatives of the second and third terms of Equations (21) and (22), we apply the finite Fourier cosine transform twice with respect to Y and Z on Equations (21) and (22) [48]. This leads to The double finite Fourier cosine transform of C i (X, Y, Z, s) is mathematically represented as where m and n are two finite Fourier cosine transform parameters for Y and Z, respectively, and Pe z . Equations (32) and (33) are solved subject to the following boundary conditions that are the double finite Fourier cosine transform of Equations (27) and (28) The standard mathematical approach for solving a second-order ODE with constant coefficient is used to obtain the solution C i (X, m, n, s) to Equations (32) and (33). Details of the derivation are omitted herein. The analytical solutions for an arbitrary member of a radionuclide decay chain can be generalized in a concise format as follows: where Concise expressions for the target radionuclide described in Equations (37)-(40) are helpful to develop a Fortran computer program code for implementing the numerical evaluation.
The solutions in the original domain are obtained by using the double finite Fourier cosine inverse transforms and the Laplace inverse transform. This results in The analytical solution in the original domain can be obtained by performing the Laplace inverse transform on Equation (41). Although it seems to be difficult to perform analytical Laplace inversion, it can be transformed into the time domain by use of a numerical Laplace inverse algorithm such as the one developed by de Hoog et al. [49]. Several studies have demonstrated that the de Hoog et al. [49] lgorithm can give an accurate and temporally continuous evaluation of the solution [9,14,50].

Verification of the Developed Analytical Solution
A computer program code was written with FORTRAN programming language to evaluate Equation (41). Prior to applying the developed analytical model to simulate the three-dimensional radioactivity transport of a radionuclide decay chain, we test the correctness of the derived semi-analytical model and the robustness of its corresponding computer program code by comparing them with an existing analytical model developed by Chen et al. [14] for simulating two-dimensional multispecies reactive transport in a finite-length groundwater system. The test example was selected from an illustrative example of an important radionuclide decay chain first described by Higashi and Pigford [51], to demonstrate the applicability of their developed multispecies analytical transport models. The illustrative example used by Higashi and Pigford [51] considers the hydrogeological transport of a four-member radionuclide chain of 238 Pu → 234 U → 230 Th → 226 Ra released from a high-level waste (HLW) repository. Table 1 summarizes the transport parameters related to our illustrative example. A homogeneous saturated aquifer of size L × W × H = 250 m × 100 m × 100 m, as shown in Figure 2, is considered. The values of retardation factors used in Table 1 were adopted from Higashi and Pigford [51] and van Genuchten [26]. The rectangular patch boundary source at the inflow boundary is arranged on 40 m ≤ y ≤ 60 m and 0 m ≤ z ≤ 100 m segments of the inflow source boundary x = 0 m. Figures 3 and 4 depict a spatial concentration profile along the groundwater flow direction (y = 50 m) and two spatial concentration profiles perpendicular to the groundwater flow direction (x = 25 m and x = 50 m) at t = 1000 years and 10,000 years, respectively, as obtained from our derived analytical model and the analytical model developed by Chen et al. [44]. Chen et al. [44] derived the analytical solutions to the two-dimensional advection-dispersion equations coupled with sequential first-order decay reactions in a finite-length groundwater system. Numerical evaluation of their solution was always time-consuming and seems computationally inefficient due to the solutions involving summations of infinite series expansion. To reduce the impact of exit boundary on the radioactivity concentration of each individual radionuclide, the spatial concentration profiles perpendicular to the groundwater flow direction were selected far away from the finite-domain exit boundary. It can be seen that the spatial concentration profiles of four radionuclides along the groundwater flow direction and the two spatial concentration profiles perpendicular to groundwater flow direction obtained from our derived analytical model and Chens et al.'s [44] solution coincide closely to each other for a wide range of radioactivity concentrations. Excellent agreement between the two solutions clearly demonstrates the correctness of the derived analytical model and the robustness of the computer code. Subsequently, we have compared the computation times required for numerical evaluation of our solutions and Chens et al.'s [44] solutions. The computation time required for our analytical model was approximately one tenth of that required for Chen et al.'s [44] model, thus demonstrating the computation efficiency of our model.

Application of the Analytical Model
After the verification exercise, the derived model was then employed to simulate the threedimensional reactive transport of the same four-member radionuclide decay chain. The application example is similar to the verification example in that it simulates in two-dimensions the plume migration of a radionuclide decay chain, except for a slight modification of the inflow source. A

Application of the Analytical Model
After the verification exercise, the derived model was then employed to simulate the three-dimensional reactive transport of the same four-member radionuclide decay chain. The application example is similar to the verification example in that it simulates in two-dimensions the plume migration of a radionuclide decay chain, except for a slight modification of the inflow source. A schematic representation of the simulation scenario is shown in Figure 5. The simulation scenario considers a homogenous saturated groundwater system with a squared patch source at the inflow boundary. The geometry of the semifinite domain is in L × W × H = ∞ m × 100 m × 100 m dimensions. The square patch area is centered at (y, z) = (50 m, 50 m) and has a width and height of 20 m and 20 m, respectively. The modeling parameters are the same as those used in the verification exercise. The application example was designed to demonstrate that the developed model can be used to compute the long-term three-dimensional transport behavior on both time and space scales. Figure 6 illustrates the spatial concentrations of the four nuclides in the radionuclide decay chain, at different depths, for a time period of t = 1000 years. It is clearly observed that the migration distances of the four nuclide plumes are quite different. The 226 Ra plume migrates a large distance because of its lower retardation factor that reflects its weaker sorption capacity as compared with the retardation factors of the other three nuclides. The fact that the difference in the behavior of the migrating plumes is significantly affected by the retardation factors of the individual nuclides shows the importance of accounting for different retardation factors for different nuclides in the modeling of the radionuclide decay chain. of 20 m and 20 m, respectively. The modeling parameters are the same as those used in the verification exercise. The application example was designed to demonstrate that the developed model can be used to compute the long-term three-dimensional transport behavior on both time and space scales. Figure 6 illustrates the spatial concentrations of the four nuclides in the radionuclide decay chain, at different depths, for a time period of t = 1000 years. It is clearly observed that the migration distances of the four nuclide plumes are quite different. The 226 Ra plume migrates a large distance because of its lower retardation factor that reflects its weaker sorption capacity as compared with the retardation factors of the other three nuclides. The fact that the difference in the behavior of the migrating plumes is significantly affected by the retardation factors of the individual nuclides shows the importance of accounting for different retardation factors for different nuclides in the modeling of the radionuclide decay chain.  The concentrations of these radionuclides as a function of time at different distances (x = 100, 250 and 500 m) from the inflow source boundary are depicted in Figure 7. The concentration of 238 Pu at x = 100 m increases in the early time period by start of decaying after about 1800 years and ultimately reaches a negligible value after about 5200 years, whereas the concentrations of 234 U, 230Th , and 226 Ra maintain the increasing trends up to 10,000 years. It is seen that 226 Ra has the highest radioactivity concentration before 3000 years, while the concentration of 234 U is greater than either 230 Th or 226 Ra up to 3000 years. The concentration at x = 250 m of 238 Pu is exceedingly low and is not shown in Figure  7b. The concentrations of the other three nuclides increase up to 10,000 years and follow the order 226 Ra > 234 U > 230 Th. The concentration at x = 500 m of 226 Ra increases up to 10,000 years, and 234 U and 230 Th reduce to 10 -25 Bq/m 3 after 3000 and 5400 years.  Safety assessment of a radioactive waste repository requires an estimate of effective dose from multiple pathways, including ingestion of contaminated agricultural products grown with irrigated contaminated water and fishes, shell fishes, mollusk, and so on, living in contaminated surface water, originated from contaminated ground water. The model developed in this study only provides predictions of radionuclide concentrations in groundwater. It should be noted that consideration of radiological exposure from ingestion of drinking water is just a part of overall dose assessment. However, the greater intake comes from drinking contaminated water in many cases. The time histories of concentrations at different distances from the inflow source boundary, shown in Figure 7 are subsequently used to assess the committed effective dose rates consumed by members of the general public through ingestion of drinking water. The committed effective dose rate is equal to the product of concentration, drinking water ingestion rate, and the ingestion dose coefficient. Thus, the committed effective dose per person from a given radionuclide decay chain through groundwater can be calculated by (43) where IR is the rate of intake [m 3   Safety assessment of a radioactive waste repository requires an estimate of effective dose from multiple pathways, including ingestion of contaminated agricultural products grown with irrigated contaminated water and fishes, shell fishes, mollusk, and so on, living in contaminated surface water, originated from contaminated ground water. The model developed in this study only provides predictions of radionuclide concentrations in groundwater. It should be noted that consideration of radiological exposure from ingestion of drinking water is just a part of overall dose assessment. However, the greater intake comes from drinking contaminated water in many cases. The time histories of concentrations at different distances from the inflow source boundary, shown in Figure 7 are subsequently used to assess the committed effective dose rates consumed by members of the general public through ingestion of drinking water. The committed effective dose rate is equal to the product of concentration, drinking water ingestion rate, and the ingestion dose coefficient. Thus, the committed effective dose per person from a given radionuclide decay chain through groundwater can be calculated by where IR is the rate of intake [m 3 /day], C i is the radioactivity concentration in groundwater of nuclide i [Bq/m 3 ], and DF i is the ingestion dose coefficient of nuclide i for the adult age group [Sv/Bq]. The ingestion dose coefficients for the adult age group for individual nuclides in the radionuclide decay chain are adopted from ICRP [52] and included in Table 1. The corresponding time history of the total dose and the doses due to individual nuclides at different distances from the origin are given in Figure 8. The magnitude of the doses for the different nuclides follows the same sequence of concentration magnitude as shown in Figure 7. For the distance of x = 100 m, 226 Ra contributes almost the total dose during the time period of 3000 years, whereas most of the total dose is contributed by 234 U after 3000 years. At a distance of x = 250 m and 500 m, most of the total dose is contributed by 226 Ra throughout the entire time period. Comparison with the applicable guidelines show that the total annual doses at x = 100 m, 250 m, and x = 500 m are all above the WHO guideline of 0.1 mSv/year for the drinking water pathway. The generalized analytical solutions can quickly and accurately predict the three-dimensional radionuclide plume migration and assess the radiological impact posed by radionuclides in the environment as a result of leakage from a nuclear waste repository or accidental discharge from a nuclear facility. The ingestion dose coefficients for the adult age group for individual nuclides in the radionuclide decay chain are adopted from ICRP [52] and included in Table 1. The corresponding time history of the total dose and the doses due to individual nuclides at different distances from the origin are given in Figure 8. The magnitude of the doses for the different nuclides follows the same sequence of concentration magnitude as shown in Figure 7. For the distance of x = 100 m, 226 Ra contributes almost the total dose during the time period of 3000 years, whereas most of the total dose is contributed by 234 U after 3000 years. At a distance of x = 250 m and 500 m, most of the total dose is contributed by 226 Ra throughout the entire time period. Comparison with the applicable guidelines show that the total annual doses at x = 100 m, 250 m, and x = 500 m are all above the WHO guideline of 0.1 mSv/year for the drinking water pathway. The generalized analytical solutions can quickly and accurately predict the three-dimensional radionuclide plume migration and assess the radiological impact posed by radionuclides in the environment as a result of leakage from a nuclear waste repository or accidental discharge from a nuclear facility.

Conclusions
A novel semi-analytical model with parsimonious mathematical expression for modeling threedimensional radioactivity transport of radionuclide decay chain is introduced in this study. The presented model can simultaneously simulate the radioactivity concentrations of each member of a radionuclide decay chain. The model verification is confirmed with close coincidences between the spatial concentration profiles obtained from the derived semi-analytical model and an existing analytical model found in the literature. An illustrative example is employed to demonstrate the applicability of our developed semi-analytical model for simulating the three-dimensional plume migration and assessing the radiological dose of a radionuclide decay chain. Despite that there are a number of time-marching numerical models that can accommodate multiple radionuclide transport of a decay chain, our analytical model, which can synchronously evaluate the three-dimensional radioactivity concentrations of all members of a radionuclide decay chain at any desired time, without the step-by-step computation, should be especially useful as a rapid and effective tool for

Conclusions
A novel semi-analytical model with parsimonious mathematical expression for modeling three-dimensional radioactivity transport of radionuclide decay chain is introduced in this study. The presented model can simultaneously simulate the radioactivity concentrations of each member of a radionuclide decay chain. The model verification is confirmed with close coincidences between the spatial concentration profiles obtained from the derived semi-analytical model and an existing analytical model found in the literature. An illustrative example is employed to demonstrate the applicability of our developed semi-analytical model for simulating the three-dimensional plume migration and assessing the radiological dose of a radionuclide decay chain. Despite that there are a number of time-marching numerical models that can accommodate multiple radionuclide transport of a decay chain, our analytical model, which can synchronously evaluate the three-dimensional radioactivity concentrations of all members of a radionuclide decay chain at any desired time, without the step-by-step computation, should be especially useful as a rapid and effective tool for assessing long-term transport behavior as well as the radionuclide impact posed by radionuclides in the environment.