Numerical Modeling of Contaminant Transport with Spatially-Dependent Dispersion and Non-Linear Chemical Reaction

Abstract. A one-dimensional advective-dispersive contaminant transport mo del with scale-dependent dispersion coefficient in the presence of a nonlinea r chemical reaction of arbitrary order is considered. Two types of variations of the dispersion coefficient with the downstream distance are considered. The first type assumes that the d ispersivity increases as a polynomial function with distance while the other assumes an exponentia llyincreasing function. Since the general problem is nonlinear and posses ses no analytical solutions, a numerical solution based on an efficient implicit iterative tri-dia gon l finitedifference method is obtained. Comparisons with previously published an alytical and numerical solutions for special cases of the main transport equation ar e performed and found to be in excellent agreement. A parametric study of all physical pa r meters is conducted and the results are presented graphically to illustrate interesting f eatures of the solutions. It is found that the chemical reaction order and rate coefficie nt have significant effects on the contaminant concentration profiles. Furthermore, the sc al -dependent polynomial type dispersion coefficient is predicted to obtain significant ch anges in the contaminant concentration at all dimensionless time stages compared with th e constant dispersion case. However, relatively smaller changes in the concentra tio level are predicted for the exponentially-increasing dispersion coefficient.


Introduction
The problem of contaminant transport in soil, groundwater and surface water has been a research subject of many recent and old theoretical and experimental investigations.This is due to increased public awareness of significant contamination of groundwater and surface water by industrial, municipal, agricultural chemicals, accidental spills and effect of soil contamination resulting from landfills and burying of hazardous materials.The principle differential equation governing solute transport and chemical reactions has been developed using mass balance and advective-dispersive principles [1] and is widely used in modeling solute transport phenomenon.A number of analytical solutions for steady-state flow and different boundary and initial conditions were given by van Genuchten and Alves [2] for problems with linear adsorption and zero-and first-order production and decay.Analytical solutions play an important role in modeling because they offer fundamental insight into governing physical processes, provide useful tools for validating numerical approaches, and are sometimes more computationally efficient [3].
Most previously published analytical solutions to advective-dispersive transport problems are obtained based on the assumption of a homogeneous porous medium [2].In reality, subsurface porous media through which the contaminant moves are seldom homogeneous and significant spatial variability of transport properties should be expected [3][4][5].
As a result of the heterogeneity of the porous media, the dispersion coefficients in all directions vary with the space coordinate and the resulting contaminant transport equation contains spatially-dependent coefficients.Limited analytical solutions for scaledependent dispersion coefficients have been reported in the literature.Yates [6,7] obtained one-dimensional solutions for uniform flow with constant concentration or constant concentration flux boundary conditions when the medium has a linearly or exponentiallyincreasing dispersion coefficient with the spatial coordinate.Huang et al. [8] also presented analytical solutions for a scale-dependent dispersion coefficient which increases linearly with the distance until some distance after which it reaches an asymptotic value.Logan [9] derived an analytical solution for the one-dimensional equations incorporating rate-limited sorption and first-order decay under time-varying boundary conditions, assuming an exponentially-increasing dispersion coefficient.Zi-ting [10] reported an analytical solution for an exponential-type dispersion process.However, the solutions given by Yates [6,7], Huang et al. [8], Logan [9] and Zi-ting [10] are complex and difficult to evaluate numerically.It is worth noting that other studies of solute transport employ time-dependent coefficients.Warrick et al. [11], Barry and Sposito [12] and Basha and El-Habel [13] reported analytical solutions to the one-dimensional advection-dispersion equation with arbitrary time-dependent dispersion and velocity coefficients.
This paper presents numerical solutions for one-dimensional contaminant transport through a semi-infinite porous medium domain in the presence of a nonlinear chemical reaction.The transport starts from a continuous contaminant source and the mechanical dispersion effect is assumed to vary with the downstream distance.The solutions include first-and second-order homogeneous irreversible chemical reactions as well as polynomial and exponentially-increasing spatially-dependent dispersion coefficients.

Formulation
Consider transient one-dimensional advective-dispersive contaminant transport in a porous medium from a continuous source with a non-linear chemical reaction.The movement of the contaminant takes place in the semi-infinite region 0 ≤ x < ∞ and the dispersion coefficient is assumed to be spatially-dependent.The governing equation for this situation can be written as: where t is time, x is the one-dimensional spatial coordinate (0 ≤ x < ∞) (or longitudinal distance), c(x, t) is the concentration, u is the uniform velocity, D(x) is the spatially variable hydrodynamic dispersion coefficient representing the sum of both the effects of molecular diffusion (D e ) and mechanical dispersion (αu where α is usually called the dispersivity), k is the chemical reaction rate coefficient and n is the order of the homogeneous irreversible chemical reaction.It should be noted that the first term on the left hand side of equation ( 1) represents the transient or accumulation effect.The second term on the left hand side of the equation represents the advection or convection effect which is defined as the transport of contaminant by the mean velocity in the flow stream.The first term on the right side of the equation accounts for the dispersion or diffusion effect which is responsible for the spreading of the contaminant in the medium.The last term of equation (1) represents the nonlinear reaction effect which may take place (depending on the nature and properties of the contaminant) between the contaminant and the medium.The initial and boundary conditions for this problem are where c 0 represents a constant continuous concentration source.Equation ( 1) can be written as It is convenient to work with dimensionless equations.This can be accomplished by using where x 0 and D 0 are characteristic longitudinal distance and dispersion constant, respectively.
Substituting equation (4) into equations ( 3) and (2) gives, respectively where are the Peclet number and dimensionless chemical reaction rate constant, respectively.

Numerical method
In its most general form, equation ( 5) is nonlinear.Therefore, an analytical solution to this equation is unlikely and a numerical procedure is required.Many existing computer codes employ a finite-difference approach for the solution of transport equations.It is logical to investigate the applicability of this methodology to equation ( 5).In the present work, an implicit iterative tri-diagonal finite-difference method similar to that discussed by Blottner [14] is employed.A two-point backward difference quotient is used to represent the dimensionless time τ derivative and three-point central difference quotients are used to represent the dimensionless space η derivatives.The computation starts at τ = 0 and marches forward in time.At each time, a system of non-linear algebraic equations must be solved to determine the η distributions of C.An iterative procedure is employed for this purpose.At each iteration, an equivalent linear system of algebraic equations (the linearization being effected by representing some quantities by their values from the previous iteration) must be solved.These equations have a tri-diagonal form and can be solved by the Potter's method variables which can be determined by a forward sweep in the h direction.Then the physical variables can be found from a corresponding backward sweep.This process avoids the need for matrix inversion.Iteration is continued until convergence is obtained at a given time.The procedure moves forward for the next time.
It is helpful to have some analytical solutions for special cases of equation ( 5) to use as standards of comparison for the numerical procedure.

Analytical solutions
Consider the special case where D * (η) = 1 (constant dispersion) and n = 1 (first-order chemical reaction) for which equations ( 5) and ( 6) are simplified to read These equations are linear and can be solved analytically by the Laplace transformation methods.Without going into detail, the solution of the above initial-value problem can be shown to be where erf and erfc are the error function and complimentary error function, respectively and θ is a dummy variable.It should be noted that equation ( 10) is consistent with and represents the dimensionless form of that reported earlier by van Genuchten and Alves [2].As Pe → ∞ (no diffusivity), equation (10) reduces to where H is the unit step function.For small amounts of diffusivity (1/Pe 1) the discontinuity exhibited by equation ( 12) at η = τ is replaced by a narrow continuous transition layer.
For the special case in which λ = 0 (no chemical reaction), equation ( 8) becomes a simple convection-diffusion equation.Equations (10) and ( 12) are valid with λ = 0.This leads to the respective results: The accuracy of the numerical method discussed above is validated by direct comparisons with the analytical results given in equations ( 10) and ( 13).These comparisons are presented in Figs. 1 and 2, respectively.It is clear from these figures that excellent agreement between the numerical and analytical results exists at all presented dimensionless times for both Pe = 1 and Pe = 100 considered in these figures.In order to check the accuracy of the numerical results for a nonlinear chemical reaction, further comparisons are performed with the work of Onyejekwe [15] who reported the solution of a single-phase isothermal flow with nonlinear kinetics involving one reactant.The governing equation for this problem in an idealized one-dimensional finite region is given by the following transport equation and conditions: The initial-value problem given by equations ( 15) and ( 16) was solved numerically by the implicit finite-difference scheme discussed above.The obtained results for the exit concentration were compared with those reported by Dale [16], Ramachandran [17] and Onyejekwe [15] for different values of n, Pe and λ.These various comparisons are shown in Tables 1 through 3. Again, these comparisons show good agreement except for high values of Pe for which Onyejekwe's [15] method seem to under predict the exit concentration considerably and blows up for Pe = 1000.These discrepancies are probably due to the inaccurate evaluation of the closed-form solutions who involve complicated functions reported by Onyejekwe's [15] for large values of Pe as evident from the fact that his solution blows up for Pe = 1000.The various favorable comparisons reported in Figs. 1 and 2 and Tables 1 through 3 lend confidence in the numerical results to be reported in the next section.

Results and discussion
Numerical solutions for the general advective-dispersive-reactive contaminant transport equation ( 5) subject to the initial and boundary conditions ( 6) are obtained for two different types of scale-dependent dispersion coefficients.These are where all of N, n 1 , a and b are dimensionless constants.It should be noted that when N = 0 in equation ( 17) and a = 0 in equation ( 18), the constant dispersion cases are recovered.It should be noted herein that Zoppou and Knight [18] derived an analytical solution for a transport problem with variable velocity and diffusivity.In the notation of this work, Zoppou and Knight [18] assumed that D * (η) = η 2 .Also, Zi-ting [10] used the exponential-type dispersion coefficient given by equation (18).The computational domain was divided up into 500 points in the η direction and 600 points in the τ direction with variable step sizes in both directions.The initial step sizes and growth factors employed in the η and τ directions were 0.001, 0.001, 1.055 and 1.03, respectively.In this case, η max = 4 × 1010 represented the condition η → ∞.These values were arrived at after performing various numerical experiments to access grid-independent results.The convergence criterion required that the difference between the current and previous iterations must be 10 −7 .
Various numerical results are obtained and a representative set of results is presented in Figs. 3 through 12.These results are chosen to illustrate the influence of the chemical reaction order n, the Peclet number Pe, the chemical reaction constant λ, and the scale-dependent dispersion constants n 1 and a.In Figs. 3 through 6, the dispersion coefficient is constant while in Figs.7 through 12, the dispersion coefficient is spatially variable.Fig. 3 presents the temporal development of the contaminant concentration profiles for various values of the chemical reaction order n (corresponding to first-, second-and third-order reactions) for the case of constant dispersion effects with Pe = 1.In this and all subsequent figures, the space coordinate is represented by a logarithmic scale so as to capture the complete transition from unsteady conditions at small time values of the dimensionless time τ to steady-state conditions at τ = ∞.In general, it is predicted that increases in the chemical reaction order increases the concentration.This is accompanied by an increase in the concentration boundary-layer thickness.The increase in the concentration field and its boundary layer appears to be more significant at larger time values especially at steady-state conditions.Physically, the increase in the concentration boundary-layer thickness as n increases means that it a solute transport with a first-order chemical reaction reaches its steady-state conditions at a faster rate than it would with a higher chemical reaction order.Fig. 4 displays the same parameters as in Fig. 3 except the value of the Peclet number which is set to 100 representing a small dispersive effect.The same general conclusion is obtained is which the contaminant concentration profile and its boundary-layer thickness tend to increase as n increases.Also, effect of increasing n is more pronounced at the steady-state conditions.However, by direct comparison with Fig. 3 for Pe = 1, it can be seen that the effect of the chemical reaction order is more for contaminant transports with higher dispersive effects.
Fig. 5 illustrates the influence of the Peclet number Pe on the concentration profile for a contaminant transport with a second-order reaction at three different time values corresponding to early time (τ = 0.6), intermediate time (τ = 52) and state-state conditions (τ = ∞).It is predicted that higher dispersive effects (low values of Pe) provide smoothing effects in the concentration profiles.Also, increasing the dispersive effects (that is, decreasing the values of Pe) has the tendency to increase the ability of the contaminant to transport easier through the porous medium.This is reflected in the increases in the values of C as Pe decreases.In addition, as Pe decreases, the concentration boundary-layer thickness increases and this seems to be more pronounced at early transport time stages.It should also be noted that at early time stages (τ = 0.6), the concentration increases as Pe decreases every where except in the immediate vicinity of the inlet boundary (η = 0).This is because as Pe → ∞ (no dispersive effects), it is expected that the concentration profile to drop sharply (a step function) to the terminal condition as η → ∞.Fig. 6 depicts the influence of the chemical reaction constant λ on the temporal development of the concentration profiles for a second-order reaction and Pe = 1.Physically, the chemical reaction term in equation ( 5) represents a concentration decay or sink term.This means that for a specific reaction order, as λ increases, the decaying effect increases causing the contaminant concentration to decrease everywhere in the flow region away from the boundaries and for all times except the early time stages.It is worth noting that the decaying effect is much more significant at the steady-state conditions than at all other time conditions.In addition, the concentration boundary-layer thickness decreases as λ increases.for n 1 = 2, the steady-state condition is reached at a faster rate than that observed for the linear case where n 1 = 1.Specifically, for n 1 = 2 the steady-state conditions are achieved already at τ = 12.Fig. 8 presents the same parameters as in Fig. 7 except for the value of n which is set to 2 representing a second-order chemical reaction.As is obvious from this figure, the same conclusion as in Fig. 7 is reached.That is, increasing the value of n 1 causes reductions in the contaminant concentration level everywhere except in the region close to the end boundary where it increases resulting in increases in the concentration boundarylayer thickness.Also, for n 1 = 2 the steady-state conditions are achieved at τ = 12.
Fig. 9 shows the effect of increasing n 1 on the contaminant concentration for a second-order chemical reaction and Pe = 100.Again, in general, the concentration   level decreases in most of the domain except near the downstream boundary where it increases as n 1 increases.The changes in the concentration profiles are more prominent at the steady-state conditions than at the earlier time stages of the flow.Figs. 10 through 12 present the effects of the constant of exponential-type dispersion coefficient a on the unsteady concentration profiles for the cases of first-order chemical reaction with Pe = 1, second-order chemical reaction with Pe = 1, and second-order chemical reaction with Pe = 100, respectively.It is observed from these figures that increasing the value of the constant a causes decreases in the contaminant concentration level and in the concentration boundary layer thickness.Also, this decrease is greater at small time stages than it is at higher time stages.In addition, as expected, the effect of increasing the constant a is very little for high values of Pe (small dispersion effects).These trends are clearly depicted in Figs. 10 through

Conclusion
A one-dimensional advective-dispersive contaminant transport model with scale-dependent dispersion coefficient in the presence of a nonlinear chemical reaction was considered.The scale-dependent dispersion coefficient was used to characterize dispersion in a heterogeneous porous medium.Two types of variations of the dispersion coefficient with the downstream distance were considered.The first type assumed that the dispersivity increased in a polynomial function with distance while the other assumed an exponentiallyincreasing function.The nonlinear chemical reaction assumed an arbitrary reaction order.Since the general problem was nonlinear and possessed no analytical solutions, a numerical solution based on an efficient implicit iterative tri-diagonal finite-difference method was obtained.The accuracy of the numerical method was validated by various favorable comparisons with known analytical solutions and reported numerical solutions for special cases of the main transport equation.Several numerical solutions based on the general model were reported assuming a uniform flow field.A parametric study was conducted and the results were presented graphically to illustrate interesting features of the solutions.It was found that the chemical reaction order and rate coefficient had significant effects on the contaminant concentration field especially at the steady-state conditions.It was predicted that as the chemical reaction order increased, the contaminant concentration increased.On the other hand, increases in the chemical reaction rate coefficient produced reductions in the contaminant concentration level.In addition, as the Peclet number was increased, the concentration level was decreased.Furthermore, the scale-dependent polynomial-type dispersion coefficient was predicted to obtain significant changes (reductions) in the contaminant concentration at all time stages compared with the constant dispersion case.However, relatively smaller changes (reductions) in the concentration level were predicted for the exponentially-increasing dispersion coefficient.

Fig. 7
Fig.7displays the effect of the dispersion coefficient power exponent n 1 for the polynomial-type dispersion coefficient on the temporal development of the contaminant concentration profiles for a first-order chemical reaction.In general, as the power exponent n 1 increases, the concentration level decreases everywhere except in the region close to the end boundary where it increases causing the concentration boundary-layer thickness to increase.This behavior takes place for almost all time stages.It should be noted that

Table 3 .
Onyejekwe (1997)merical and approximate exit concentration for a single phase reactor withOnyejekwe (1997)with λ = 2.5 and n = 0.5 12.Fig.12. Effects of exponential dispersion constant on the temporal development of concentration profiles.