Groundwater pollution may occur due to human activities, industrial effluents, cemeteries, mine spoils, etc. This paper deals with one-dimensional mathematical modeling of solute transport in finite aquifers. The governing equation for solute transport by unsteady groundwater flow is solved analytically by the Laplace transform technique. Initially, the aquifer is subjected to the spatially dependent source concentration with zero-order production. One end of the aquifer receives the source concentration and is represented by a mixed-type boundary condition in the splitting time domain. The concentration gradient at the other end of the porous media is assumed to be zero. The temporally dependent velocity and the dispersion coefficients are considered. A numerical solution is obtained by using an explicit finite difference scheme and compared with the analytical result. Accuracy of the solution is discussed by using the root mean square error method. Truncation error is also explored for the parameters like numerical dispersion and velocity terms. The impact of Peclet number is examined. For graphical interpretation, unsteady velocity expressions (i.e., such as exponential, sinusoidal, asymptotic, and algebraic sigmoid) are considered. The work may be used as a preliminary predictive tool for groundwater resource and management.

INTRODUCTION

In India, many aquifers are being contaminated by a host of human activities, such as sewage disposal, refuse dumps, pesticides and chemical fertilizer contamination, industrial effluent discharges, and toxic waste disposal (Rausch et al. 2005; Batu 2006). The traditional advection dispersion equation is based on the conservation of mass and Fick's law of diffusion (Fried & Combarnous 1971; Bear & Verruijt 1987; Chrysikopolous et al. 1990) and constitutes the basis of solute transport models that are used for predicting the movement of contaminants in groundwater systems.

There has been some research on solute transport in groundwater systems. Hunt (1978) analyzed one-, two-, and three-dimensional solutions for instantaneous, continuous, and steady state pollution sources in uniform groundwater flow. Freeze & Cherry (1979) provided a relation between dispersion and groundwater velocity in which the dispersion is proportional to a power of the velocity and experimentally observed that the power ranged between 1 and 2. van Genuchten (1981) explored derivations of analytical solutions using the Laplace transform for the solute transport equation with zero-order production and/or first-order decay subjected to first and third type boundary conditions. Zoppou & Knight (1994) evaluated analytical solutions that are still useful for validating numerical schemes for solving the advection–diffusion equation with spatially variable coefficients. Logan (1996) obtained analytical solutions for a scale-dependent dispersion coefficient increasing exponentially with position up to some constant limiting values. Hantush & Marino (1998) developed analytical solutions using the Laplace and Fourier transform methods and superposition principle for the first-order rate model in an infinite porous medium.

Using Peclet and Courant numbers and a new sink/source dimensionless number, Ataie-Ashtiani et al. (1999) discussed truncation errors associated with finite difference solutions of the advection–dispersion equation (ADE) with first-order reaction. Bedient et al. (1999) presented a mathematical model of the ADE for describing the migration and fate of pollutants in groundwater. Neville et al. (2000) presented semi-analytical solutions for a multi-process non-equilibrium model for describing contaminant concentration distribution patterns. Balla et al. (2002) presented a computational case study using a transport model for pollution of underground water due to damage of the waterproofing system in a waste material depository or sewage sludge composting plant.

Lowry & Li (2005) discussed an improved finite analytical solution method for solving the time-dependent ADE that does not discretize the derivative terms rather solving the equation analytically in the space–time domain. Smedt (2006) presented an analytical solution for solute transport in rivers, including the effect of transient storage and first-order decay. Tkalich (2006) explored the derivation of high-order advection–diffusion schemes by employing the interpolation polynomial method. Hill et al. (2007) proposed upscaling models of solute transport in porous media through genetic programming in heterogeneous porous media. Chen (2007) presented an analytical solution of two-dimensional non-axisymmetrical solute transport in a radially convergent flow tracer test with a diffusion coefficient increasing with travel distance. Yeh & Yeh (2007) derived solutions of the transport equation with a point-source term considered as the point-source solution under the condition that the solute was introduced into the flow system from the boundary that was considered as the boundary-source solution. Kumar et al. (2008) also described transport through a heterogeneous porous medium with a time-dependent dispersivity in solute transport modeling. Zhan et al. (2009) explored two-dimensional solute transport in an aquifer–aquitard system by maintaining the mass conservation at the aquifer–aquitard interface. Gao et al. (2012) explored a mobile–immobile model with an asymptotic dispersivity function of travel distance with the concept of scale-dependent dispersion during solute transport in finite heterogeneous porous media. Rezaei et al. (2013) derived a semi-analytical solution to the two-dimensional conservative solute transport in an aquifer bounded by thin aquitards in the Laplace domain. Singh & Das (2015) explored the analytical and numerical solutions for one-dimensional scale-dependent solute transport in heterogeneous media in which analytical and numerical solutions were compared and found very good agreement among them. The root mean square (RMS) error analysis was made to check the accuracy of the solution.

In the present work, we focus on one-dimensional solute transport modeling using the ADE in a finite aquifer with first-order decay and zero-order production. To simplify ADE, different transformations were applied. Non-dimensional parameters were employed for reducing the number of parameters of the ADE. To predict the pattern of contaminant concentration, different types of unsteady velocities, such as sinusoidal, exponentially decreasing, asymptotic, and algebraic sigmoid, were considered. They helped describe the nature of contaminant concentration in time and space.

MATHEMATICAL FORMULATION

The contaminant transport phenomena in a groundwater layered medium or reservoir, i.e., aquifer, can be predicted by analytical solute transport modeling of the ADE. The modeling depends on the distance as well as time. The dominant process of mass transfer in groundwater is the advection and diffusion; in contrast, it refers to the solute transport through the action of random motions in the groundwater reservoir. The one-dimensional ADE for an isotropic homogeneous finite aquifer with zero-order production and first-order decay can be expressed as: 
formula
1
Here, is the longitudinal dispersion coefficient (i.e., representing longitudinal dispersion), is the volume averaged dispersing solute concentration in the liquid phase, is the volume averaged dispersing solute concentration in the solid phase, is the unsteady uniform downward pore seepage velocity, is the longitudinal direction of flow, is time, is the zero order production rate coefficients for solute production in the liquid phase, is the first-order decay rate coefficient in the liquid phase, is the bulk density of the porous medium, and n is the porosity of geological formation.
As the contaminant goes from the solid phase into the liquid phase under the linear isotherm condition, we can express mathematically: 
formula
2
where is the distribution coefficient, defined as the concentration adsorbed by the solid phase to the liquid phase into the groundwater reservoir, i.e., aquifer.
Equation (1) was solved analytically with appropriate initial and boundary conditions. As the aquifer is not solute free initially, i.e., at , a linear combination of initial concentration and a function of space-dependent along with zero-order production term is considered. The loss or gain of solute mass occurs due to chemical components within the liquid, radioactive decay and biodegradation, and growth due to bacterial activities. The loss of gain term is usually described through the sources of solute. The solute can grow in liquid phase and the solid phase. Its growth is frequently expressed by a zero-order production. In this study, initially the growth of solute along the space was a linear combination of the initial concentration taken into consideration. This can be written as: 
formula
3
Due to the increasing human activity at the earth surface, the solute concentration in groundwater increases in time. Hence, a mixed type boundary condition in the splitting time domain at the source is considered as follows: 
formula
4
Due to no mass flow at the other end of the domain, a flux type boundary condition is considered as follows: 
formula
5
Using Equation (2), Equation (1) can be written as: 
formula
6
where 
formula
7
The dispersion theory, proposed by Freeze & Cherry (1979), was employed here. As dispersion is directly proportional to the seepage velocity, i.e., , and so, , where, A is constant that depends upon the pore geometry of the aquifer. Let and , where, and are the initial seepage velocity and initial dispersion coefficient, respectively. Also, let and , where, and are the initial first-order decay rate coefficient and the initial zero-order production rate coefficient for solute production in the liquid phase, respectively. is the non-dimensional expression where is the flow resistance coefficient. Using this concept, Equation (6) can be written as: 
formula
8
A new time variable (Crank 1975): 
formula
9
is introduced and Equation (8) can be written as: 
formula
10
The following non-dimensional variables: 
formula
11
are used to reduce the number of parameters of Equation (10) and one can then write: 
formula
12
where 
formula
13
The initial and boundary conditions given in Equation (3), Equation (5), and Equation (6) can be written in non-dimensional form, respectively, as follows: 
formula
14
 
formula
15
 
formula
16
Now, the following transformation: 
formula
17
was used to remove the advection term from Equation (12). The initial and boundary conditions accordingly were transformed and the Laplace transform technique was employed. The analytical solution can be obtained as follows (see Appendix, available with the online version of this paper): 
formula
18
 
formula
19
where 
formula
20a
 
formula
20b
 
formula
20c
 
formula
20d
 
formula
20e
 
formula
21a
 
formula
21b
 
formula
21c
 
formula
21d
 
formula
21e
 
formula
21f
 
formula
21g
 
formula
21h
 
formula
21i
 
formula
21j
 
formula
21k
 
formula
21l
 
formula
21m
 
formula
22a
 
formula
22b
 
formula
22c
 
formula
22d
 
formula
22e
 
formula
22f
 
formula
22g
 
formula
22h
 
formula
22i
 
formula
22j
 
formula
22k
 
formula
22l
 
formula
23a
 
formula
23b
 
formula
23c
 
formula
23d
 
formula
23e
 
formula
23f
 
formula
23g
 
formula
23h
 
formula
23i
 
formula
23j
 
formula
23k
 
formula
24a
 
formula
24b
 
formula
24c
 
formula
24d
 
formula
24e
 
formula
24f
 
formula
24g
 
formula
24h
 
formula
24i
 
formula
24j
 
formula
24k
 
formula
25a
 
formula
25b
 
formula
25c
 
formula
25d
 
formula
25e
 
formula
25f
 
formula
25g
 
formula
25h
 
formula
25i
 
formula
25j
 
formula
25k
 
formula
25l
 
formula
25m
 
formula
26a
 
formula
26b
 
formula
26c
 
formula
26d
 
formula
26e
 
formula
26f
 
formula
26g
 
formula
26h
 
formula
26i
 
formula
27a
 
formula
27b
 
formula
27c
 
formula
27d
 
formula
27e
 
formula
27f
 
formula
27g
 
formula
27h
 
formula
28a
 
formula
28b
 
formula
28c
 
formula
28d
 
formula
28e
 
formula
28f
 
formula
28g
 
formula
29
 
formula
30
 
formula
31

NUMERICAL SOLUTION

The finite-difference technique has been applied in various works for numerically modeling the solute transport phenomena in groundwater systems (i.e., Molenkamp 1968; van Genuchten & Gray 1978; Richard 1988; Hogarth et al. 1990; Dudley et al. 1991; Moldrup et al. 1996; Zheng & Bennett 2002; Meerschaer & Tadjeran 2004; Ataie-Ashtiani & Hosseini 2005; Wang et al. 2005; Mohebbi & Abbaszadeh 2013). In order to solve numerically the ADE (Equation (12)) together with the initial and boundary condition transformed in domain , a suitable space transform was used as: 
formula
32
Equation (12) together with the initial and boundary condition may be written as: 
formula
33
 
formula
34
 
formula
35
 
formula
36
A finite difference technique was derived by using Taylor's series expansion (Mickley et al. 1957; Carnahan et al. 1969). The explicit finite difference scheme is commonly used, even though it may require an extended computing time because of its restrictive stability criterion. Here, we used a general form of the explicit finite difference approximation with forward time and central space forward difference scheme to approximate Equations (33)–(36). We obtained as follows: 
formula
37
 
formula
38
 
formula
39
 
formula
40
where superscript j refers to time, subscript i refers to space, is the time increment, and is the space increment.
The space domain Z and time domain T are discretized by a rectangular grid of points with mesh size and , respectively. So, one can write as follows: 
formula
 
formula
The contaminant concentration at a point in space with the jth subinterval of time T was defined by .

STABILITY CONDITION

The finite difference scheme is convergent if the discretization error approaches zero as the grid spacing and tends to zero (Carnahan et al. 1969). The stability test of the finite difference scheme was proposed by a matrix method (Smith 1978) and this was used by Notodarmojo et al. (1991,). Kwok (1992) investigated the stability properties of the various two-level, six-point finite difference schemes for the approximation of the convection–diffusion equation. The solution was convergent, subject to the satisfaction of the stability criterion. The finite difference scheme of the governing partial differential Equation (37) can be written as follows: 
formula
41
where 
formula
42
 
formula
43
and 
formula
44
Equation (41) was expressed in matrix form as: 
formula
45
where A contains the entire constants.

The difference approximation equation was stable if the eigenvalues of A had modulus values less than or equal to unity, i.e., , where was the eigenvalue of matrix A.

To find the bounds of the eigenvalues of matrix A on applying the Gerschgorin circle method, the stability criterion for the time step was as follows: 
formula
46

TRUNCATION ERROR

Numerical dispersion was first quantified by Lantz (1971). Ataie-Ashtiani et al. (1999) explored the expansion of the Taylor series of solute concentration along the ADE used for determining the truncation error in one dimension. Chaudhari (1971) investigated a second-order error through the examination of the truncated Taylor series approximation with explicit finite difference solution of the one-dimensional ADE. We also explored the truncation error for the various parameters, such as dispersion, seepage velocity, first-order decay, and zero-order production term.

From Equation (37), we obtained as follows: 
formula
47
From the Taylor series expansion of each term of Equation (47), we obtained as follows: 
formula
48
 
formula
49
 
formula
50
where is simply denoted by C.
After imposing the Taylor series expansion on Equation (47), we got the truncation error of the finite difference approximation of order . The transport parameters were constant within each combination of time and space increments in finite difference calculations. The second-order temporal derivative of C was written in terms of the spatial derivative of C by using the partial differential equation obtained from Equation (47). We obtained as follows: 
formula
51
Then the partial differential equation obtained from Equation (47) can be written as: 
formula
52

Now, after comparison between Equation (52) and the original partial differential equation, we found different forms of the truncation error, as discussed by Ataie-Ashtiani et al. (1996). These errors can be identified as follows:

  • The second-order truncation error or numerical dispersion was 
    formula
    53
  • The first-order truncation error or numerical seepage velocity was 
    formula
    54
  • The zero-order truncation error or numerical first-order decay was 
    formula
    55

Constant error or numerical zero-order production was 
formula
56
After removing the induced numerical errors from the finite difference model, Equation (52) can be written as follows: 
formula
57
where 
formula
58
 
formula
59
 
formula
60
 
formula
61

ACCURACY OF THE SOLUTION

The accuracy of the solution was obtained by comparison of the numerical result with the analytical one. In the numerical solution, the accuracy is the degree of closeness of concentration values of the numerical result obtained with various methods to those of the analytical result. Towler & Yang (1979) adopted a criterion of comparison that was more systematic and consistent: (1) the RMS error and (2) the absolute maximum error between the analytical solution and the numerical solution at all grid points. Roberts & Selim (1984) used the RMS method to calculate the average error at each nodal point of the grid. Singh & Das (2015) explored the accuracy of solution of the solute transport equation in comparison of the analytical result with numerical one. For testing the accuracy of solutions in this paper, we used the RMS error which is the most appropriate method to check the accuracy of the solution (Chai & Draxler 2014). The RMS method was used to calculate the average error at each point which is defined by: 
formula
62
where 
formula
The difference between the analytical and numerical concentration values is denoted by , and N is the number of data which were used to evaluate the accuracy of the solution.

RESULTS AND DISCUSSION

The analytical solution, obtained by Equations (18) and (19), was computed for the following data (Singh & Kumari 2014): 
formula
The average porosity of the different geological formation was considered as n = 0.37 (sand), 0.55 (clay) (Freeze & Cherry 1979). We considered four different time-dependent forms of velocity expressions that can be written as follows:
  1. Exponential decreasing form of velocity: 
    formula
    63
  2. Sinusoidal form of velocity: 
    formula
    64
  3. Asymptotic form of velocity: 
    formula
    65
  4. Algebraic form of velocity: 
    formula
    66
where m is the flow resistance coefficient and is the constant parameter. The first and third ones have been used by Aral & Liao (1996), the second one by Singh et al. (2009), and the last one is based on the properties of the algebraic sigmoid function which include the error function. It starts to progress from a small beginning, accelerates in the rainy season, and then reaches a limit over a period of time.
From Figure 1 it can be seen that the maximum contaminant concentration is observed, i.e., 0.8 in the case of aquitard (i.e., clay) and 0.6 in the case of aquifer (i.e., sand) which are subsequently lowered down to minimum concentration tends to zero at the far end of the domain with respect to the distance, i.e., and , respectively. The concentration values increases with time at each of the positions in both the media.
Figure 1

Concentration distribution pattern for exponentially decreasing velocity for time with average porosity of the medium.

Figure 1

Concentration distribution pattern for exponentially decreasing velocity for time with average porosity of the medium.

Figure 2 depicts the contaminant concentration pattern for the unsteady sinusoidal form of velocity with zero order production parameter. The effect of this parameter is predicted with respect to the aquifer (i.e., sand). The concentration level decreases with the increasing value of the zero-order production parameter, but the contaminant concentration increases with increasing time. The peak of the contaminant concentration is lower at the source due to the increasing zero-order production parameter. The effective parameters, like zero-order production, first-order decay, etc., in tropical regions are more significant for transport modeling for groundwater bodies.
Figure 2

Concentration distribution pattern for sinusoidal velocity for time with zero-order production in the aquifer.

Figure 2

Concentration distribution pattern for sinusoidal velocity for time with zero-order production in the aquifer.

Figure 3 compares the contaminant concentration patterns for exponentially decreasing type unsteady velocity expression in the time domain between aquifer (i.e., sand) and aquitard (i.e., clay). This is the case when the source has been removed from the aquifer. From this figure it is observed that after removing the source, some initial concentration exists at the origin. As a result, the effect of contaminant concentration increases with distance that attains the maximum peak in the case of the aquitard as compared to the aquifer and then decreases up to a harmless concentration with respect to the distance. The concentration pattern initially decreases with the increasing value of time and after covering some distance it takes the reverse pattern with respect to time in the aquifer and aquitard and ultimately goes to a minimum harmless concentration with respect to distance. From this figure it is observed that the concentration pattern is high in the aquitard as compared to the aquifer.
Figure 3

Concentration distribution for the exponential decreasing form of unsteady form of velocity for time with average porosity of the medium.

Figure 3

Concentration distribution for the exponential decreasing form of unsteady form of velocity for time with average porosity of the medium.

Figure 4 shows a comparison of the analytical result against the numerical one with averaging porosity of sand (0.37). The concentration distribution pattern follows its decreasing nature with respect to the distance in both the results. Initially, the concentration values are minimum in the numerical solution, but after covering some distance the numerical result attains slightly higher values in comparison to the analytical one. Both the results attain the minimum level at a particular point, after that it takes the reverse pattern in the case of exponentially decreasing form of the velocity pattern.
Figure 4

Concentration distribution pattern for the exponentially decreasing form of the velocity pattern for the sand medium.

Figure 4

Concentration distribution pattern for the exponentially decreasing form of the velocity pattern for the sand medium.

The contaminant concentration pattern for the asymptotic form and algebraic sigmoid form of velocity patterns are almost similar and are shown in Figures 5 and 6, respectively. The contaminant concentration increases slowly with respect to time and decreases rapidly with respect to distance. The effect of the Peclet number is shown in Figure 7. The Peclet number physically measures the relative magnitude of advection versus dispersion. The contaminant concentration increases with the increasing Peclet number and decreases with the decreasing Peclet number. Both patterns with different Peclet numbers reach minimum concentration with distance. The effect of the Peclet number on the solute concentration is observed with respect to the distance. For high Peclet number, the concentration level takes a minimum distance to reach its minimum concentration, whereas for low Peclet number it takes more distance for the same value. The above discussion is for time domain for the analytical solution given in Equation (18).
Figure 5

Concentration distribution pattern for asymptotic velocity with average porosity of the medium.

Figure 5

Concentration distribution pattern for asymptotic velocity with average porosity of the medium.

Figure 6

Concentration distribution pattern for algebraic sigmoid velocity with average porosity of the medium.

Figure 6

Concentration distribution pattern for algebraic sigmoid velocity with average porosity of the medium.

Figure 7

Concentration distribution pattern for algebraic sigmoid velocity with different values of Peclet number with the average porosity of the sand medium.

Figure 7

Concentration distribution pattern for algebraic sigmoid velocity with different values of Peclet number with the average porosity of the sand medium.

A similar type of concentration pattern was found for the case of the sinusoidal form of unsteady velocity pattern, as shown in Figure 8. The concentration is high in the aquitard as compared to the aquifer and ultimately goes to a minimum harmless concentration with respect to distance and time.
Figure 8

Concentration distribution for sinusoidal velocity for the time with average porosity of the medium.

Figure 8

Concentration distribution for sinusoidal velocity for the time with average porosity of the medium.

Figures 9 and 10 compare the concentration patterns for the asymptotic form and the algebraic sigmoid form of unsteady velocity patterns, respectively. From these figures it is observed that the concentration pattern attains a maximum peak value at a certain distance and thereafter it begins to decrease with respect to distance. The concentration pattern initially decreases with the increasing value of time and after covering some distance of approximately 0.005 for the case of aquifer (i.e., sand) and 0.012 for the case of aquitard (i.e., clay) the concentration takes a reverse pattern with respect to time and goes to a minimum harmless concentration. The concentration pattern takes on a lesser peak value in the case of asymptotic form and algebraic sigmoid form than for the exponentially decreasing and sinusoidal form of the velocity pattern for time domain .
Figure 9

Concentration distribution for asymptotic velocity for time with average porosity of the medium.

Figure 9

Concentration distribution for asymptotic velocity for time with average porosity of the medium.

Figure 10

Concentration distribution for the algebraic sigmoid velocity for time with average porosity of the medium.

Figure 10

Concentration distribution for the algebraic sigmoid velocity for time with average porosity of the medium.

Figure 11 exhibits contaminant concentration for different values of the Peclet number with the algebraic sigmoid form of unsteady velocity pattern for time domain . The solute concentration value increases with the increasing value of the Peclet number and the peak of the concentration is reduced with the increasing time.
Figure 11

Concentration distribution pattern for algebraic sigmoid velocity for time with different values of Peclet numbers with average porosity of the sand medium.

Figure 11

Concentration distribution pattern for algebraic sigmoid velocity for time with different values of Peclet numbers with average porosity of the sand medium.

Figure 12 predicts for the sinusoidal form of the velocity pattern for the averaging porosity of sand. It also follows the same type of nature as in Figure 4. The concentration value attains its minimum level in a short distance for the case of sinusoidal form of the velocity pattern in comparison to exponentially decreasing velocity pattern. From Figures 4 and 12, we observed that the decreasing nature of the concentration distribution patterns is faster for the sinusoidal form of velocity pattern (which shows the nature of groundwater contamination in the tropical region in which the fluctuation behavior of groundwater recharge is shown) in comparison to the exponentially decreasing form of velocity (Singh & Singh 2001; Jain et al. 2007) in which contaminant concentration follows a decreasing nature.
Figure 12

Concentration distribution pattern for the sinusoidal form of the velocity pattern for the sand medium.

Figure 12

Concentration distribution pattern for the sinusoidal form of the velocity pattern for the sand medium.

Figures 13 and 14 predict for the averaging porosity of the clay. After the rainy season the contamination fluctuation in groundwater is shown from the asymptotic and algebraic sigmoid forms of velocity patterns (Aral & Liao 1996; Singh & Kumari 2014). The concentration distribution pattern for the asymptotic form of the velocity pattern for the clay is shown in Figure 13. The concentration values are minimum for the numerical result in comparison to the analytical one and beyond some distance it takes a reverse pattern. Both the analytical and numerical concentration values attain their minimum level at a particular position. A similar type of nature of the contaminant distribution pattern was observed for the algebraic sigmoid form of the velocity pattern depicted in Figure 14.
Figure 13

Concentration distribution pattern for the asymptotic form of the velocity pattern for the clay medium.

Figure 13

Concentration distribution pattern for the asymptotic form of the velocity pattern for the clay medium.

Figure 14

Concentration distribution pattern for the algebraic sigmoid form of the velocity pattern for the clay medium.

Figure 14

Concentration distribution pattern for the algebraic sigmoid form of the velocity pattern for the clay medium.

In this paper, the RMS error was used to check the validity of numerical solution against the analytical one, as shown in Tables 1 and 2. The two parameters, and , play an important role to investigate the performance of the numerical solution. In the explicit finite difference scheme is restricted under the stability condition. Thus, in this present study, the accuracy was investigated by selecting different mesh sizes. The RMS error was investigated for , 0.05, and 0.07 for the particular time period 20 years in the time domain for the exponential decreasing and the asymptotic form of the velocity patterns, and 30 years with within the sinusoidal and algebraic sigmoid form of the velocity patterns in sand and clay media, respectively.

Table 1

The RMS values for the sand medium at particular 20 year in 0 < TTP

Distance Analytical result Numerical result
 
ΔZ = 0.02 ΔZ = 0.05 ΔZ = 0.07 
Case i: For exponential decreasing form of the velocity pattern 
 0.0019 0.5125 0.0119 0.0149 0.0169 
 0.0055 0.3188 0.0159 0.0248 0.0308 
 0.0093 0.1736 0.0198 0.0348 0.0448 
 0.0130 0.0837 0.0238 0.0448 0.0588 
 0.0167 0.0376 0.0278 0.0548 0.0727 
 RMS error  0.2719 0.2665 0.2634 
Case ii: For asymptotic form of the velocity pattern 
 0.0019 0.3613 0.0120 0.0150 0.0170 
 0.0055 0.0687 0.0160 0.0250 0.0310 
 0.0093 0.0130 0.0200 0.0350 0.0449 
 0.0130 0.0094 0.0240 0.0449 0.0589 
 0.0167 0.0091 0.0280 0.0549 0.0729 
 RMS error  0.1583 0.1585 0.1596 
Distance Analytical result Numerical result
 
ΔZ = 0.02 ΔZ = 0.05 ΔZ = 0.07 
Case i: For exponential decreasing form of the velocity pattern 
 0.0019 0.5125 0.0119 0.0149 0.0169 
 0.0055 0.3188 0.0159 0.0248 0.0308 
 0.0093 0.1736 0.0198 0.0348 0.0448 
 0.0130 0.0837 0.0238 0.0448 0.0588 
 0.0167 0.0376 0.0278 0.0548 0.0727 
 RMS error  0.2719 0.2665 0.2634 
Case ii: For asymptotic form of the velocity pattern 
 0.0019 0.3613 0.0120 0.0150 0.0170 
 0.0055 0.0687 0.0160 0.0250 0.0310 
 0.0093 0.0130 0.0200 0.0350 0.0449 
 0.0130 0.0094 0.0240 0.0449 0.0589 
 0.0167 0.0091 0.0280 0.0549 0.0729 
 RMS error  0.1583 0.1585 0.1596 
Table 2

RMS values for the clay medium with their averaging porosity at particular 30 years for T > Tp

Distance Analytical result Numerical result
 
ΔZ = 0.02 ΔZ = 0.05 ΔZ = 0.07 
Case i: For the sinusoidal form of the velocity pattern 
 0.0028 0.1094 0.0118 0.0148 0.0168 
 0.0083 0.2425 0.0158 0.0248 0.0307 
 0.0138 0.2524 0.0198 0.0347 0.0447 
 0.0193 0.1869 0.0238 0.0447 0.0586 
 0.0248 0.1135 0.0278 0.0546 0.0725 
 RMS error  0.1726 0.1596 0.1514 
Case ii: For the algebraic sigmoid form of the velocity pattern 
 0.0028 0.1001 0.0119 0.0149 0.0169 
 0.0083 0.1608 0.0159 0.0249 0.0309 
 0.0138 0.0868 0.0199 0.0349 0.0449 
 0.0193 0.0301 0.0239 0.0449 0.0589 
 0.0248 0.0126 0.0279 0.0549 0.0728 
 RMS error  0.0818 0.0780 0.0774 
Distance Analytical result Numerical result
 
ΔZ = 0.02 ΔZ = 0.05 ΔZ = 0.07 
Case i: For the sinusoidal form of the velocity pattern 
 0.0028 0.1094 0.0118 0.0148 0.0168 
 0.0083 0.2425 0.0158 0.0248 0.0307 
 0.0138 0.2524 0.0198 0.0347 0.0447 
 0.0193 0.1869 0.0238 0.0447 0.0586 
 0.0248 0.1135 0.0278 0.0546 0.0725 
 RMS error  0.1726 0.1596 0.1514 
Case ii: For the algebraic sigmoid form of the velocity pattern 
 0.0028 0.1001 0.0119 0.0149 0.0169 
 0.0083 0.1608 0.0159 0.0249 0.0309 
 0.0138 0.0868 0.0199 0.0349 0.0449 
 0.0193 0.0301 0.0239 0.0449 0.0589 
 0.0248 0.0126 0.0279 0.0549 0.0728 
 RMS error  0.0818 0.0780 0.0774 

In both the tables, was fixed. Tables 1 and 2 were tabulated for the RMS error in the aquifer (i.e., sand) and aquitard (i.e., clay) for four different types of the velocity patterns, respectively. The RMS error decreases with the increasing grid space for the exponential decreasing form of the velocity pattern for the sand medium, which was observed from Case i of Table 1. In the asymptotic form of the velocity pattern the RMS error increases with respect to the increasing grid space in Case ii of Table 1. The RMS error attains its minimum value with the increasing mesh size for Cases i and ii in Table 2 for the clay medium. The RMS error was evaluated for the accuracy of solution for the sinusoidal form of the velocity pattern tabulated in Case i in Table 2, and for the algebraic sigmoid form of the velocity pattern tabulated in Case ii of Table 2. In both the velocity patterns the result is more accurate for the maximum value of the mesh size, except in the case of the asymptotic form of the velocity pattern where the result is more accurate in the case of minimum value of the mesh size.

VALIDATION OF THE MODEL WITH EXISTING SOLUTION OF LIU ET AL. (1998) 

One-dimensional ADE in multilayer porous media was solved analytically using generalized integral transform technique by Liu et al. (1998), where the analytical solution was derived under arbitrary initial and boundary conditions. In this present paper, the authors have shown the validation of the model equation with the existing research work done by Liu et al. (1998). The analytical solution obtained in Equations (18) and (19) was computed for the same set of input data, except some parameters which had an effect on the solute transport modeling, which have been compared with the input values taken by Liu et al. (1998) as follows: 
formula
 
formula
The following inputs are taken for validation purposes in this present paper. The concentration distribution pattern for Input (i) and (ii) are predicted for the different geological formations with their averaging porosity and shown in Figures 15 and 16 for different time domains. The concentration distribution pattern for a particular time, 10 years, for sand and clay medium are predicted and shown in Figure 15, and for 25 years is also predicted and shown in Figure 16. The concentration values in each of the positions in Input (i) are higher as compared to Input (ii) in both the media as observed in Figure 15. The clay medium attains maximum concentration values as compared to the sand medium, but both patterns ultimately follow minimum concentration values with their respective distance. After removing the source, contaminant concentration increases with distance up to a certain limit, then decreases with respect to the distance in both the medium for Input (i) and (ii), as seen in Figure 16. The concentration pattern attains its maximum peak in the case of the aquitard (i.e., clay) as compared to the aquifer (i.e., sand) in both Input (i) and (ii). Initially, minimum concentration values are attained for Input (i) as compared with Input (ii), but after covering a certain distance the reverse pattern is observed. At the end of the position, minimum concentration values are attained by Input (ii), which have been taken in the present paper for illustration.
Figure 15

Concentration distribution pattern with different input for time with average porosity of the medium.

Figure 15

Concentration distribution pattern with different input for time with average porosity of the medium.

Figure 16

Concentration distribution with the different input for time with average porosity of the medium.

Figure 16

Concentration distribution with the different input for time with average porosity of the medium.

SUMMARY AND CONCLUSIONS

Employing the concept of linear isotherms, analytical solutions for the ADE with respect to the solid and liquid phase are derived. The mixed type boundary condition is employed at the source in the splitting time domain. The contaminant concentration patterns for different types of velocity patterns are evaluated. The following conclusions can be drawn from this study:

  1. The impact of contaminant concentration for linear isotherms with the distribution coefficient is significantly observed in the splitting time domain for different velocity patterns, such as exponentially decreasing, sinusoidally varying, algebraic, and asymptotic forms.

  2. The contaminant concentration values depend upon the decreasing or increasing values of the zero-order production term and first-order decay rate coefficient.

  3. The contaminant concentration distribution behavior is predicted for different geological formations in two time domains, i.e., and .

  4. Comparison of the analytical result with the numerical result is taken into account. Accuracy of the solution is significantly observed by using RMS error. Truncation error of various parameters is also explored, which causes the inconsistency among analytical and numerical results.

  5. The validation of the model is made with the result of an existing solution given by Liu et al. (1998) and the same trend for contaminant concentration was found.

ACKNOWLEDGEMENTS

The authors are grateful to the Indian School of Mines, Dhanbad for providing financial support to PhD candidate under the ISMJRF scheme. The authors are also grateful to the editor and reviewers for the comments which helped improve the quality of the paper.

REFERENCES

REFERENCES
Ataie-Ashtiani
B.
Lockington
D. A.
Volker
R. E.
1996
Numerical correction for finite-difference solution of the advection-dispersion equation with reaction
.
J. Contam. Hydrol.
23
,
149
156
.
Ataie-Ashtiani
B.
Lockington
D. A.
Volker
R. E.
1999
Truncation errors in finite difference models for solute transport equation with first-order reaction B
.
J. Contam. Hydrol.
35
,
409
428
.
Balla
K.
Keri
G.
Rapcsak
T.
2002
Pollution of underground water – a computational case study using a transport model
.
J. Hydroinform.
4
(
4
),
255
263
.
Batu
V.
2006
Applied Flow and Solute Transport Modeling in Aquifers: Fundamental Principles and Analytical and Numerical Methods
.
CRC Press
,
Boca Raton, FL
,
USA
.
Bear
J.
Verruijt
A.
1987
Modeling Groundwater Flow and Pollution
.
Reidel
,
Dordrecht
,
The Netherlands
.
Bedient
P. B.
Rifai
H. S.
Newell
C. J.
1999
Ground Water Contamination: Transport and Remediation
,
2nd edn
.
Prentice-Hall PTR
,
Upper Saddle River, NJ
,
USA
.
Carnahan
B.
Luther
H. A.
Wilkes
J. O.
1969
Applied Numerical Methods
.
Wiley
,
New York
,
USA
, pp.
429
530
.
Chrysikopolous
C. V.
Roberts
P. V.
Kitanidis
P. K.
1990
One-dimensional solute transport in porous media with partial well-to-well recirulation: application to field experiment
.
Water Resour. Res.
26
(
6
),
1189
1195
.
Crank
J.
1975
The Mathematics of Diffusion
,
2nd edn
.
Oxford University Press
,
London
,
UK
.
Freeze
R. A.
Cherry
J. A.
1979
Groundwater
.
Prentice-Hall
,
Upper Saddle River, NJ
,
USA
.
Fried
J. J.
Combarnous
M. A.
1971
Dispersion in porous media
.
Adv Hydrosci.
7
,
169
281
.
Gao
G.
Zhan
H.
Feng
S.
Fu
B.
Huang
G.
2012
A mobile-immobile model with an asymptotic scale-dependent dispersion function
.
J. Hydrol.
424–425
,
172
183
.
Hill
D. J.
Minsker
B. S.
Valocchi
A. J.
Babovic
V.
Keijzer
M.
2007
Upscaling models of solute transport in porous media through genetic programming
.
J. Hydroinform.
9
(
4
),
251
266
.
Hunt
B.
1978
Dispersion sources in uniform groundwater flow
.
J. Hydraul. Div.
HY1
,
74
85
.
Jain
S. K.
Agarwal
P. K.
Singh
V. P.
2007
Hydrology and Water Resources of India
.
Springer
,
Dordrecht
,
The Netherlands
, p.
87
.
Mickley
H. S.
Sherwood
T. K.
Reed
C. E.
1957
Applied Mathematics in Chemical Engineering.
McGraw-Hill
,
New York
,
USA
.
Neville
C. J.
Ibaraki
M.
Sudicky
E. A.
2000
Solute transport with multiprocess nonequilibrium: a semi analytical solution approach
.
J. Cont. Hydrol.
44
(
2
),
141
159
.
Notodarmojo
S.
Ho
G. E.
Scott
W. D.
Davis
G. B.
1991
Modelling phosphorus transport in soils and groundwater with two-consecutive reactions
.
Water Res.
25
(
10
),
1205
1216
.
Rausch
R.
Schafer
W.
Therrien
R.
Wagner
C.
2005
Solute Transport Modelling: An Introduction to Models and Solution Strategies
.
Gebruder Borntrager
,
Berlin
,
Germany
.
Richard
A. S.
Jr
1988
Small grid testing of finite difference transport schemes
.
J. Hydraul. Eng.
114
(
3
),
329
336
.
Singh
M. K.
Kumari
P.
2014
Contaminant concentration prediction along unsteady groundwater flow
.
Modell. Simul. Diff. Process.
XII
,
257
276
.
Singh
P.
Singh
V. P.
2001
Snow and Glacier Hydrology
.
Kluwer Academic Publishers
,
Amsterdam
,
The Netherlands
, p.
78
.
Smith
G. D.
1978
Numerical Solution of Partial Differential Equations: Finite Difference Methods
,
2nd edn
.
Oxford University Press
,
Oxford
,
UK
.
Zheng
C.
Bennett
G. D.
2002
Applied Contaminant Transport Modeling
,
2nd edn
.
Wiley
,
New York
,
USA
.

Supplementary data