A New Numerical Model Coupling Modified Method of Characteristics and Galerkin Finite Element Method for Simulation of Solute Transport in Groundwater Flow System

A New Numerical Model Coupling Modified Method of Characteristics and Galerkin Finite Element Method for Simulation of Solute Transport in Groundwater Flow System. ABSTRACT A new numerical model coupling the modified method of characteristics (MMOC) with the Galerkin finite element method (GFE) is proposed for the assessment of groundwater pollution in regional groundwater flow system. The MMOC-GFE solves governing solute transport equation which involves the simulation of advection and dispersion parts by MMOC and GFE, respectively. This model allows the use of large time steps (large Courant numbers) and large spatial steps (large Peclet numbers) with stable and convergent solutions. The proposed model is successfully implemented for a test case which includes comparison of the model results with reported solutions of other numerical models. It is found that MMOC-GFE model is quite adequate in simulating solute transport in heterogeneous aquifer with combined pumping and injection schemes.


INTRODUCTION
Assessment of groundwater pollution using numerical models has the challenge of getting accurate and stable numerical solutions especially in highly heterogeneous media. The proposed MMOC-GFE model in this study attempts to meet this challenge by using Eulerian-Lagrangian formulation. The newly developed model controls numerical dispersion and provides stable solutions by solving advection part and hydrodynamic dispersion part by MMOC and GFE methods respectively. It is shown that for chosen test case (Chiang et al., 1989) even with high Courant and Peclet numbers the model solutions are in close agreement with reported solutions.
An adaptive Eulerian-Lagrangian formulation (Neuman et al., 1984) solves the advective component of steep concentration front by single step reverse particle tracking of moving particles clustered around each front and the dispersive transport component is solved by Lagrangian formulation on a fixed grid. Chiang et al. (1989) discussed the model combining modified method of characteristics and mixed finite element method for solute transport simulation in groundwater flow system. The model allows arbitrary placement of injection and pumping wells within or on the boundaries of the domain together with different boundary conditions. It is found that the use of particle clusters at steep concentration fronts only caused numerical dispersion. Illangasekare et al. (1989) developed a discrete kernel approach for generating the velocity field solving the groundwater flow equation. Further the two-dimensional transient solute transport in water table aquifers is simulated by employing a method of characteristics. However, the simulation of sharp concentration fronts at injection well and divergent velocity fields at pumping well caused numerical smearing of the solutions. Garder et al. (1964) proposed a method of characteristic technique is for the solution of the solute transport equation which uses imaginary moving particles to represent solute mass and to calculate the concentration change at a given point in the domain for each time step. Goode (1990) modified the conventional USGS MOC, 1988 model by using bilinear interpolation schemes for velocity interpolation. But the heterogeneities at the cell interfaces require the use of higher order interpolation schemes to accurately simulate advective transport by method of characteristics. Zhang et al. (1993) proposed a computationally efficient Eulerian-Lagrangian method for solving advection-dispersion equation in both steady and transient velocity field. The method uses single step reverse particle tracking technique for steep concentration fronts and separate weighting factors which relate to grid Peclet and Courant numbers are used for upstream and downstream region of the advection front. In transient velocity field the model determines the weighting factors automatically based on the mass balance errors. This model is found to be more suitable for the solution of the advection dominated problems. Frolkovic (2002) discussed a flux based method of characteristics for simulation of contaminant transport in flowing groundwater. This method combines advantages of numerical discretizations by finite volume methods and by methods of characteristics. This method is further extended for complex transport problems on multidimensional unstructured computational grids. It is reported that this method is difficult to implement for general computational grids and also it suffers from choice of time steps. (Lichtner, 2002) applied the new form of the dispersion tensor involving axisymmetric media utilizing particle tracking techniques. It is demonstrated that for the case of spatially variable flow the drift term of dispersion must generally be included in the particle tracking algorithm to obtain accurate results. Banas (2004) solved non-linear convectiondiffusion equation using method of characteristics. The features of the scheme are demonstrated for onedimensional numerical experiment. Pinder and Gray (1977) discussed Galerkin finite element formulation to solve dispersion part of the solute transport equation. Kinzelbach and Frind (1986) used two-dimensional Galerkin finite element model with linear elements for solution of dispersion equation and effect of design of grid orientation on model accuracy is also investigated. Sheu and Chen (2002) solved the unsteady advection-diffusion equation using finite element model which employs a quadratic basis function to approximate the contaminant concentration. The development of a weighted residuals finite element model involves constructing a biased test function to retain the scheme stability for wide ranges of values of the physical coefficients. Igboekwe (2014) examined the movement of water solute from the surface of the earth to the aquifer and also to assess the impact of existing or proposed activity on the quality and quantity of groundwater through the use of groundwater flow and solute transport models. Finite element method is used for finding approximate solutions of partial differential equations as well as that of integral equations. The solution approach is based on either eliminating the differential equation completely or rendering the partial differential equation into an approximating system of ordinary differential equation. It is found that the order of the interpolation function has greater effect on solute transport solutions than groundwater flow solutions.
The objectives of this study are i) to develop and verify a computationally efficient numerical model coupling modified method of characteristics and Galerkin finite element method for solving advection and dispersion parts of the solute transport equation, respectively, ii) to verify the accuracy and stability of model results by comparing with reported solutions and iii) to analyze the effects of porosity, dispersivity and pumping rate on solute concentration distribution.

MATHEMATICAL FORMULATION Groundwater Flow Equation
The governing equation of two-dimensional and transient groundwater flow in heterogeneous unconfined aquifer with pumping and injection wells is given as (Illangasekare et al., 1989). Equation (1) is subject to the following initial condition which is given as where ℎ 0 is the initial head [ ] over the aquifer domain Ω [ 2 ]. Equation (1) is subject to the following Dirichlet and Neumann boundary conditions where ℎ 1 is the specified head over aquifer boundary Γ 1 , [ ], is the specified flux across boundary Γ 2 [ −1 ], [[ ℎ] ∇ℎ( , , )] is the groundwater flux due to head gradient [ −1 ] and is the normal unit vector in outward direction. The x-and y-components of groundwater flow velocity are given as: where θ is the effective porosity of the aquifer and , are velocity components [ −1 ].

Solute Transport Equation
The governing equation of solute transport in two-dimensional transient groundwater flow system with pumping/injection wells can be given as (Illangasekare et al., 1989), The hydrodynamic dispersion coefficients in the tensor form can be given as (Zhang et al., 1993).
where α and α are longitudinal and transverse dispersivities, respectively [ ] and | | is the magnitude of the groundwater flow velocity [ −1 ].

NUMERICAL METHODS
The flow chart in Figure 1 describes the various steps involved in MMOC-GFE Model.

Modified Method of Characteristics
The modified method of characteristics (MMOC) is proposed as a variant of the conventional USGS-MOC model (Konikow et al., 1978). Unlike the conventional USGS-MOC model, the MMOC method solves the advection part of solute transport equation using real fluid particles. Following steps are involved in MMOC model: Step 1. At the start of the simulation, a set of fluid particles are generated in the block centered finite difference grid. Each fluid particle has assigned the solute concentration equal to the concentration of the respective grid to which it belongs.
Step 2. The Equation 5 is rearranged in the following form using material derivative term as, Step 3. The left hand side term of the above equation is solved using modified method of characteristics which involves the tracking of particles backward in time along characteristic lines. Thus the particle location at the base of the characteristic curve represents the location of the drifting particle at the previous time level, the slope of characteristic curves representing velocity components are used to compute the advective transport as follows: where δ , and δ , are the advective displacements in x-and y-directions, respectively [L].
Step 4. After each time step, the concentration in the grid is updated by computing the concentration change due to advection alone as follows: Step 5: Using bilinear interpolation scheme (Goode, 1990), from the quadrantal location of particles, the concentration at the grid point is updated as follows:

Galerkin Finite Element Method
The Galerkin finite element method (GFE) is used to compute the change in nodal concentration due to hydrodynamic dispersion as follows: Step 1: The second-order hydrodynamic dispersion term in Equation (5) is approximated using GFE method. It uses Green function to yield the following system of algebraic linear equations (Pinder and Gray, 1977), where [A] is the global conductance matrix which is formed by assembling the elemental conductance matrices [ ] which is given as, [B] is the global storage matrix which is assembled from elemental storage matrices [ ] which is given as, {d L } is the global load vector which is assembled from elemental load vectors { } which is given as, {g} is the global boundary flux vector which is assembled from the elemental boundary flux vectors { } which is given as, Step 2: From the concentration distribution obtained using Equation (13) the change in concentration due to hydrodynamic dispersion distribution i.e. , +Δ is obtained by Equation (14).
Step 3: After every time step, the final concentration at a node is computed by adding the change in concentration due to advection and dispersion.

Numerical Stability and Accuracy
The stability of the transport model solutions is checked by the criteria of Courant and Peclet numbers which are given as (Daus et al., 1986), where , and , are Courant and Peclet numbers, respectively. The accuracy of the transport model solutions is verified by the criteria of mass balance error which is given as, where is the average mass balance error [percent], Δ is the net solute flux [ppm], Δ is the change in solute mass stored [ppm] and 0 is the initial solute mass [ppm].

Test Case
For the test case an example problem (Chiang et al., 1989) is chosen as shown in Figure 2. The purpose of this test case to study the effect of porosity, dispersivity and pumping rate on solute concentration distribution. Test case involves a two-dimensional square aquifer of size 9.14 m as shown in  Figure 3 shows the velocity field of the test case described above by solving Equation (1) using Galerkin finite element method. The velocity field is obtained for chosen hydraulic conductivity filed. From Figure 3 the stream lines show that the convergent and divergent flow field is developed at the location of injection and pumping wells respectively. It is found that at the source and sink point for same pumping and injectin rates the resultant velocity is observed to be 1.7 m/d. Figure 4 shows the comparison of solute concentration breakthrough curves obtained by the MMOC-GFE model with the reported solution (Chiang et al., 1989). It is found that both the breakthrough curves match very closely during the simulation period from 12 to 17 days. However they deviate in initial and final stages of simulation up to 10% due to numerical dispersion. The 56% injected solute reaches at observation point after 20 days.    Figure 6 shows the comparison of the concentration breakthrough curves simulated with the three values of the transverse dispersivity. The magnitude of transverse dipersivity is varied as one order less and one order higher respectively compared to its base value (0.00914) i.e. 0.0914 and 0.000914 m, respectively. It is found that the concentration breakthrough simulated with the dispersivity of 0.0914 deviates from the breakthrough curve simulated with the initially chosen dispersivity value by 7%. The concentration breakthrough simulated with the dispersivity of 0.000914 m deviates from the breakthrough curve simulated with the initially chosen value of the dispersivity by 22%. Figure 7 shows the comparison of the three concentration breakthrough curves at the pumping well which is considered as observation well. The initial pumping rate of 0.279 m 3 /d varied to its half (0.1395 m 3 /d) and double value of 0.418 m 3 /d. Due to high pumping rate the concentration goes down and due to low pumping rate the less concentration mass is taken out from the contaminated aquifer resulting into the occurrence of higher concentrations at the observation well. This physical condition is adequately simulated by numerical models. The concentration breakthrough simulated with increased pumping rate drops down to 46% of that of the