A New Space-Time CE / SE Numerical Tracking of Contaminant Transport in Fractured Stratified Geologic Profiles

To date, efficient numerical simulation of contaminant transport in geologic porous media is challenged by parametric jumps resulting from stratification and the use of ideal initial/boundary conditions. Thus, to resolve some contaminant hydrology problems, this work presents the development of the Space-Time Conservation Element/Solution Element (CE/SE) scheme for advection-dispersion-reaction a-d-r transport in geologic media. The CE/SE method derives from the native form of Gauss conservation law. Therefore, it is able to effectively handle non-trivial discontinuities that may exist within the problem domain. In freshwater aquifer, stratification and other parametric jumps are examples of such discontinuity. To simulate the Nigerian experience of nitrate pollution of freshwater aquifers; the a-d-r contaminant transport model is herein solved under a time periodic nitrate fertilizer loading condition on farmlands. Results show that this approach is able to recover the wellknown field pattern of nitrate profiles under farmlands. Cyclic loading impacts more on the dispersivity of an aquifer. Hence, dispersion coefficient modulates the response of aquifers to loading frequency. However, aquifers with conductivity less than 10 m/day are almost insensitive to periodic loads. The CE/SE method is able to sense slight (i.e. order of 10) variation in hydrological parameters. Also, CE/SE computes contaminant concentration and its flux simultaneously. Thus, it facilitates a better understanding of some reported phenomena such as contaminant accumulation and localized reverse transport at the interface between fracture and matrix in geologic medium. Clearly, CE/SE is an efficient and admissible tool into the family of numerical methods available for tracking contaminant transport in porous media.


Introduction
The Space-Time Conservation Element/Solution Element (CE/SE) method is gradually emerging as a reliable tool for tackling some problems in engineering and mathematical physics.In particular, it has been used to resolve some challenging problems in aero-acoustics and related shock flows.These include the shock tube and multicomponent combustion problems as treated in [1].However, at the other end of this class of problems are some issues in environmental systems.For example, in percolation problem, process rates are slow and sometimes of geological time scale.To assess the efficacy of CE/SE method in this regime, we look at the transport of contaminants through geological profile.Although, there are many analytic, semi-analytic and numerical methods for modeling and simulation in this regime, to date, none of these methods is able to completely capture all the intricate features of the actual flow pattern as met in practice.
In approach, two basic concepts, namely; the percolation and hydro dispersive schemes are generally used to model transport of contaminants through geologic systems.However, in recent time, attempts have been made to hybridize these models.For example, [2] is an innovative work in this regard.The microscopic approach considers transport on the basis of individual pores in the medium.Therefore it models the boundaries of an aquifer as an infinite network of interconnected pores of various sizes.This represents a skeletal conglomerate of solid element having variable porosity and permeability.
On the other hand, the macroscopic framework of contaminant transport in geologic systems differentiates between coordinated and uncoordinated pore spaces.The coordinated pore spaces determine some macroscopic properties of the aquifer e.g.conductivity and permeability.However, uncoordinated pore spaces determine other properties such as storativity, adsorption and retardation.Under certain flow condition, an additional term is used to model geo-chemical reaction or physical transformation.For example, [3] is a recent formulation of hydrodispersive scheme that pays special attention to geo-chemical reaction.To simplify analysis, the macroscopic model resolves contaminant flow in a geological medium into two major zones; namely the momentum transport zone, where fluid flow is large and a region of small fluid flow dominated by other transport processes.The overall process is usually described by the parabolic or hyperbolic partial differential equations in two or higher spatial dimensions.In principle, this equation is an expression of conservation laws.However, in an ideal situation, conservation laws are best expressed as integral equations.This is more important when the domain of the problem is discontinuous.However, most existing works in this field are based on differential equations.The underlying assumption is that the domain is smooth.In natural groundwater reservoirs, this assumption does not hold.Real aquifers are known to consist of various geological formations with sharp differences in hydrological properties.
Traditionally, the Finite Difference, Finite Elements and Finite Volume methods are the well-known numerical tools that are used to solve contaminant transport problems.These tools take on the differential form of the conservation laws.Hence, they dissipate flux in the neighborhood of discontinuity.In such regions, the computational grid is usually refined to minimize flux dissipation.However, the need to satisfy the Courant, Friedrich, and Lewy (CFL) stability conditions sometimes overshoots the cost of this strategy.In the case of the FD, some robust schemes have been recently developed with enhanced flux conservation, stability and accuracy features.These schemes are built on the concepts of Total Variation Diminishing (TVD) and localized flux conservation, [4].However, TVD schemes are also computationally expensive.
On the other hand, the Finite Element Method (FEM) obtains an approximate solution of the transport equation using the trial solution approach.The individual steps involved in FEM analysis are quite simple.However, the method becomes numerically tedious especially at the point where the elementary solutions are assembled to build the global solution.According to [5], at this stage, FEM has a tendency to divorce the modeler from the physics of the problem.In addition, for the trial solution to be adequate in regions of rapid changes, higher order basis functions or very fine grids of complex geometry are usually required.This inflates computational cost.Furthermore, [6] observed that FEM hardly evaluates dependent variable at the borderline of two adjacent elements.Hence, FEM may not be best suited for problems with non-trivial discontinuities.Although FEM solves some form of integral equations, these do not derive from the original form of conservation law.They are equivalence of the differential form of conservation law under Reynolds transport theorem.Other limitations of FEM in this regime are reiterated in [7].Despite these shortcomings; FEM remains an efficient numerical tool in computational physics, structural mechanics and environmental research.For example, using exponential element approach, [8] presents a novel FEM analysis of multiple porosity contaminant transport in porous media.However, most other works are based on available commercial FEM packages for contaminant tracking in geologic profiles.
In addition to domain discontinuity, there are instances in contaminant hydrology problems when the model equation alternates between parabolic and hyperbolic states.This depends on the instantaneous order of magnitude of the advection and diffusion speeds.This instability in physics often leads to severe numerical oscillations.In such cases, standard numerical schemes that are built on the differential forms of conservation law are known to experience some difficulties.To handle such difficulties, extensive studies have been carried out either to customize standard methods, or design new ones.However, most of the emerging schemes are designed to address specific types of interface, or resolve definite numerical problems.Notable on this list is the Euler-Lagrangian based schemes.They include localized adjoint method and Eulerian-Lagrangian localized adjoint method developed in [9].Although, the adjoint method and its derivatives are able to optimize flux dissipation, they do not enforce mass conservation and are deficient in ability to handle general boundary conditions.
Spatial and temporal scale-effects are additional difficulties with numerical simulation of contaminant transport in geologic media.They induce errors in contaminants profiles in aquifers.In a finite element analysis, [10] noted that scale effects are amplified at the interface between two strata even when the axis of flow is parallel to stratification.To handle scale-effect, [11] developed a stochastic mixed finite element technique that coarsens the flow equation in space.This is followed by piecewise parameter up-scaling using adaptive sparse grid collocation method.Clearly, the method does not completely account for parametric jumps that are common in real aquifers.On the same note, [12] used a strategic porous medium with high fracture frequency to present a numerical solution for colloid facilitated contaminant transport in aquifers.This minimized scale effect, but the approach is constrained by the high cost of computing.
For the foregoing reasons, we herein present the Space-Time Conservation Element/Solution Element numerical method as an efficient tool for contaminant tracking in fractured porous media.The CE/SE method is a collection of innovative numerical tools originated by [13,14] for solving conservation laws.It was developed from a perspective that is quite different in concept and methodology from traditional numerical schemes.Its development was motivated by the desire to build a general and coherent numerical framework that avoids many of the limitations of traditional numerical methods.
Rather than the differential form of the conservation laws, the CE/SE method derives its solution from the native Gauss integral form of the conservation law.As a result, local and global flux conservation is intrinsic to CE/SE formulations.Therefore, it is able to overcome most limitations encountered by numerical schemes that are built on the differential form of conservation laws.In effect, CE/SE solutions are naturally compliant with the physical model.Like other numerical schemes, the CE/ SE method uses an approximating function i.e. the Taylor series to describe numerical approximation of the exact solution on a computational mesh.The CE/SE mesh consists of Conservation Elements (CE) and Solution Elements (SE).A CE is a small region in two-dimensional Euclidean space E 2 , within which flux conservation is enforced.This ensures global flux conservation.On the other hand, a SE is another small region in E 2 within which unknown variables are approximated.

Transport Equations
In this work, we consider the hydro dispersive scheme of contaminant transport in two dimensions.Thus, three transport processes namely advection, dispersion and reaction are combined to describe the concentration profile i.e.C (y,z,t) of an aqueous contaminant in natural freshwater aquifer.The effects of sorption and retardation are also considered to give the governing differential equation as: where Here, u z and u y are the longitudinal and transverse velocities of fluid.Similarly, D z and D y are the corresponding dispersion coefficients, while R y and R z are the coefficients of retardation of the contaminant along y and z directions respectively. is the rate constant of the attendant geo-chemical and allied reactions. y and  z are coefficients of dispersivity, while D  is the diffusivity index of the medium.

The Physical Model
Next we simplify Equation (1) by using the dual porosity concept to delineate contaminant pathways into two regions.The momentum transport zone called the channel.Flow in the channel is coupled with transport processes such as sorption, diffusion, and reaction in the adjoining matrix.Thus, using some simplifying assumptions as detailed in [15], an iconic model for this system of flow is illustrated as in Figure 1.
We consider channels of thickness , bounded by vertical lines at − and the origin.The centerline of the channel is at −/2.The matrix separating two channels proceeds from the origin along the horizontal axis to 2.Therefore, the matrix centerline is at y = .Consequently, for this flow system, the equivalence of Equation (1) for transport in the channel and matrix is given by Equations (3) and (4) respectively. where ;0 ; with and

Model Simplification
For further simplification, we normalize flow parameters in Equations ( 3) and ( 4) with their corresponding scale factors.Then we group the resulting parameters into standardized dimensionless forms.This is followed by order of magnitude analysis following the approach of [ 15,16].Hence, Equations ( 3) and ( 4) are reduced to Equations ( 5) and ( 6) respectively.
Equation ( 6) derives from Equation ( 4) under the assumption that the Peclet number only applies to the channel.In this case, the reaction term is treated as a first order geochemical reaction.In addition, the following non-dimensional variables apply

Contaminant Loading and Boundary Conditions
Next, we consider the boundary and initial conditions needed to complement Equations ( 5) and ( 6).Here, we attempt to model the Nigerian experience of nitrate pollution of groundwater resources.Thus, we introduce a realistic and coherent boundary condition for the associated groundwater pollution problem.Precisely, we simulate what obtains in normal agricultural practice.In this case, seasonal loading of inorganic fertilizer with an ar-bitrary nitrate concentration C o is normally applied on topsoil.The surface behavior of the applied nitrate fertilizer is herein modeled as an integral Dirichelet/Neumann condition.This gives an exponentially decaying function of the form C o t e   .The function accounts for the progressive decay of the contaminant on topsoil as a result of some physico-chemical processes.With repeated seasonal (i.e.periodic) applications, the corresponding nitrate profile on topsoil can be represented as shown in Figure 2.
Here, β is the rate of disappearance of applied fertilizer on the soil surface.This may be due to percolation, evaporation, runoff etc. n is the number of cycles of fertilizer applications.On extension, the contaminant loading profile described above has widespread practical applications.With slight modification, it can be used to handle various types of physical loading conditions that are of interest to contaminant hydrologists and environmentalists.To be sure, the periodic, exponentially decaying function truly models actual physical conditions better than the impulsive, uniform and step function loadings that are well considered in literature.For instance, a leaking canister of toxic waste in a geological repository can be considered as having periodic effect on the surrounding aquifer, if for instance a similar canister has leaked in the region in the past or will most probably leak in the future.Such periodic loading profiles act to alter the hydrological properties of the aquifer and the initial concentration of the latest spill.
Our analysis is now reduced to solving simultaneous partial differential equations posed in Equations ( 5) and ( 6) to track the breakthrough curves of the contaminants in time and depth of an aquifer, under the normalized bo ndary condition: u Equation (8a) is a Fourier series representation of the loading pattern in fracture and matrix has  (10) Using this decompos n, the fracture flow Equation ( 5) splits into the form; itio

Development of CE/SE Scheme
By invoking Gauss divergence theorem, the conservation sed in law expres Equa (11a) can also be written as; tion and be sure, the integral equations expressed in Equations (12) and (13)  .j = 0, ±1/2, ±1, ±3/2, … ±J; k = 0, ±1/2, ±1, ±3/2, … ±K, and n = 0, 1/ 3/2, … N. At y = 0, j = 0 so we have the corresponding computational r mesh poi (0, k, n), as shown in Figure 3 In Figure 3(a), t hombus CD and ADEF constitute a pair of Conservation Elements namely CE + and CE − respectively.The common/adjoining boundaries of these conservation elements for the mesh point (0, k, n) define their solution ore, the mesh points are staggered so that the difference between the values of k and n at adjacent nodes alternates between whole numbers and half integers.

Derivation of the One-Dimensional Advection-Dispersion Scheme
Flux conservation demands that the total flux across a conservation element e.g.ABCD (CE + ) and in There is no temporal change on AB , so that In this regard, ( Using this procedure, the total flux across the bound ries of CE is given as g r g r r g r g r (16) Therefore, the condition for flux balance on t servation element is given in Equation ( 17 Normalizing Equation ( 17 The complementary part of Equation (18a) is obtained by enforcing flux conservation criterion on the c nant across the boundaries of CE − .This evaluates to ontami- Thus the solution; where Equation ( 18) is a one-dimensional CE/SE advectiondispersion (ad) scheme.a is the a and dvection speed 1 z Pe  ion (5 erica is the coefficient of dispersion.In uat ), a = 1.Special cases of the ad sc pu view of Eqheme are the re advection z Pe   , and the pure dispersion a = 0 schemes.Using the analysis of [13], the a-d scheme is num lly stable provided z  ≤ 1.However, the pure dispe n scheme is unconditionally stable.
we recall Equation ( 6) and introduce an arbitrary advection term with sp rsio

Construction of CE/SE Solution for Dispersion Reaction Transport
To build a CE/SE solution for the reaction dispersion equation in the matrix, and to streamline developm eed ent, y a , thus we have In Equation (21a) Setting 0 y a  in Equation ( 20) to obtain the dispersion reactio ort, we have n transp

Construction of CE/SE Solution for Fracture and Matrix Flow
 is the numerical analogue of the two dimendistribution of the contaminant in the fracture.The numerical equivalence of Equation (11) for the fracture flow can be represented as; Similarly the concentration distribution across the fracture-matrix interface is given as; and the concentration distribution in the matrix satisfies where , 1 2 , , , Here, we have chosen Δy the step size acr fracture-matrix interface as to stretch computation in the region.Followi results of [17], Equa-0) and ( 23

Validation of CE/SE Numerical Sch
To perform an indirect check on the accuracy and validity el em of n j H order in ac eme of CE/SE algorithms dev oped in Section 3, we consider as a test case the probl the transient distribution of charge carriers q(x, t) on the normalized base length of an electrical transistor satisfying the following transport model; Here  and k are material constants.They are the drift velocity and the diffusion coefficient of charge carriers along the cross section of the semi-conductor.Equation (29) (a, b and c) has a simple closed form analytical solution given in Equation ( 30) below.
Comparative performance of the CE/SE scheme relative to the exact solution and the quickest five point's finite difference scheme are shown on Figures 4 and 5.
The chosen values of a, k, x, and t are identical to the parameters employed for simulations in this work.
In Figure 4, there is a very good agreement between CE/SE and Finite Difference approximations with analytical solution.In this case, absolute errors are evaluated as 0.16% and 0.14% respectively.The corresponding Courant and diffusion numbers for the simulation are 2.56  10 −4 and 7.50  10 −3 respectively.However, in Figure 5, it is seen that the quickest five point finite difference scheme presented in [18] is marginally stable in this regime.The notation b  10 n implies b e  n.

Discussion of Results
Numerical simulations of fracture/porous matrix flow conditions are performed over an average of 20,000 static conservation elements.Except for redefining the computational step size to stretch the fracture matrix interface, no other mesh refinement algorithm is implemented hroughout the physical domain.Our simulation results t a re the time dependent concentration profiles of contaminant in the fracture and matrix.Simultaneously, we also obtain the time dependent flux variation of the contaminant.This flux variation result is new and it explains me phenomena identified in literature.The CE/SE emes described in this work are explicit solvers, hence it requires minimal computational resources.This is because it does not require global matrix assembling and inversion processes.

Response of Profiles in Fracture to Cyclic
Loading (Homogeneous Aquifer) Figures 6 and 7 show the concentration profile of the contaminant in silty clay and fissured limestone aquifers.These distributions are due to one, two and three cycles of fertilizer applications per annum for a period of 30 years.
As shown in Table 1, for these aquifers, we have used identical values of conductivities but different coefficients of dispersion.In line with expectation, concentration profiles in the fracture increase with the frequency o cycles per annum is less pronounced than rate one and two r that aquifer response to so sch f loading.On both figures, the difference between loading two and three the difference between loading at the cycles per annum.Thus, we infe high frequency of loading is similar to the asymptotic effect of constant loading.Also, these figures show that the sensitivity of the aquifers to periodic loading is more of a direct function of their dispersivity.Conversely, it can be deduced that the dispersion coefficient of a contaminant in an aquifer depends on its loading history.These figures further show that an increase in the dispersion coefficients in the fracture will accelerate solute arrival time.Theoretical basis for the wavy patterns of the conc tration profiles is provided in [19].In addition, these p correctly simulate the patterns of field of nitrate profile under agricultural land published in [20].This is the first work that recovers the field profile of nitrate under agricultural soil.
However, Figure 8 shows that for consolidated formations with mean conductivity less than 10 −6 , cyclic loading has  the other hand, when compared with Figure 7; Figure 9 illustrates the effect of loading duration on contaminant font in a relatively permeable fissured limestone aquifer.Sustained periodic loading rectifies the wavy pattern of nitrate profile and marginally increase its breakthrough depth.

Corresponding Response of the Matrix
Figure 10 shows representative concentration distribution at various depths in the matrix block separating tw fractur t two ycles per annum for thirty years.o es in a typical loose clay aquifer loaded a c Similarly, Figure 11 shows the corresponding concentration distribution at various depths in the matrix of limestone formation loaded at the same frequency over a period of hundred years.
The profiles in the matrix indicate storage in regions adjacent to the fracture.This delays the progress of the  depth of the aquifer.This suggests that the matrix has the tendency to feed the fracture at certain locations where advective transport is slower than diffusive transport or when flow in the fracture is off season.Thus, Figures 12 and 13 explain the observation made by [21] associated with the difficulty of finite element method to evaluate contaminant flux at the interface between fracture and matrix.Equally important remark is that the solute flux vanished with distance into the matrix.This behavior satisfies the no flux condition specified on the matrix centerline y =  by [13].This condition is also prescribed in our complementary boundary condition Equation (7d).These profiles further validate the consistency and efficacy of the CE/SE in this regime. .However, the lower consolidated formation acts to protect aquifers deeper than 37.5 meters from contamination.Loading frequency does not affect the amount of solute crossing the interface.Rather, the contaminant builds up in the upper stratum.This effect is well pronounced in Figure 14, which shows the time response for Loose Clay/Consolidated Clay stratified system.The presence of an accumulation zone at the interface between the two strata is evident.This effect is pronounced when the ratio of the conductivities or dispersion coefficients of the upper to the lower formation is of the order of 10  of the profiles at the interface between two strata.Clearly, the CE/SE method has an exceptional ability to handle stratification.The angle of refraction depends on the relative conductivity and dispersion coefficient of the neighboring strata.Also, mild variation in either the conductivity or dispersion coefficient of a geological system is sufficient to modulate stratification.

Effects of Stratification
On the other hand, Figure 16 illustrates the effects of the geological log on the breakthrough curves of contaminant in a stratified aquifer.The loading frequency for this study is two cycles per annum for duration of 50 years.In the 20A/80B composition, the upper stratum is heavily polluted.Consequently, water wells that are less than 20 meters in depth in such an aqu domestic consumption even whe ifer are unsafe for n the concentration of onding to 80A/20B and 20 for drinking so long as th la y depends on the ordering of the layers.Figure 17 is the outcome of numerical exercise to re-nitrate on topsoil is just 10.0 ppm.Also, if the solute concentration on topsoil is greater than 18.0 ppm, wells that are less than 37.5 m in depths will be unsafe for drinking using WHO standards of maximum of 10 ppm of nitrate in drinking water.
Similarly, in the 80A/20B composition the minimum depth to safe water is 60.0 meters if nitrate concentrations in excess of 50.0 ppm on topsoil.Clearly, the composition of a stratified aquifer determines its susceptibility to pollution.Next, we look at the effects of ordering of layers on relative susceptibility of stratified aquifers.The breakthrough curves corresp B/80A arrangement of layers in a binary system are shown in Figure 17.While depth to safe water lies below 60.0 meters for nitrate concentration in excess of 50.0 ppm on topsoil in the 80A/20B system, wells in the 20B/80A arrangement are safe eir bottom is deeper than 10.0 meters into the aquifer.Evidently, contaminant profile in multi-layered systems rgel place the 80A/20B and the 20B/80A systems with an equivalent aquifer whose hydraulic parameters are the weighted averages and mean hydrological parameters of constituents.Clearly, this approach is misleading due to sharp variations in properties.Similar results are obtained even when the conductivities are comparable but with significant difference in coefficien Similarly, Figure 18 describes the concentration distribution in cretaceous clay (u = 9.6e − 5 m/day, D = 5.43e − 3 m 2 /day) underlain with fissured limestone.The figure shows that mean value or weighted average approach to system's simplification in stratified aquifers can be used to optimize computational cost provided the conductivities and dispersion coefficients of adjacent layers are relatively close.

Conclusion
In recognition of the importance of protecting freshwater aquifers from aqueous phase liquid contaminants, this work presents the development and application of an a gure 18.Effects of parameters averaging in formations dvection-dispersion-reaction a-d-r Space-Time Conservation Element/Solution Element (CE/SE) numerical scheme for contaminant tracking in fractured porous media.Using flux splitting approach, a two-dimensional CE/SE scheme is developed for simulating the evolution and fate of aqueous phase liquid contaminant in stratified fractured porous media.Thus, the CE/SE scheme is herein extended to handle geological time scale problems.The developed scheme is able to resolve some of th lude the use of ideal initial/boundary conditions and parametric jumps that may be due to stratification or scale effects.To simulate the Nigerian experience of nitrate pollution of freshwater aquifers; we have deployed the developed scheme to solve the advection-dispersionreaction a-d-r equation in geologic media under a time periodic Dirichelet type boundary condition.This models the actual pattern of aqueous phase contaminant on arable agricultural lands.It is established that the CE/SE method is a viable and efficient tool for tracking contaminant transport in geologic profiles.The scheme is able to sense order of 10 −3 variation in hydrological properties of aquifers.It is also shown that attempts at systems simplification in heterogeneous reservoirs through the use of mean or weighted average values of hydrological parameters are in error unless the hydro-geological properties of neighboring geological formations ar is computationally inexpensive and easy programmed.

Figure 1 .
Figure 1.Flow geometry of contaminant motion in fractured porous medium.

Figure 2 .
It describes the time v tio   a) arian of nitrate on the soil surface due to periodic farming pattern.Here 2πk w   ;  is the dimensionless loading period which i n t is case as 365/T.T is the characteristic time ant transport through the aq-s given i h of contamin c.Boundary layer between the uniform concentration d.Symmetry along the y =  line sugges uifer to the water table.The complementary boundary conditions associated with the contaminant loading profile in Equation (8a) are as follows: b.Contaminant concentration decays with depth in the aquifer as defined by components of the contamin stribution in the fractur the vertical and horizon irections respectively.Here, to facili eric duce the ion s e along tate num al development, we have re-introadvect peed a.However, Equation (6) retains its original form.It is applied only in the range 0 y     .

Figure 3 (
a) should satisfy Equation (12 quently, we compute the flux of the contamin oundaries of the conservation elements.

Figure 4 .
Figure 4. Comparison of CE/SE and finite difference approximations.

Figure 5 .
Figure 5.Comparison of CE/SE and finite difference approximations.

in and 15 Figure 14
Figure14represents the response of a 50% Loose/50% Consolidated Clay system to 1, 2 and 3 cycles/annum loadings over 50 years duration.The behavior of the upper stratum is identical to what obtains in homogeneous loose clay aquifer; shown in Figure6.However, the lower consolidated formation acts to protect aquifers deeper than 37.5 meters from contamination.Loading frequency does not affect the amount of solute crossing the interface.Rather, the contaminant builds up in the upper stratum.This effect is well pronounced in Figure14, which shows the time response for Loose Clay/Consolidated Clay stratified system.The presence of an accumulation zone at the interface between the two strata is evident.This effect is pronounced when the ratio of the conductivities or dispersion coefficients of the upper to the lower formation is of the order of 10 p with p  1.0.notherimportant feature of the breakthrough curves A Figures 14

Figure 15 .
Figure 15.Effect of loading in stratified aquifer.

Figure 16 .Figure 17 .
Figure 16.Effects of stratification of aquifer on profile of breakthrough curve.
numerical solution of groundwater contaminant hydrology problems.Such difficulties inc e h very close.The developed model was validated wit analytical and the quickest five points finite difference cheme.CE/SE s il are different from the Eulerian La-