Advanced
The Modified Eulerian-Lagrangian Formulation for Cauchy Boundary Condition Under Dispersion Dominated Flow Regimes: A Novel Numerical Approach and its Implication on Radioactive Nuclide Migration or Solute Transport in the Subsurface Environment
The Modified Eulerian-Lagrangian Formulation for Cauchy Boundary Condition Under Dispersion Dominated Flow Regimes: A Novel Numerical Approach and its Implication on Radioactive Nuclide Migration or Solute Transport in the Subsurface Environment
Journal of Soil and Groundwater Environment. 2015. Apr, 20(2): 10-21
Copyright © 2015, Korean Society of Soil and Groundwater Environment
  • Received : March 19, 2015
  • Published : April 30, 2015
Download
PDF
e-PUB
PubReader
PPT
Export by style
Article
Author
Metrics
Cited by
TagCloud
About the Authors
K.V., Sruthi
Groundwater Department, Korea Institute of Geoscience and Mineral Resources, Daejeon, Republic of Korea
Heejun, Suk
Groundwater Department, Korea Institute of Geoscience and Mineral Resources, Daejeon, Republic of Korea
sxh60@kigam.re.kr
Elango, Lakshmanan
Department of Geology, Anna University, Chennai 600025, India
Byung-Gon, Chae
Groundwater Department, Korea Institute of Geoscience and Mineral Resources, Daejeon, Republic of Korea
Hyun-su, Kim
Department of Earth and Environmental Sciences, Chonbuk National University, Jeonju, Republic of Korea

Abstract
The present study introduces a novel numerical approach for solving dispersion dominated problems with Cauchy boundary condition in an Eulerian-Lagrangian scheme. The study reveals the incapability of traditional Neuman approach to address the dispersion dominated problems with Cauchy boundary condition, even though it can produce reliable solution in the advection dominated regime. Also, the proposed numerical approach is applied to a real field problem of radioactive contaminant migration from radioactive waste repository which is a major current waste management issue. The performance of the proposed numerical approach is evaluated by comparing the results with numerical solutions of traditional FDM (Finite Difference Method), Neuman approach, and the analytical solution. The results show that the proposed numerical approach yields better and reliable solution for dispersion dominated regime, specifically for Peclet Numbers of less than 0.1. The proposed numerical approach is validated by applying to a real field problem of radioactive contaminant migration from radioactive waste repository of varying Peclet Number from 0.003 to 34.5. The numerical results of Neuman approach overestimates the concentration value with an order of 100 than the proposed approach during the assessment of radioactive contaminant transport from nuclear waste repository. The overestimation of concentration value could be due to the assumption that dispersion is negligible. Also our application problem confirms the existence of real field situation with advection dominated condition and dispersion dominated condition simultaneously as well as the significance or advantage of the proposed approach in the real field problem.
Keywords
1. Introduction
Ever increasing environmental concerns have spawned additional research interests in solute transport through porous media because of the importance in the areas of underground natural resource recovery, waste storage, soil physics, and environmental remediation ( Rutquist, 2012 ; Zhao and Valiappan, 1994 ). Currently, subsurface disposal of radioactive waste is one of the preferred options because it provides effective long-term isolation of waste through two main barriers: engineering and natural barriers ( National Research Council, 2005 ). But the environmental impact will be significant, if the contaminants from the subsurface repository move through the engineering barrier and enter the geological medium (the natural barrier). Groundwater flow will be the main pathway for contaminant migration to the biosphere ( Bradbury and Green, 1985 ; Patera et al., 1990 ; Eddebbarh et al., 2003 ). Extend of advection or dispersion processes in the groundwater favors the radionuclide transport through subsurface environment ( Freeze and Cherry, 1979 ). In many circumstances, complex geochemical reactions, such as sorption, biodegradation, oxidation/reduction precipitation/dissolution and decay play major role in movement of contaminants through the groundwater system and mixing with the ambient water ( Barry, 1992 ). Therefore, the systematic hydrogeological characterization and radioactive nuclide migration behavior are essential. Analysis using numerical models on a series of hypothetical cases that incorporate hydrogeological features, similar to actual radioactive repository sites can determine the factors of greatest importance in controlling the extent of contaminant migration ( Chen et al., 1999 ; Yaouti et al., 2008 ; Radu et al., 2011 ). The widely known computer codes are HST3D, FEFLOW, FEMWATER, FEMWASTE, MODFLOW, MODPATH, MT3D, TOUGH 2, and FEHM which are applicable to modeling radionuclide transport in groundwater and improving the understanding of governing process in groundwater systems ( Kipp, 1987 ; Pruess, 1991 ; Diersch, 1997 ; Birdsell et al, 2000 ; Qian et al., 2001 ; Zuo et al., 2009 ; Bergvall et al., 2011 ; Yi et al., 2012 ; Jacques et al., 2008 ). This type of numerical models and analysis will prove useful in identifying hydrogeologic parameters necessary for characterizing actual field sites, for safety assessment programmes and for selection of suitable remediation technology ( Smith et al., 1997 ; Birdsell et al., 2000 ). Therefore it is no doubt that the development of an efficient and accurate numerical approach is imperative for the purpose of appropriate management and control of radioactive waste repository.
In the context of radioactive substance or solute transport in groundwater, advection-diffusion equation has been solved analytically and numerically using appropriate initial and boundary conditions. The solute-transport equation is difficult to solve numerically because its mathematical character can vary from parabolic to hyperbolic depending on the relative strength of advection to dispersion ( Konikow et al., 2007 ). Most of the numerical methods solving advection-dispersion equation apply Eulerian, Lagrangian or mixed Lagrangian-Eulerian approaches ( Neuman, 1984 ). The Eulerian method solves the transport equation on a fixed grid by FDM or Finite Element Method (FEM), while Lagrangian method solves on either deforming grid or a fixed grid in deforming co-ordinates ( Varoglu and Finn, 1980 ). Previous studies reported that FDM or FEM is appropriate only for dispersion dominated conditions at low Peclet Numbers whereas under high Peclet Number or advection dominated situations is effectively handled by Lagrangian methods ( Cooley, 1971 ; Freeze, 1971a ; Freeze, 1971b ; Lynch and Neil, 1980 ; Neil, 1981; Neuman, 1981 ; Neuman, 1984 ; Zienkiewicz et al., 2013 ; Koch and Nowak, 2014 ). Although, the Lagrangian methods suffers from excessive deformation of grid system for long term simulation and is prone to numerical instability ( Yeh and Chou, 1981 ). Consequently, by combining the simple fixed eulerian grid with the computational power of Lagrangian method, mixed Lagrangian-Eulerian approach is evolved ( Neuman, 1981 ; Douglass and Russell, 1982 ; Neuman and Sorek, 1982 ; Neuman, 1993 ).
However, above numerical approach of mixed Lagrangian-Eulerian method ( Sorek and Braester, 1988 ; Sorek, 1988 ) is based on the assumption that the dispersion differential operator term of an advective dependent variable should be much less than advection/dispersion partial differential operator term of dispersive dependent variable ( Sorek and Braester, 1988 ; Sorek, 1988 ). This assumption is effective only at advection dominated regimes with high Peclet Numbers. Even though this assumption has been applied for the decomposition of advection and diffusion equation into two decoupled partial differential equations, the previous studies didn’t perform or show any investigation on numerical error especially when the region is dispersion dominated regime under Cauchy boundary condition. Hence, it can be expected that the numerical error may propagate and may be accumulated on the Cauchy boundary condition in a dispersion dominated regime ( Suk, 2015 ).
The aim of the present study was thus to investigate the effect of Cauchy boundary condition under low Peclet Numbers on solution accuracy of existing Lagrangian-Eulerian approaches. Here, we propose a novel numerical approach to handle the low range of Peclet Number problems when Cauchy boundary condition is assigned. Finally the importance of proposed numerical scheme is illustrated by applying to a problem considering the radionuclide transport through nuclear waste repository which is a current waste management issue.
2. Theory
The partial differential equation describing one dimensional advection-dispersion under heterogeneous condition is,
PPT Slide
Lager Image
where c is concentration of species [ ML −3 ]; x is spatial coordinate defined relative to characteristic length [ L ]; and t is time [ T ];
PPT Slide
Lager Image
D ( x , t ) is hydrodynamic dispersion coefficient tensor [ L 2 T −1 ];
PPT Slide
Lager Image
is Darcy velocity [ LT −1 ];
PPT Slide
Lager Image
R ( x , t ) is retardation factor; θ ( x ) is porosity of the medium; x L and x R are left and right boundary points; (0, Γ] is time interval of interest.
The initial conditions of Eq. (1) can be written as,
PPT Slide
Lager Image
Subjected to a Cauchy boundary condition along inflow boundary
PPT Slide
Lager Image
and Neuman boundary condition along outflow boundaries,
PPT Slide
Lager Image
Here c 0 ( x ) is initial concentration, Q ( t ) is prescribed mass flux [ ML −2 T −1 ] on Cauchy boundary condition and q is prescribed dispersive flux along Neuman boundary. Here, proposed method will describe only the Cauchy boundary condition with prescribed mass flux at inlet because traditional Eulerian-Lagrangian method has no difficulty to deal with Dirichlet boundary condition ( Neuman, 1981 ).
- 2.1. Neuman Approach
According to the Neuman approach, the Eq. (1) is solved for c by splitting into two sets of equations. One in terms of pure advection
PPT Slide
Lager Image
called as Lagrangian concentration and other in terms of dispersion,
PPT Slide
Lager Image
called as residual concentration where,
PPT Slide
Lager Image
In order to explain procedure of the Neuman approach, the problem domain is discretized into N number of nodes. First step is to obtain the Lagrangian concentration,
PPT Slide
Lager Image
reverse particle tracking method is applied by sending a fictitious particle from each node, i backward to the point
PPT Slide
Lager Image
Where, n denotes the nth discrete time level t n , Δ t is the time step size,
PPT Slide
Lager Image
is the position of particle at time level t n , N C is the truncated integer value of the Courant number which is
PPT Slide
Lager Image
Δ x is grid size.
Eqn (6) means that particle leaving
PPT Slide
Lager Image
at t n will reach the location, x i exactly at time, t n+1 and thus becomes
PPT Slide
Lager Image
For simplicity, we assume that Courant number is less than one. In Neuman method, initial concentration of
PPT Slide
Lager Image
is set to that of c as followings.
PPT Slide
Lager Image
Thus according to Eqs (5) and (8), initial condition of
PPT Slide
Lager Image
is as follows,
PPT Slide
Lager Image
Also, in Neuman method, Lagrangian concentration at Cauchy boundary was calculated as following.
PPT Slide
Lager Image
where c D is the prescribed concentration. The α is known as mass transfer coefficient [ LT −1 ], which relates between mass flux and difference in concentration, c D and
PPT Slide
Lager Image
. The concentration difference acts as driving force for the mass flux ( Van Genuchten and Wierenga, 1976 ; Herr et al., 1989 ; Cho et al., 1990 ; Imhoff et al., 1994 ). Also α controls the type of boundary condition ( Neuman, 1981 , 1984 ; Sorek and Braester, 1988 ; Diersch, 2013 ). If α →∞, Eq. (10) represents a prescribed concentration condition; if α = 0 , it is a prescribed mass flux condition; otherwise, it is mixed condition.
Generally, when dispersive flux is dominant over advective flux at Cauchy boundary condition, the Lagrangian concentration at the boundary cannot be equal to incoming concentration. While in the advection dominated case, it is practically acceptable that Lagrangian concentration can be approximated to incoming concentration. But this approximation might induce significant errors during highly dominant dispersive flux over advective flux at the Cauchy boundary condition. The proposed algorithm identified this problem and overcome it by combining finite difference method on the boundary and traditional Eulerian-Lagrangian method on interior node. The proposed method will be described in detail later.
Since the Lagrangian concentrations at all nodes are computed using Eqs (7) and (10), residual concentrations at all nodes,
PPT Slide
Lager Image
are computed as follows. Here the hydrodynamic derivative of c ,
PPT Slide
Lager Image
can be defined as,
PPT Slide
Lager Image
According to eqn (1),
PPT Slide
Lager Image
Substituting eqn (5) in eqn (12) results,
PPT Slide
Lager Image
Expand eqn (13) as follows,
PPT Slide
Lager Image
Here, it should be noted that
PPT Slide
Lager Image
( Varoglu and Finn, 1980 ). Therefore, eqn (14) can be written as follows,
PPT Slide
Lager Image
The first term of RHS indicates the dispersion dependent residual concentration and second term shows the dispersion dependent Lagrangian concentration. It should be noted that Lagrangian concentration,
PPT Slide
Lager Image
is known function since it already solved with Eqs (7) and (10). Therefore, the second term in RHS of Eq. (15) is known. Now it is needed to solve only the residual concentration part,
PPT Slide
Lager Image
. In order to solve Eq. (15), initial condition of Eq. (9) are used. The boundary condition for
PPT Slide
Lager Image
can be obtained by using Eqs (3), (4), (5), and (10) as follows.
Substituting eqn (5) in eqn (3) yields,
PPT Slide
Lager Image
Expand the eqn (16)
PPT Slide
Lager Image
When α = 0 in eqn (10) of Neuman method, it represents the Cauchy Boundary condition which yields the following eqn,
PPT Slide
Lager Image
Substitution of eqn (18) in eqn (17) results in the following equation for Cauchy boundary condition of residual concentration determination,
PPT Slide
Lager Image
In outflow boundary where particles are leaving the flow field, it has no effect on
PPT Slide
Lager Image
and are irrelevant. Also in the outflow field, it is convenient to set α →∞ ( Neuman, 1981 ), which results residual concentration at outflow boundary as below,
PPT Slide
Lager Image
This choice is not necessary and is always possible due to the insensitivity of
PPT Slide
Lager Image
to conditions prevailing at outflow boundaries.
It should be noted that Lagrangian concentration,
PPT Slide
Lager Image
is known function since it already solved with Eqs (7) and (10). Since we have governing Eq. (15), initial condition of Eq. (9), and boundary condition of Eq. (19), (20), residual concentration,
PPT Slide
Lager Image
can be solved using the finite element method or finite difference method. The present study employed central finite difference approximation for spatial terms and fully implicit backward finite difference approximation for temporal terms.
The Neuman approach needs
PPT Slide
Lager Image
at Cauchy boundary to obtain concentration at all nodes. The
PPT Slide
Lager Image
at Cauchy inflow boundary is obtained by assumption that the effect of dispersive flux is negligible. Therefore,
PPT Slide
Lager Image
at Cauchy boundary becomes exactly equal to the incoming concentration. While in the practical situation, the influence of dispersive flux on mass flux at Cauchy boundary may or may not be negligible according to the condition. Therefore, neglecting the dispersive flux term in the Cauchy boundary can induce error; these errors can be propagated to the calculation of concentration at interior nodes also. Especially during dispersion dominated problems, the induced error could be severe. To avoid the inheritance of error due to the above assumption, the proposed algorithm applies a novel numerical concept where Lagrangian concentration at Cauchy boundary is not needed to obtain the concentration at interior nodes. The proposed method is explicitly described in the following section.
- 2.2. Proposed Approach
In the proposed approach, solute transport problem with Cauchy boundary condition can be solved using Lagrangian form of equation as following ( Neuman, 1981 )
PPT Slide
Lager Image
A fully implicit central finite difference approximation for spatial terms can be written in the proposed approach as,
PPT Slide
Lager Image
Where, n denotes the discrete time level t n , when the soluntion is known,
PPT Slide
Lager Image
is concentration at node i +1 at time level t n +1 ,
PPT Slide
Lager Image
is concentration at the node i at discrete time level t n +1 ,
PPT Slide
Lager Image
is concentration at node i –1 at discrete time level t n +1 , Δ x i is the grid size to location x i , Δ x i+1 is grid size to location x i+1 , Δ x i –1 is grid size to location x i –1 ,
PPT Slide
Lager Image
is the internodal dispersion coefficient between nodes i and i +1,
PPT Slide
Lager Image
is the internodal dispersion coefficient between nodes i and i –1 .
Temporal variations in concentration can be approximated using fully implicit backward Euler time approximation for the first term of left hand side of Eq. (21)
PPT Slide
Lager Image
Where Δ t is time step size and
PPT Slide
Lager Image
is concentration at node i . The finite difference expressions for spatial and temporal derivatives using Eqs (22) and (23), are rearranged in agreement with Eq. (21), by collecting all the unknowns on the LHS and all the known on the RHS. The resulting equation is as follows,
PPT Slide
Lager Image
Where,
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
Here, the Lagrangian concentration is calculated only at interior nodes using the reverse particle tracking method as that of the Neuman method. Internodal dispersion coefficient can be obtained by arithmetic average between dispersion coefficients of adjacent nodes.
PPT Slide
Lager Image
Eq. (24) has N unknowns but linear algebraic equations of
PPT Slide
Lager Image
it means that
PPT Slide
Lager Image
equations are needed to solve concentrations at the time level of n + 1. Since we assume that courant number is less than 1, one additional equation should be formulated. This equation can be made by applying mass balance principle at the Cauchy boundary cell ( Fig. 1 ). According to the mass balance principle, the incoming mass flux at Cauchy boundary minus outward flux of mass leaving at the right of boundary element will be equal to mass change rate within the boundary cell. It can be written mathematically as.
PPT Slide
Lager Image
Schematic representation of proposed approach for Cauchy boundary condition.
PPT Slide
Lager Image
Where A is area of cell facing to the flow
PPT Slide
Lager Image
( Fig. 1 ), V is the volume of boundary cell
PPT Slide
Lager Image
( Fig. 1 ), CB indicates a boundary cell adjacent to Cauchy boundary ( Fig. 1 ). xL and xR are left and right points of boundary cell ( Fig. 1 ). Substituting Eq. (3) in Eq. (30) results,
PPT Slide
Lager Image
In order to obtain the concentration at inflow boundary, Eq. (31) is discretized using fully implicit central finite difference approximation for spatial terms and fully implicit backward Euler approximation for temporal terms as following
PPT Slide
Lager Image
Rearranging the Eq. (32) for boundary node by collecting all the unknowns on the LHS and all the known on the RHS as,
PPT Slide
Lager Image
Where,
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
Assembling Eqs (24) and (33) into global matrix, it becomes simultaneous N linear algebraic equation set with N unknowns. Therefore, it can be solved for N concentrations by applying iterative solvers such as Newton’s method or Picard’s method.
It should be pointed out that the proposed approach combines two schemes such as finite difference method and traditional Eulerian-Lagrangian method. Finite difference method is applied only to boundary cell using mass balance principle which directly includes Cauchy boundary condition without any assumption or approximation. While, the traditional Eulerian-Lagrangian method is applied to all interior cells except boundary cell. Since the proposed method doesn’t compute Lagrangian concentration or advection concentration
PPT Slide
Lager Image
at inlet Cauchy boundary condition, therefore it doesn’t involve any significant errors at inlet Cauchy boundary with dispersion dominated regime.
3. Results and Discussion
Based on the novel numerical concept, a one-dimensional, mixed finite difference method and Eulerian-Lagrangian method has been developed for dispersion dominated problems under Cauchy boundary condition and the results are compared with the analytical solution ( Van Genuchten and Alves, 1982 ). The illustrated examples show the performance of novel numerical approach to handle dispersion dominated problems with inlet Cauchy boundary condition as well as inability of Neuman approach to deal with these type of problem. An application oriented problem is explained to prove the significance of proposed approach in the real field situation.
- 3.1. Comparison of the Proposed Approach with Analytical Solution, FDM, and Neuman Approach
The proposed numerical approach is implemented to compare the numerical results with analytical solution, FDM and Neuman approach ( Neuman, 1981 ). The initial and boundary condition for different approaches are shown in Table 1 and the input parameters for different problems are given in Table 1 .
Input parameters for case 1
PPT Slide
Lager Image
Input parameters for case 1
Here, the study illustrates four problems with different Peclet Numbers ranging from 0.01 to 2. The comparison between results of different numerical methods is depicted in Fig. 2 . The case 1 is extremely dispersion dominated condition with a Peclet Number of 0.01, with Cauchy boundary condition. As shown in the Fig. 2a , it is obvious that the Neuman approach severely deviates from the analytical solution, while the proposed numerical approach and FDM results show a good agreement with the analytical solution. It should be noted that FDM performs well for low Peclet Number, which has been already reported ( Neuman, 1981 , 1984 ).
PPT Slide
Lager Image
Comparison between the simulation result of FDM method, proposed approach, Neuman approach of single time stepnode scheme with analytical solution for different Peclet Number at different time interval. (a) Pe = 0.01, (b) Pe = 0.1, (c) Pe = 1, (d) Pe = 2.
The important observation is that Neuman approach produces significant errors because Lagrangian concentration at Cauchy boundary wasn’t estimated properly and errors propagated into the interior region. This error originates from the assumption inherent in Neuman method that the dispersion differential operator term of an advective dependent variable is much less than advection/dispersion partial differential operator term of a dispersive dependent variable ( Sorek, 1988 ). On the other hand the proposed approach could solve the problem without any error. It is due to the fact that this approach does not need to calculate Lagrangian concentration at inlet Cauchy boundary as Neuman method.
The case 2 adopts same input parameters as that of case 1 except for dispersion coefficient; the dispersion coefficient is 5 m 2 d −1 and the Peclet Number is 0.1. The comparison of the numerical results at t = 0.5, 2, 5 and 10 days are exhibited in Fig. 2b . The result shows that the proposed numerical approach and FDM schemes are well matched with the analytical solution, while the numerical results of Neuman approach significantly deviates from the analytical solution. However, the magnitude of deviation from analytical solution is comparatively lesser than the case 1. Therefore, from the case 1 and 2, it can be revealed that the Neuman approach is not applicable to solve dispersion dominated problems corresponding to Peclet Numbers of 0.01 to 0.1 with Cauchy boundary condition whereas the proposed approach is suitable to deal the dispersion dominated problems with Cauchy boundary condition. The case 3 represents a problem with Peclet Number of 1. The input parameters are same as case 1 except dispersion coefficient. The dispersion coefficient is 0.5 m 2 d −1 . The numerical results are given in Fig. 2c . As expected, the numerical results of FDM, proposed approach and Neuman approach show an excellent match with analytical solution. Therefore, it reveals that when the Peclet Number increases or advection dominates, Neuman approach tends to be in well accordance with the analytical solution. Also the following cases can confirm that Neuman approach is best for advection dominated problems rather than dispersion dominated cases with Cauchy boundary condition. The case 4 deals with a problem of Peclet Number of 2 and the dispersion coefficient is 0.25 m 2 d −1 while the other input parameters are same as case 1. The numerical results of the FDM, proposed approach and Neuman approach yields very close results with analytical solution ( Fig. 2d ). Also it can be observed from the results that as Peclet Number increases, the numerical results of the Neuman approach become closer to the analytical solution. It has to be noted that the proposed method produces very close solution with analytical solution for wide range of Peclet Numbers, which represent from dispersion dominated condition to advection dominated condition.
Therefore the novel numerical approach addresses better both the dispersion dominated problem and advection dominated problem simultaneously with Cauchy boundary condition which cannot be accurately handled by Neuman approach.
- 3.2. Application Problem
In order to validate the applicability of proposed numerical code in the real field, a conceptual model is developed. In this conceptual model, radioactive waste disposal in the near surface is considered. Geological repositories for disposal of nuclear wastes generally rely on a multi-barrier system to isolate radioactive wastes from the biosphere. It typically consists of a geological barrier system (GBS), including repository host rock and its surrounding subsurface environment, and an engineering barrier system (EBS). EBS represents the man-made, engineered materials placed within a repository. The main component used for EBS is very low permeable Bentonite clay to limit any advection transfer ( Delage et al., 2010 ; Bianchi et al., 2014 ). In this real situation, if the radionuclide transport occurs or engineering barrier layer fails, initially it transports through low hydraulic conductivity medium of bentonite clay (diffusion dominated) and it enters to the GBS later onwards gradually it enters to the geological medium (advection dominated) ( Bianchi et al., 2014 ). The conceptual model for the application problem has shown in Fig. 3 . Here the EBS layer is made by bentonite clay, the GBS is composed of alluvial clay and the geological medium consist moderately weathered granite as well as highly weathered granite. The hydraulic conductivity and dispersion coefficient of each medium varies drastically from one layer to another. The seepage velocity and dispersion coefficient for EBS, geological barrier and geological medium are adopted from previous studies which have shown in Table 2 ( Delage et al., 2010 ; Yi et al., 2012 ). The peclet Number (Pe) of EBS and geological media were calculated using equation,
PPT Slide
Lager Image
by applying adopted hydrological parameters. Here the peclet Number of EBS is 0.02 and GBS is 0.003. The Peclet Number in the geological medium which is composed of moderately weathered granite and highly weathered granite varies from 0.03-34.5. It is obvious that the Peclet Number of the system varies drastically. Therefore, it is clear that this type of real situation experiences diffusion dominated condition and advection dominated condition simultaneously. The present study already revealed that Neuman approach fails to handle dispersion dominated problems. Therefore, this type of real situation can be handled only by the proposed approach. The numerical results of Neuman approach and proposed approach are compared each other for application problem ( Fig. 4 ). From the results it is clear that Neuman approach exaggerates or overestimates the concentration value in the order of 100 compared to the numerical results of proposed approach ( Fig. 4 ). Also the similar type of behavior is exhibited by Neuman approach in numerical result of dispersion dominated problem ( Fig. 2a ). The comparison of break through curve from Neuman approach and proposed approach at a distance of 1km reveals that Neuman approach overestimate the concentration of nuclides ( Fig. 5 ). According to Neuman approach, the radionuclide concentration becomes extremely high level (~2600 mg/L) at 10,000 years, while the proposed approach yields a concentration value of less than 6 ppm at 10,000 years ( Fig. 5 ). This type of overestimation in Neuman approach could be due to the assumption that dispersion is negligible, whereas the proposed approach does not impose any assumption. Therefore, it is clear that Neuman approach can lead to the misinterpretation of radioactive migration in real field situations, which can mislead the systematic management of radioactive waste repository. Therefore, this application problem confirms the existence of real situation of varying diffusion dominated condition to advection dominated condition simultaneously as well as the significance of the proposed approach to handle these real field problems especially during radioactive nuclide migration from waste repository.
PPT Slide
Lager Image
Conceptual model for application problem showing domain, boundary condition and initial condition.
Input parameters for application problem
PPT Slide
Lager Image
Input parameters for application problem
PPT Slide
Lager Image
Comparison of numerical results of application problem for proposed approach and Neuman approach after different time intervals.
PPT Slide
Lager Image
Break through curve at a location x = 1000 m for application problem.
4. Conclusions
A novel numerical approach is introduced for addressing the dispersion dominated problems with Cauchy boundary condition in the Eulerian-Lagrangian framework. The significance of the proposed numerical approach is that it can handle dispersion dominated cases with Cauchy boundary condition with extreme numerical accuracy. The novel numerical approach is developed by combining finite difference method at the boundary and traditional Eulerian-Lagrangian method at the interior. And it does not assign any assumptions and is not involved Lagrangian concentration at inlet Cauchy boundary to obtain the concentration at inlet boundary node, which eliminates the error propagation. We compared the numerical solutions of traditional FDM, Neuman approach, and the proposed numerical approach against the analytical solution under Cauchy boundary condition with various Peclet Numbers. The numerical results of proposed approach exhibited a well accordance with the analytical solution for the dispersion dominated problem while the Neuman approach significantly deviates from analytical solution. Therefore, the present study could reveal the limitation or inability of traditional Neuman approach to deal the dispersion dominated problem with Cauchy boundary condition. In addition, the proposed approach is validated against a real field problem of migration of radioactive nuclide from radioactive waste repository by adopting input parameters from previous studies. The application problem reveals the ability of the proposed approach to handle dispersion dominated and advection dominated condition simultaneously, where the Neuman approach fails to handle this situation. The Neuman approach overestimates the concentration value during these types of real field problems and it can lead to the misinterpretation of radioactive nuclide migration from waste repository and issues in the management of the repository. Therefore the proposed method is more reliable in these types of real field situations. Thus, we conclude that the developed novel numerical approach can produce very reliable results especially for dispersion dominated problems with Cauchy boundary condition. Moreover, the extension of the proposed numerical approach to multidimensional problems does not pose any conceptual difficulty.
Acknowledgements
This work was supported by the Radioactive Waste Management of the Korea Institute of Energy Technology Evaluation and Planning (KETEP) grant funded by the Korea government Ministry of Knowledge Economy (2011T100 200152).
References
Barry D.A. , Ghadiri H. , Rose C. W 1992 Modelling chemical transport in soil: natural and applied contaminants Lewis Publishers Boca Raton, Florida 105 - 144
Bergvall M. , Grip H. , Sjostrom J. , Laudon H. 2011 Modeling subsurface transport in extensive glaciofluvial and littoral sediments to remediate a municipal drinking water aquifer Hydrol. Earth Syst. Sci. 15 2229 - 2244    DOI : 10.5194/hess-15-2229-2011
Bianchi M. , Liu H.H. , Birkholzer J. 2014 Radionuclide Transport Behavior in a Generic Geological Radioactive Waste Repository Groundwater. 1 - 12
Birdsell K.H. , Wolfsberg A.V. , Hollis D. , Cherry T.A. , Bower K.M. 2000 Groundwater Flow and Radionuclide Transport Calculations for a Performance Assessment of a Low-Level Waste Site J. Contam. Hydrol. 46 99 - 129    DOI : 10.1016/S0169-7722(00)00129-7
Bradbury M.H. , Green A. 1985 Measurement of important parameters determining aqueous phase diffusion rates through crystalline rock matrices J. Hydrol. 82 39 - 55    DOI : 10.1016/0022-1694(85)90045-9
Chen D.W , Carsel R.F. , Moeti L. , Vona B. 1999 Assessment and prediction of contaminant transport and migration at a Florida superfund site Environ. Monit. Assess. 57 291 - 299    DOI : 10.1023/A:1006046009829
Cho J.H. , Jaffe P.R. 1990 The volatilization of organic compounds in unsaturated porous media during infiltration J. Contam. Hydrol. 6 387 - 410    DOI : 10.1016/0169-7722(90)90036-G
Cooley R.L. 1971 A finite difference method for unsteady flow in variably saturated porous media: application to a single pumping well Water. Resour. Res. 7 1607 - 1625    DOI : 10.1029/WR007i006p01607
Delage P. , Cui Y.J. , Tang A.M. 2010 Clays in radioactive waste disposal J. Rock Mech. Geotech. Eng. 2 111 - 123
Diersch H.J.G. 1997 Interactive, graphics-based finite element simulation system FEFLOW for modeling groundwater flow, contaminant mass and heat transport processes, FEFLOW User’s Manual Version 4.7 WASY Institute for Water Resources Planning and Systems Research Berlin
Diersch H.J.G. 2013 Finite Element Modeling of Flow, Mass and Heat Transport in Porous and Fractured Media Springer New York
Douglass J. , Russell T.F. 1982 Numerical methods for convection dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures SIAM J. Numer. Anal. 19 871 - 885    DOI : 10.1137/0719063
Eddebbarh A.A. , Zyvoloski G.A. , Robinson B.A. , Kwicklis E.M. , Arnold B. , Corbet T. , Kuzio S. , Faunt C. 2003 The saturated zone at Yucca Mountain: an overview of the characterization and assessment of the saturated zone as a barrier to potential radionuclide migration J. Contam. Hyd. 62-63 477 - 493    DOI : 10.1016/S0169-7722(02)00154-7
Freeze R.A. 1971 Three dimensional transient, saturated-unsaturated flow in a groundwater basin Water. Resour. Res. 7 347 - 366    DOI : 10.1029/WR007i002p00347
Freeze R.A. 1971 Influence of the unsaturated flow domain on seepage through earth dams Water. Resour. Res. 7 929 - 941    DOI : 10.1029/WR007i004p00929
Freeze R.A. , Cherry J.A. 1979 Groundwater Prentice-Hall Englewood Cliffs, New Jersey
Herr M. , Schafer G. , Spitz K. 1989 Experiments studies of mass transport in porous media with local heterogeneities J. Contam. Hydrol. 4 127 - 137    DOI : 10.1016/0169-7722(89)90017-X
Imhoff P.T. , Jaffe P.R. 1994 Effect of Liquid Distribution on Gas-water phase mass transfer in unsaturated sand during infiltration J. Contam. Hydrol. 16 359 - 380    DOI : 10.1016/0169-7722(94)90044-2
Jacques D. , Simu nek J. , Mallants D. , van Genuchten M.Th. 2008 Modelling coupled water flow, solute transport and geochemical reactions affecting heavy metal migration in a podzol soil Geoderma. 145 449 - 461    DOI : 10.1016/j.geoderma.2008.01.009
Kipp K.L. 1987 HST3D: a computer code for simulation of heat and solute transport in three dimensional groundwater flow systems. US Geological Survey, WRIR report Denver, Colorado, USA 86 - 4095
Koch J. , Nowak W.A. 2014 Method for implementing Dirichlet and third-type boundary conditions in PTRW simulations Water Resour. Res. 50 1374 - 1395    DOI : 10.1002/2013WR013796
Konikow L.F. , Reilly T.E. , Barlow P.M. , Voss C.I. , Delleur J 2007 The Handbook of Groundwater Engineering CRC Press Boca Raton Chapter 23, Groundwater modeling 54 -
Lynch D. , O’Neill K. , Wang S.Y. , Alonso C.V. , Brebbia C.A. , Gray W.G. , Pinder G.F. 1980 Finite Elements in Water Resources Oxford, Mississippi Elastic grid deformation for moving boundary problems in two space dimensions
2005 Risk and decisions about disposition of transuranic and high-level radioactive waste National Academy Washington, DC
Neuman S.P. 1981 An Eulerian-Lagrangian numerical scheme for the dispersion-convection equation using conjugate spacetime grids J. Comput. Phys. 41 270 - 294    DOI : 10.1016/0021-9991(81)90097-8
Neuman S.P. , Sorek S. , Holz K.P. 1982 Finite Elements in Water Resources Springer Newyork Eulerian-Lagrangian methods for advection-dispersion
Neuman S.P. 1984 Adaptive Eulerian-Lagrangian finite element method for advection-dispersion Int. J. Numer. Meth. Eng. 20 321 - 337    DOI : 10.1002/nme.1620200211
Neuman S.P. 1993 Eulerian-Lagrangian theory of transport in space-time nonstationary velocity fields; exact formalism by conditional moments and weak approximation Water. Resour. Res. 19 633 - 645
O’Neill K. 1981 Highly efficient, oscillation free solution of the transport equation over long times and large space Water. Resour. Res. 17 1665 - 1675    DOI : 10.1029/WR017i006p01665
Patera E.S. , Hobart D.E. , Meijer A. , Rundberg R.S. 1990 Chemical and physical processes of radionuclide migration at Yucca mountain, Nevada J. Radioanal. Nucl. Chem. 142 331 - 337    DOI : 10.1007/BF02039472
Pruess K. 1991 TOUGH2: a general-purpose numerical simulator for multiphase fluid and heat flow. Lawrence Berkeley Laboratory Report Series Berkeley, California, USA
Qian T.W. , Tang H.X. , Chen J.J. , Wang J.S. , Liu C.L. , Wu G.B. 2001 Simulation of the migration of 85Sr in Chinese loess under artificial rainfall condition Radiochim. Acta. 89 403 -
Radu F.A. , Suciu N. , Hoffmann J. , Vogel A. , Kolditz O. , Park C.H. , Attinger S. 2011 Accuracy of numerical simulations of contaminant transport in heterogeneous aquifer: a comparative study Adv. Water Resour. 34 (1) 47 - 61    DOI : 10.1016/j.advwatres.2010.09.012
Rutqvist J. 2012 The geomechanics of CO2 storage in deep sedimentary formations Geotech. Geol. Eng. 30 525 - 51    DOI : 10.1007/s10706-011-9491-0
Smith P.A. , Gautschi A. , Vomvoris S. , Zuidema P. , Mazurek M. 1997 The development of a safety assessment model of the geosphere for a repository sited in the crystalline basement of northern Switzerland J. Contam. Hydrol. 26 309 - 324    DOI : 10.1016/S0169-7722(96)00078-2
Sorek S. 1988 Eulerian-Lagrangian method for solving transport in aquifers Adv. Water Resour. 11 67 - 73    DOI : 10.1016/0309-1708(88)90039-5
Sorek S. , Braester C. 1988 Eulerian-Lagrangian Formulation of the Equations for Groundwater Denitrification Using Bacterial Activity Adv. Water Resour. 11 162 - 169    DOI : 10.1016/0309-1708(88)90029-2
Suk H. 2015 Modified mixed Lagrangian- Euleraian method based on numerical framework of MT3DMS on Cauchy boundary Groundwater (under review)
Van Genuchten M. , Wierenga P.J. 1976 Mass transfer studies in sorbing porous media. I. Analytical solutions Soil Sci. Soc. Am. J. 40 473 - 481    DOI : 10.2136/sssaj1976.03615995004000040011x
Van Genuchten M.T. , Alves W.J. 1982 Analytical solutions of the one-dimensional convective-dispersive solute transport equation, Technical Bulletin United States Department of Agriculture 1661 -
Varoglu E. , Finn W.D.L. 1980 Finite elements incorporating characteristics for one-dimensional diffusion-convection equation J. Comput. Phys. 34 371 - 389    DOI : 10.1016/0021-9991(80)90095-9
Yaouti F. , Mandour A. , Khattach D. , Kaufmann O. 2008 Modeling groundwater flow and advective contaminant transport in the Bou-Areg unconfined aquifer (NE Morocco) J. Hydro-environ. Res. 2 192 - 209    DOI : 10.1016/j.jher.2008.08.003
Yeh G.T. , Chou F.K. 1981 Moving boundary numerical surge model J. Waterw. Port Coastal Ocean Div. Am. Soc. Civ. Eng. 107 34 - 37
Yi S. , Ma H. , Zheng C. , Zhu X. , Wang H. , Li X. , Hu X. , Qin J. 2012 Assessment of site conditions for disposal of low- and intermediate-level radioactive wastes: A case study in southern China Sci. Total Environ. 414 624 - 631    DOI : 10.1016/j.scitotenv.2011.10.060
Zhao C. , Valliappan S. 1994 Numerical modeling of transient contaminant migration problems in infinite porous fractured media using finite/infinite element technique, Part I: Theory Int J .Numer. Anal. Met. 18 523 - 541    DOI : 10.1002/nag.1610180802
Zienkiewicz O.C. , Taylor R.L. , Zhu J.Z. 2013 The Finite Element Method its Basis and Fundamentals 7th ed. Butterworth-Heinemann, Oxford
Zuo R. , Teng Y. , Wang J. 2009 Sorption and retardation of strontium in fine-particle media from a VLLW disposal site J. Radioanal. Nucl. Chem. 279 893 - 899    DOI : 10.1007/s10967-008-7394-1