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

- Received : March 19, 2015
- Published : April 30, 2015

Download

PDF

e-PUB

PubReader

PPT

Export by style

Article

Metrics

Cited by

TagCloud

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.
where
c
is concentration of species [
ML
^{−3}
];
x
is spatial coordinate defined relative to characteristic length [
L
]; and
t
is time [
T
];
D
(
x
,
t
) is hydrodynamic dispersion coefficient tensor [
L
^{2}
T
^{−1}
];
is Darcy velocity [
LT
^{−1}
];
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,
Subjected to a Cauchy boundary condition along inflow boundary
and Neuman boundary condition along outflow boundaries,
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
).
c
by splitting into two sets of equations. One in terms of pure advection
called as Lagrangian concentration and other in terms of dispersion,
called as residual concentration where,
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,
reverse particle tracking method is applied by sending a fictitious particle from each node,
i
backward to the point
Where,
n
denotes the nth discrete time level
t
_{n}
, Δ
t
is the time step size,
is the position of particle at time level
t
_{n}
,
N
_{C}
is the truncated integer value of the Courant number which is
Δ
x
is grid size.
Eqn (6) means that particle leaving
at
t
_{n}
will reach the location,
x
_{i}
exactly at time,
t
_{n+1}
and thus becomes
For simplicity, we assume that Courant number is less than one. In Neuman method, initial concentration of
is set to that of c as followings.
Thus according to Eqs (5) and (8), initial condition of
is as follows,
Also, in Neuman method, Lagrangian concentration at Cauchy boundary was calculated as following.
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
. 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,
are computed as follows. Here the hydrodynamic derivative of
c
,
can be defined as,
According to eqn (1),
Substituting eqn (5) in eqn (12) results,
Expand eqn (13) as follows,
Here, it should be noted that
(
Varoglu and Finn, 1980
). Therefore, eqn (14) can be written as follows,
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,
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,
. In order to solve Eq. (15), initial condition of Eq. (9) are used. The boundary condition for
can be obtained by using Eqs (3), (4), (5), and (10) as follows.
Substituting eqn (5) in eqn (3) yields,
Expand the eqn (16)
When
α
= 0 in eqn (10) of Neuman method, it represents the Cauchy Boundary condition which yields the following eqn,
Substitution of eqn (18) in eqn (17) results in the following equation for Cauchy boundary condition of residual concentration determination,
In outflow boundary where particles are leaving the flow field, it has no effect on
and are irrelevant. Also in the outflow field, it is convenient to set
α
→∞ (
Neuman, 1981
), which results residual concentration at outflow boundary as below,
This choice is not necessary and is always possible due to the insensitivity of
to conditions prevailing at outflow boundaries.
It should be noted that Lagrangian concentration,
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,
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
at Cauchy boundary to obtain concentration at all nodes. The
at Cauchy inflow boundary is obtained by assumption that the effect of dispersive flux is negligible. Therefore,
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.
A fully implicit central finite difference approximation for spatial terms can be written in the proposed approach as,
Where,
n
denotes the discrete time level
t
_{n}
, when the soluntion is known,
is concentration at node
i
+1 at time level
t
_{n +1}
,
is concentration at the node
i
at discrete time level
t
_{n +1}
,
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}
,
is the internodal dispersion coefficient between nodes
i
and
i
+1,
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)
Where Δ
t
is time step size and
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,
Where,
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.
Eq. (24) has
N
unknowns but linear algebraic equations of
it means that
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.
Schematic representation of proposed approach for Cauchy boundary condition.
Where
A
is area of cell facing to the flow
(
Fig. 1
), V is the volume of boundary cell
(
Fig. 1
), CB indicates a boundary cell adjacent to Cauchy boundary (
Fig. 1
).
x_{L}
and
x_{R}
are left and right points of boundary cell (
Fig. 1
). Substituting Eq. (3) in Eq. (30) results,
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
Rearranging the Eq. (32) for boundary node by collecting all the unknowns on the LHS and all the known on the RHS as,
Where,
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
at inlet Cauchy boundary condition, therefore it doesn’t involve any significant errors at inlet Cauchy boundary with dispersion dominated regime.
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
).
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.
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.
Conceptual model for application problem showing domain, boundary condition and initial condition.
Input parameters for application problem
Comparison of numerical results of application problem for proposed approach and Neuman approach after different time intervals.
Break through curve at a location x = 1000 m for application problem.

Dispersion dominated Cauchy boundary condition
;
Neuman approach
;
Novel numerical approach
;
Radioactive waste repository

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

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

- 2.1. Neuman Approach

According to the Neuman approach, the Eq. (1) is solved for
PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

- 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

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

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

PPT Slide

Lager Image

- 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

PPT Slide

Lager Image

Input parameters for application problem

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

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

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**

Citing '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
'

@article{ JGSTB5_2015_v20n2_10}
,title={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}
,volume={2}
, url={http://dx.doi.org/10.7857/JSGE.2015.20.2.010}, DOI={10.7857/JSGE.2015.20.2.010}
, number= {2}
, journal={Journal of Soil and Groundwater Environment}
, publisher={Korean Society of Soil and Groundwater Environment}
, author={Sruthi, K.V.
and
Suk, Heejun
and
Lakshmanan, Elango
and
Chae, Byung-Gon
and
Kim, Hyun-su}
, year={2015}
, month={Apr}