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

*n*is the porosity of geological formation.

*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: A new time variable (Crank 1975): is introduced and Equation (8) can be written as: The following non-dimensional variables: are used to reduce the number of parameters of Equation (10) and one can then write: where The initial and boundary conditions given in Equation (3), Equation (5), and Equation (6) can be written in non-dimensional form, respectively, as follows: Now, the following transformation: 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): where

## NUMERICAL SOLUTION

*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: Equation (12) together with the initial and boundary condition may be written as:

*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: where superscript

*j*refers to time, subscript

*i*refers to space, is the time increment, and is the space increment.

## STABILITY CONDITION

*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: where

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*.

## 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.

*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:

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:

## ACCURACY OF THE SOLUTION

*N*is the number of data which were used to evaluate the accuracy of the solution.

## RESULTS AND DISCUSSION

*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:

*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.

*et al.*2007) in which contaminant concentration follows a decreasing nature.

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.

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 |

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)

*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: 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.

## 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:

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.

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

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

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.

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.