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

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.

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,

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 [L2T−1]; is Darcy velocity [LT−1]; R(x, t) is retardation factor; θ(x) is porosity of the medium; xL and xR 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 c0(x) is initial concentration, Q(t) is prescribed mass flux [ML−2T−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 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 tn, Δt is the time step size, is the position of particle at time level tn, NC is the truncated integer value of the Courant number which is Δx is grid size.

Eqn (6) means that particle leaving at tn will reach the location, xi exactly at time, tn+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 cD is the prescribed concentration. The α is known as mass transfer coefficient [LT−1], which relates between mass flux and difference in concentration, cD 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.

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)

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 tn, when the soluntion is known, is concentration at node i+1 at time level tn +1, is concentration at the node i at discrete time level tn +1, is concentration at node i–1 at discrete time level tn +1 , Δxi is the grid size to location xi, Δxi+1 is grid size to location xi+1, Δxi –1 is grid size to location xi –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.

Fig. 1.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). xL and xR 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.

 

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.

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

Fig. 2.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 m2d−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 m2d−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 m2d−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, 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.

Fig. 3.Conceptual model for application problem showing domain, boundary condition and initial condition.

Table 2.Input parameters for application problem

Fig. 4.Comparison of numerical results of application problem for proposed approach and Neuman approach after different time intervals.

Fig. 5.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.