Next Article in Journal
Environmental Profile of the Manufacturing Process of Perovskite Photovoltaics: Harmonization of Life Cycle Assessment Studies
Previous Article in Journal
Uncovering Household Carbon Footprint Drivers in an Aging, Shrinking Society
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Graduate Institute of Applied Geology, National Central University, Jhongli District, Taoyuan 32001, Taiwan
2
Department of Nursing, Fooyin University, Kaohsiung 83101, Taiwan
3
Fuel Cycle and Materials Administration, Atomic Energy Council, Taipei 23452, Taiwan
*
Authors to whom correspondence should be addressed.
Energies 2019, 12(19), 3740; https://doi.org/10.3390/en12193740
Submission received: 24 August 2019 / Revised: 23 September 2019 / Accepted: 27 September 2019 / Published: 30 September 2019
(This article belongs to the Section B: Energy and Environment)

Abstract

:
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.

Graphical Abstract

1. 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 [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23]. 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.

2. 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
D x 2 C 1 ( x , y , z , t ) x 2 + D y 2 C 1 ( x , y , z , t ) y 2 + D z 2 C 1 ( x , y , z , t ) z 2 v C 1 ( x , y , z , t ) x λ 1 R 1 C 1 ( x , y , z , t ) = R 1 C 1 ( x , y , z , t ) t
D x 2 C i ( x , y , z , t ) x 2 + D y 2 C i ( x , y , z , t ) y 2 + D z 2 C i ( x , y , z , t ) z 2 v C i ( x , y , z , t ) x λ i R i C i ( x , y , z , t ) + λ i 1 R i 1 C i 1 ( x , y , z , t ) = R i C i ( x , y , z , t ) t ,   i = 2 , , N
where C i ( x , y , z , t ) represents the radioactivity concentration of nuclide i [Bq/m3]; 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.
It is assumed that the aquifer is free of radionuclides before their release from the source area and can be stated as
C i ( x , y , z , t = 0 ) = 0 , i = 1 , , N
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]:
d M 1 ( t ) d t = λ 1 M 1 ( t ) γ 1 M 1 ( t )
d M i ( t ) d t = λ i 1 M i 1 ( t ) λ i M i ( t ) γ i M i ( t ) , i = 2 , , N
M i ( t = 0 ) = M i 0 , i = 1 , , N
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/m2]; γ i is the proportionality constant for nuclide i that quantifies the release rate from the repository, and M i 0 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
φ v C i ( x = 0 , y , z , t ) φ D x C i ( x = 0 , y , z , t ) x = { γ i M i ( t ) y 1 y y 2 ,   z 1 z z 2 0 elsewhere
where φ is the effective porosity.
Boundary conditions for all species are set at infinity as follows:
C i ( x , y , z , t ) = 0
The other boundary conditions required for obtaining unique solutions to Equations (1) and (2) are
C i ( x , y = 0 , z , t ) y = C i ( x , y = W , z , t ) y = 0
C i ( x , y , z = 0 , t ) z = C i ( x , y , z = H , t ) z = 0
Using X = x L , Y = y W , Z = z H , T = t L v = v t L , we can express Equations (1)–(10) in dimensionless form as
1 P e x 2 C 1 ( X , Y , Z , T ) X 2 + η y 2 P e y 2 C 1 ( X , Y , Z , T ) Y 2 + η z 2 P e z 2 C 1 ( X , Y , Z , T ) Z 2 C 1 ( X , Y , Z , T ) X μ 1 R 1 C 1 ( X , Y , Z , T ) = R 1 C 1 ( X , Y , Z , T ) T
1 P e x 2 C i ( X , Y , Z , T ) X 2 + η y 2 P e y 2 C i ( X , Y , Z , T ) Y 2 + η z 2 P e z 2 C i ( X , Y , Z , T ) Z 2 C i ( X , Y , Z , T ) X μ i R i C i ( X , Y , Z , T ) + μ i 1 R i 1 C i 1 ( X , Y , Z , T ) = R i C i ( X , Y , Z , T ) T
C i ( X , Y , Z , T = 0 ) = 0 , i = 1 , , N
d M 1 ( T ) d T = μ 1 M 1 ( T ) δ 1 M 1 ( T )
d M i ( T ) d T = μ i 1 M i 1 ( T ) μ i M i ( T ) δ i M i ( T ) ,   i = 2 , , N
M i ( T = 0 ) = M i 0 ,   i = 1 , , N
C i ( X = 0 , Y , Z , T ) 1 P e C i ( X = 0 , Y , Z , T ) X = { γ i v φ M i ( T ) Y 1 Y Y 2 ,   Z 1 Z Z 2 0 elsewhere
C i ( X 0 , Y , Z , T ) = 0
C i ( X , Y = 0 , Z , T ) Y = C i ( X , Y = 1 , Z , T ) Y = 0
C i ( X , Y , Z = 0 , T ) Z = C i ( X , Y , Z = 1 , T ) Z = 0
where P e x = v L D x , P e y = v L D y , P e y = v L D y , η y 2 = L 2 W 2 , η z 2 = L 2 H 2 , μ i = λ i L v and δ i = γ i L v .
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
1 P e x 2 C 1 ¯ ( X , Y , Z , s ) X 2 + η y 2 P e y 2 C 1 ¯ ( X , Y , Z , s ) Y 2 + η z 2 P e z 2 C 1 ¯ ( X , Y , Z , s ) Z 2 C 1 ¯ ( X , Y , Z , s ) X ( s R 1 + μ 1 R 1 ) C 1 ¯ ( X , Y , Z , s ) = 0
1 P e x 2 C i ¯ ( X , Y , Z , s ) X 2 + η y 2 P e y 2 C i ¯ ( X , Y , Z , s ) Y 2 + η z 2 P e z 2 C i ¯ ( X , Y , Z , s ) Z 2 C i ¯ ( X , Y , Z , s ) X ( s R i + μ i ) C i ¯ ( X , Y , Z , s ) = μ i 1 R i 1 C i 1 ¯ ( X , Y , Z , s ) ,   i = 2 , , N
M 1 ¯ ( s ) = M 1 0 s + ψ 1
M i ¯ ( s ) = μ i 1 M i 1 ¯ ( s ) + M i 0 s + ψ i ,   i = 2 , , N
where ψ i = μ i + δ i .
The Laplace transform of C i ¯ ( X , Y , Z , s ) and M i ¯ ( s ) is defined as
C i ¯ ( X , Y , Z , s ) = 0 e s T C i ( X , Y , Z , T ) d T
M i ¯ ( s ) = 0 e s T M i ( T ) d T
where s represents the Laplace transform parameter.
The transformed boundary conditions from Equations (17)–(20) turn into
C i ¯ ( X = 0 , Y , Z , s ) 1 P e x C i ¯ ( X = 0 , Y , Z , s ) X = { γ i v φ M i ¯   ( s ) Y 1 Y Y 2 ,   Z 1 Z Z 2 0 elsewhere
C i ¯ ( X , Y , Z , s ) = 0
C i ¯ ( X , Y = 0 , Z , s ) Y = C i ¯ ( X , Y = 1 , Z , s ) Y = 0
C i ¯ ( X , Y , Z = 0 , s ) Z = C i ¯ ( X , Y , Z = 1 , s ) Z = 0
Solving for M i ¯ ( s ) in sequence using Equations (23) and (24), one has
M i ¯ ( s ) = M i 0 s + λ i + k = 1 k = i 1 M i k 0 Π j =   1 j = k μ i j Π j = 0 j = k s + λ i j
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
1 P e x d 2 C 1 ¯ ¯ ( X , m , n , s ) d X 2 d C 1 ¯ ¯ ( X , m , n , s ) d X Θ 1 ( m , n , s ) C 1 ¯ ¯ ( X , m , n , s ) = 0
1 P e x d 2 C i ¯ ¯ ( X , m , n , s ) d X 2 d C i ¯ ¯ ( X , m , n , s ) d X Θ 2 ( m , n , s ) C i ¯ ¯ ( X , m , n , s ) = μ i 1 C i 1 ¯ ¯ ( X , m , n , s )
The double finite Fourier cosine transform of C i ¯ ( X , Y , Z , s ) is mathematically represented as
C i ¯ ¯ ( X , m , n , s ) = F Y Z [ C i ¯ ( X , Y , Z , s ) ] = 0 1 [ 0 1 C i ( X , Y , Z , s ) cos ( m π Y ) d Y ] cos ( n π Z ) d Z
where m and n are two finite Fourier cosine transform parameters for Y and Z , respectively, and Θ i ( m , n , s ) = s R i + μ i + η y 2 m 2 π 2 P e y + η z 2 n 2 π 2 P e 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)
C i ¯ ¯ ( X = 0 , m , n , s ) 1 P e x d C i ¯ ¯ ( X = 0 , m , n , s ) d X = γ i ψ ( m ) ζ ( n ) M i ¯ ( s ) v φ
C i ¯ ¯ ( X , m , n , s ) = 0
where ψ ( m ) = { Y 2 Y 1 m = 0 sin ( m π Y 2 ) m sin ( m π Y 1 ) m m = 1 , 2 , and ζ ( n ) = { Z 2 Z 1 n = 0 sin ( n π Y 2 ) n sin ( n π Z 1 ) n n = 1 , 2 , .
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:
C i ¯ ¯ ¯ ( X , m , n , s ) = H i e α i X + n = 1 n = i 1 P i , i n e α i n X
where
α i = P e x P e x 2 4 P e x   Θ i ( m , n , s ) 2
H i = γ i v φ ψ ( m ) ζ ( n ) M i ¯ ( s ) n = 1 n = i 1 P i , i n ( 1 α i n P e x ) 1 α i P e x
P i , i j = { μ i 1 H i 1 α i 1 2 P e x α i 1 Θ i ( m , n , s ) j = i 1 μ i 1 P i 1 , j α j 2 P e x α j Θ i ( m , n , s ) j i 1
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
C i ( X , Y , Z , T ) = L 1 [ C i ¯ ( X , Y , Z , s ) ] = L 1 [ F Y Z 1 [ C i ¯ ¯ ( X , m , n , s ) ] ] = [ L 1 [ C i ¯ ¯ ( X , m = 0 , n = 0 , s ) ] + 2 n = 1 n = L 1 [ C i ¯ ¯ ( X , m = 0 , n , s ) ] cos ( n π Z ) ] + m = 1 m = [ L 1 [ C i ¯ ¯ ( X , m , n = 0 , s ) ] + 2 n = 1 n = L 1 [ C i ¯ ¯ ( X , m , n , s ) ] cos ( n π Z ) ] cos ( m π Y )
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].

3. Results and Discussion

3.1. 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 P u 234 U 230 T h 226 R a 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 . Figure 3 and Figure 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.

3.2. 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 226Ra 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 238Pu 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 234U, 230Th, and 226Ra maintain the increasing trends up to 10,000 years. It is seen that 226Ra has the highest radioactivity concentration before 3000 years, while the concentration of 234U is greater than either 230Th or 226Ra up to 3000 years. The concentration at x = 250 m of 238Pu 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 226Ra > 234U > 230Th. The concentration at x = 500 m of 226Ra increases up to 10,000 years, and 234U and 230Th reduce to 10−25 Bq/m3 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
Committed   effective   dose = i = 1 i = N I R × C i × D F i
where I R is the rate of intake [m3/day], C i is the radioactivity concentration in groundwater of nuclide i [Bq/m3], and D F 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, 226Ra contributes almost the total dose during the time period of 3000 years, whereas most of the total dose is contributed by 234U after 3000 years. At a distance of x = 250 m and 500 m, most of the total dose is contributed by 226Ra 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.

4. 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.

Author Contributions

C.-P.L. and J.-S.C. conceived and designed the study, J.-S.C. and C.-H.C. contributed to the development of the analytical model, C.P. designed the application example and radiological dose calculation, C.-P.L., J.-S.C. and M.-H.W. interpreted the results, C.-P.L. and J.S. wrote the manuscript and made a critical revision.

Acknowledgments

The authors are grateful to Ministry of Science and Technology, Republic of China, for their financial support of this research under contract MOST 107-2623-E-008-003 -NU.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Yeh, G.T. AT123D: Analytical Transient One-, Two-, and Three- Dimensional Simulation of Waste Transport in the Aquifer System, ORNL-5602; Oak Ridge National Laboratory: Oak Ridge, TN, USA, 1981. [Google Scholar]
  2. Van Genuchten, M.T.; Alves, W.J. Analytical Solutions of the One-Dimensional Convective-Dispersive Solute Transport Equation; Technical Bulletin No. 1661; US Department of Agriculture: Washington, DC, USA, 1982. [Google Scholar]
  3. Batu, V. A generalized two-dimensional analytical solution for hydrodynamic dispersion in bounded media with the first-type boundary condition at the source. Water Resour. Res. 1989, 25, 1125–1132. [Google Scholar] [CrossRef]
  4. Batu, V. A generalized two-dimensional analytical solute transport model in bounded media for flux-type finite multiple sources. Water Resour. Res. 1993, 29, 2881–2892. [Google Scholar] [CrossRef]
  5. Batu, V. A generalized three-dimensional analytical solute transport model for multiple rectangular first-type sources. J. Hydrol. 1996, 174, 57–82. [Google Scholar] [CrossRef]
  6. Leij, F.J.; Skaggs, T.H.; Van Genuchten, M.T. Analytical solution for solute transport in three-dimensional semi-infinite porous media. Water Resour. Res. 1991, 27, 2719–2733. [Google Scholar] [CrossRef]
  7. Leij, F.J.; Toride, N.; van Genuchten, M.T. Analytical solutions for non-equilibrium solute transport in three-dimensional porous media. J. Hydrol. 1993, 151, 193–228. [Google Scholar] [CrossRef]
  8. Park, E.; Zhan, H. Analytical solutions of contaminant transport from finite one-, two, three-dimensional sources in a finite thickness aquifer. J. Contam. Hydrol. 2001, 53, 41–61. [Google Scholar] [CrossRef]
  9. Chen, J.S.; Liu, C.W.; Liao, C.M. A novel analytical power series solution for solute transport in a radially convergent flow field. J. Hydrol. 2002, 266, 120–138. [Google Scholar] [CrossRef]
  10. Chen, J.S.; Ni, C.F.; Liang, C.P.; Chiang, C.C. Analytical power series solution for contaminant transport with hyperbolic asymptotic distance-dependent dispersivity. J. Hydrol. 2008, 362, 142–149. [Google Scholar] [CrossRef]
  11. Chen, J.S.; Ni, C.F.; Liang, C.P. Analytical power series solutions to the two-dimensional advection-dispersion equation with distance-dependent dispersivities. Hydrol. Process. 2008, 22, 4670–4678. [Google Scholar] [CrossRef]
  12. Chen, J.S.; Chen, J.T.; Liu, C.W.; Liang, C.P.; Lin, C.M. Analytical solutions to two-dimensional advection–dispersion equation in cylindrical coordinates in finite domain subject to first- and third-type inlet boundary conditions. J. Hydrol. 2011, 405, 522–531. [Google Scholar] [CrossRef]
  13. Chen, J.S.; Liu, Y.H.; Liang, C.P.; Liu, C.W.; Lin, C.W. Exact analytical solutions for two-dimensional advection-dispersion equation in cylindrical coordinates subject to third-type inlet boundary condition. Adv. Water Resour. 2011, 34, 365–374. [Google Scholar] [CrossRef]
  14. Chen, J.S.; Hsu, S.Y.; Li, M.H.; Liu, C.W. Assessing the performance of a permeable reactive barrier-aquifer system using a dual-domain solute transport model. J. Hydrol. 2016, 543, 849–860. [Google Scholar] [CrossRef]
  15. Chen, J.S.; Li, L.Y.; Lai, K.H.; Liang, C.P. Analytical model for advective-dispersive transport involving flexible boundary inputs, initial distributions and zero-order productions. J. Hydrol. 2017, 554, 187–199. [Google Scholar] [CrossRef]
  16. Zhan, H.; Wen, Z.; Gao, G. An analytical solution of two-dimensional reactive solute transport in an aquifer-aquitard system. Water Resour. Res. 2009, 45, W10501. [Google Scholar] [CrossRef]
  17. Pérez Guerrero, J.S.; Skaggs, T.H.; van Genuchten, M.T. Analytical solution for multi-species contaminant transport in finite media with time-varying boundary condition. Transp. Porous Med. 2010, 85, 171–188. [Google Scholar] [CrossRef]
  18. Gao, G.; Zhan, H.; Feng, S.; Fu, B.; Ma, Y.; Huang, G. A new mobile-immobile model for reactive solute transport with scale-dependent dispersion. Water Resour. Res. 2010, 46, W08533. [Google Scholar] [CrossRef]
  19. Gao, G.; Zhan, H.; Feng, S.; Huang, G.; Fu, B. A mobile-immobile model with an asymptotic scale-dependent dispersion function. J. Hydrol. 2012, 424, 172–183. [Google Scholar] [CrossRef]
  20. Gao, G.; Fu, B.; Zhan, H.; Ma, Y. Contaminant transport in soil with depth-dependent reaction coefficients and time-dependent boundary conditions. Water Res. 2013, 47, 2507–2522. [Google Scholar] [CrossRef]
  21. Chen, J.S.; Liu, C.W. Generalized analytical solution for advection-dispersion equation in finite spatial domain with arbitrary time-dependent inlet boundary condition. Hydrol. Earth Syst. Sci. 2011, 15, 2471–2479. [Google Scholar] [CrossRef] [Green Version]
  22. Pérez Guerrero, J.S.; Pontedeiro, E.M.; van Genuchten, M.T.; Skaggs, T.H. Analytical solutions of the one-dimensional advection–dispersion solute transport equation subject to time-dependent boundary conditions. Chem. Eng. J. 2013, 221, 487–491. [Google Scholar] [CrossRef]
  23. Liang, C.P.; Hsu, S.Y.; Chen, J.S. An analytical model for solute transport in an infiltration tracer test in soil with a shallow groundwater table. J. Hydrol. 2016, 540, 129–141. [Google Scholar] [CrossRef]
  24. Nair, R.N.; Sunny, F.; Manikandan, S.T. Modelling of decay chain transport in groundwater from uranium tailings ponds. Appl. Math. Model. 2010, 34, 2300–2311. [Google Scholar] [CrossRef]
  25. Cho, C.M. Convective transport of ammonium with nitrification in soil. Can. J. Soil Sci. 1971, 51, 339–350. [Google Scholar] [CrossRef]
  26. Van Genuchten, M.T. Convective–dispersive transport of solutes involved in sequential first-order decay reactions. Comput. Geosci. 1985, 11, 129–147. [Google Scholar] [CrossRef]
  27. Lunn, M.; Lunn, R.J.; Mackay, R. Determining analytic solution of multiple species contaminant transport with sorption and decay. J. Hydrol. 1996, 180, 195–210. [Google Scholar] [CrossRef]
  28. Sun, Y.; Clement, T.P. A decomposition method for solving coupled multi-species reactive transport problems. Transp. Porous Med. 1999, 37, 327–346. [Google Scholar] [CrossRef]
  29. Sun, Y.; Peterson, J.N.; Clement, T.P. A new analytical solution for multiple species reactive transport in multiple dimensions. J. Contam. Hydrol. 1999, 35, 429–440. [Google Scholar] [CrossRef]
  30. Sun, Y.; Petersen, J.N.; Clement, T.P.; Skeen, R.S. Development of analytical solutions for multi-species transport with serial and parallel reactions. Water Resour. Res. 1999, 35, 185–190. [Google Scholar] [CrossRef]
  31. Srinivasan, V.; Clememt, T.P. Analytical solutions for sequentially coupled one-dimensional reactive transport problems-Part I: Mathematical derivations. Adv. Water Res. 2008, 31, 203–218. [Google Scholar] [CrossRef]
  32. Srinivasan, V.; Clement, T.P. Analytical solutions for sequentially coupled one-dimensional reactive transport problems-Part II: Special cases, implementation and testing. Adv. Water Res. 31, 219–232. [CrossRef]
  33. Pérez Guerrero, J.S.; Skaggs, T.H.; van Genuchten, M.T. Analytical solution for multi-species contaminant transport subject to sequential first-order decay reactions in finite media. Transp. Porous Med. 2009, 80, 373–387. [Google Scholar] [CrossRef]
  34. Chen, J.S.; Lai, K.H.; Liu, C.W.; Ni, C.F. A novel method for analytically solving multi-species advective-dispersive transport equations sequentially coupled with first-order decay reactions. J. Hydrol. 2012, 420, 191–204. [Google Scholar] [CrossRef]
  35. Chen, J.S.; Liu, C.W.; Liang, C.P.; Lai, K.H. Generalized analytical solutions to sequentially coupled multi-species advective-dispersive transport equations in a finite domain subject to an arbitrary time-dependent source boundary condition. J. Hydrol. 2012, 456–457, 101–109. [Google Scholar] [CrossRef]
  36. Bauer, P.; Attinger, S.; Kinzelbach, W. Transport of a decay chain in homogeneous porous media: Analytical solutions. J. Contam. Hydrol. 2001, 49, 217–239. [Google Scholar] [CrossRef]
  37. Montas, H.J. An analytical solution of the three-component transport equation with application to third-order transport. Water Resour. Res. 2003, 39, 1036. [Google Scholar] [CrossRef]
  38. Quezada, C.R.; Clement, T.P.; Lee, K.K. Generalized solution to multi-dimensional multi-species transport equations coupled with a first-order reaction network involving distinct retardation factors. Adv. Water Res. 2004, 27, 507–520. [Google Scholar] [CrossRef]
  39. Sudicky, E.A.; Hwang, H.T.; Illman, W.A.; Wu, Y.S. A semi-analytical solution for simulating contaminant transport subject to chain-decay reactions. J. Contam. Hydrol. 2013, 144, 20–45. [Google Scholar] [CrossRef]
  40. Kreft, A.; Zuber, A. Comment on “Flux averaged and volume averaged concentrations in continuum approaches to solute transport”. Water Resour. Res. 1986, 22, 1157–1158. [Google Scholar] [CrossRef]
  41. Parker, J.C.; van Genuchten, M.T. Flux-averaged and volume-averaged concentrations in continuum approaches to solute transport. Water Resour. Res. 1984, 20, 866–872. [Google Scholar] [CrossRef]
  42. Van Genuchten, M.T.; Parker, J.C. Boundary conditions for displacement experiments through short laboratory soil columns. Soil Sci. Soc. Am. J. 1984, 48, 703–708. [Google Scholar] [CrossRef]
  43. Parlange, J.Y.; Barry, D.A.; Starr, J.L. Comments on “Boundary conditions for displacement experiments through short laboratory soil columns”. Soil Sci. Soc. Am. J. 1985, 49, 1325. [Google Scholar] [CrossRef]
  44. Chen, J.S.; Liang, C.P.; Liu, C.W.; Li, L.Y. An analytical model for simulating two-dimensional multispecies plume migration. Hydrol. Earth Sys. Sci. 2016, 20, 733–753. [Google Scholar] [CrossRef] [Green Version]
  45. Huang, J.; Goltz, M.N. Solutions to equations incorporating the effect of rate-limited contaminant mass transfer on vadose zone remediation by soil vapor extraction. Water Resour. Res. 1999, 35, 879–883. [Google Scholar] [CrossRef]
  46. Lai, K.H.; Liu, C.W.; Liang, C.P.; Chen, J.S.; Sie, B.R. A novel method for analytically solving a radial advection-dispersion equation. J. Hydrol. 2016, 542, 532–540. [Google Scholar] [CrossRef]
  47. Chen, J.S.; Ho, Y.C.; Liang, C.P.; Wang, S.W.; Liu, C.W. Semi-analytical model for coupled multispecies advective-dispersive transport subject to rate-limited sorption. J. Hydrol. 2019, 579, 124–164. [Google Scholar] [CrossRef]
  48. Sneddon, I.H. The Use of Integral Transforms; McGraw-Hill: New York, NY, USA, 1972. [Google Scholar]
  49. De Hoog, F.R.; Knight, J.H.; Stokes, A.N. An improved method for numerical inversion of Laplace transforms. SIAM J. Sci. Stat. Comput. 1982, 3, 357–366. [Google Scholar] [CrossRef]
  50. Chen, J.S.; Chen, C.S.; Chen, C.Y. Analysis of solute transport in a divergent flow tracer test with scale-dependent dispersion. Hydrol. Process. 2007, 21, 2526–2536. [Google Scholar] [CrossRef]
  51. Higashi, K.; Pigford, T. Analytical models for migration of radionuclides in geological sorbing media. J. Nucl. Sci. Technol. 1980, 17, 700–709. [Google Scholar] [CrossRef]
  52. ICRP. Compendium of Dose Coefficients based on ICRP Publication 60; ICRP Publication: Ottawa, ON, Canada, 2012. [Google Scholar]
Figure 1. Schematic view of the three-dimensional plume migration of a radionuclide decay chain.
Figure 1. Schematic view of the three-dimensional plume migration of a radionuclide decay chain.
Energies 12 03740 g001
Figure 2. Schematic diagram of the verification scenario: radionuclide source is illustrated as square patch with dimensions of 40 m ≤ y ≤ 60 m and 0 m ≤ z ≤ 100 m at x = 0 m.
Figure 2. Schematic diagram of the verification scenario: radionuclide source is illustrated as square patch with dimensions of 40 m ≤ y ≤ 60 m and 0 m ≤ z ≤ 100 m at x = 0 m.
Energies 12 03740 g002
Figure 3. Comparison between spatial concentration profiles of a two-dimensional four-member decay chain migration at t = 1000 years obtained from the derived analytical model and Chen et al.’s [14] analytical model at t = 1000 years: (a) along the groundwater flow direction (y = 50 m); (b) perpendicular to the groundwater flow direction (x = 25 m); (c) perpendicular to the groundwater flow direction (x = 50 m).
Figure 3. Comparison between spatial concentration profiles of a two-dimensional four-member decay chain migration at t = 1000 years obtained from the derived analytical model and Chen et al.’s [14] analytical model at t = 1000 years: (a) along the groundwater flow direction (y = 50 m); (b) perpendicular to the groundwater flow direction (x = 25 m); (c) perpendicular to the groundwater flow direction (x = 50 m).
Energies 12 03740 g003aEnergies 12 03740 g003b
Figure 4. Comparison between spatial concentration profiles of a two-dimensional four-member decay chain migration at t = 10,000 years obtained from the derived analytical model and Chen et al.’s [14] analytical model (a) along the groundwater flow direction (y = 50 m); (b) perpendicular to the groundwater flow direction (x = 25 m); (c) perpendicular to the groundwater flow direction (x = 50 m).
Figure 4. Comparison between spatial concentration profiles of a two-dimensional four-member decay chain migration at t = 10,000 years obtained from the derived analytical model and Chen et al.’s [14] analytical model (a) along the groundwater flow direction (y = 50 m); (b) perpendicular to the groundwater flow direction (x = 25 m); (c) perpendicular to the groundwater flow direction (x = 50 m).
Energies 12 03740 g004aEnergies 12 03740 g004b
Figure 5. Schematic diagram of the simulation scenario: radionuclide source is illustrated as square patch with dimensions of 40 m ≤ y ≤ 60 m and 40 m ≤ z ≤ 60 m at x = 0 m.
Figure 5. Schematic diagram of the simulation scenario: radionuclide source is illustrated as square patch with dimensions of 40 m ≤ y ≤ 60 m and 40 m ≤ z ≤ 60 m at x = 0 m.
Energies 12 03740 g005
Figure 6. The concentration contours of each member of a four-member decay chain at t = 1000 years. (a) 238 Pu ; (b) 234 U ; (c) 230 Th ; (d) 226 Ra .
Figure 6. The concentration contours of each member of a four-member decay chain at t = 1000 years. (a) 238 Pu ; (b) 234 U ; (c) 230 Th ; (d) 226 Ra .
Energies 12 03740 g006aEnergies 12 03740 g006b
Figure 7. Time history of each member of a four-member decay chain at different distances from the inflow source boundary: (a) x = 100 m; (b) x = 250 m; (c) x = 500 m.
Figure 7. Time history of each member of a four-member decay chain at different distances from the inflow source boundary: (a) x = 100 m; (b) x = 250 m; (c) x = 500 m.
Energies 12 03740 g007aEnergies 12 03740 g007b
Figure 8. The total and individual effective dose of a four-member decay chain acquired through drinking groundwater pathway at different distances from the inflow source boundary: (a) x = 100 m; (b) x = 250 m; (c) x = 500 m.
Figure 8. The total and individual effective dose of a four-member decay chain acquired through drinking groundwater pathway at different distances from the inflow source boundary: (a) x = 100 m; (b) x = 250 m; (c) x = 500 m.
Energies 12 03740 g008aEnergies 12 03740 g008b
Table 1. Transport parameters used for verification involving two-dimensional reactive transport of a four-member radionuclide decay chain.
Table 1. Transport parameters used for verification involving two-dimensional reactive transport of a four-member radionuclide decay chain.
ParameterUnitValue
Domain length, L m250
Domain width, W m100
Domain height, H m100
Groundwater velocity, v m/year100
Dispersion coefficient, Dxm2/year1000
Dispersion coefficient, Dym2/year100
Dispersion coefficient, Dzm2/year100
Retardation factor, Ri
238Pu 10,000
234U 14,000
230Th 50,000
226Ra 500
Radioactive decay constant, λi
238Puyear−10.0079
234Uyear−10.0000028
230Thyear−10.0000087
226Rayear−10.00043
Initial amount of nuclide, M i 0
238PuBq/m21.5×1015
234U 0
230Th 0
226Ra 0
Proportionality constant, γi
238Puyear−10.001
234Uyear−10.001
230Thyear−10.001
226Rayear−10.001
Ingestion dose coefficient, DFi
238PuSv/Bq2.3 × 10−7
234USv/Bq4.9 × 10−8
230ThSv/Bq2.1 × 10−7
226RaSv/Bq2.8 × 10−7

Share and Cite

MDPI and ACS Style

Chen, J.-S.; Liang, C.-P.; Chang, C.-H.; Wan, M.-H. Simulating Three-Dimensional Plume Migration of a Radionuclide Decay Chain through Groundwater. Energies 2019, 12, 3740. https://doi.org/10.3390/en12193740

AMA Style

Chen J-S, Liang C-P, Chang C-H, Wan M-H. Simulating Three-Dimensional Plume Migration of a Radionuclide Decay Chain through Groundwater. Energies. 2019; 12(19):3740. https://doi.org/10.3390/en12193740

Chicago/Turabian Style

Chen, Jui-Sheng, Ching-Ping Liang, Cheng-Hung Chang, and Ming-Hsien Wan. 2019. "Simulating Three-Dimensional Plume Migration of a Radionuclide Decay Chain through Groundwater" Energies 12, no. 19: 3740. https://doi.org/10.3390/en12193740

APA Style

Chen, J. -S., Liang, C. -P., Chang, C. -H., & Wan, M. -H. (2019). Simulating Three-Dimensional Plume Migration of a Radionuclide Decay Chain through Groundwater. Energies, 12(19), 3740. https://doi.org/10.3390/en12193740

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop