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.

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.

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

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:

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.

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:
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:whereThe 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.
The analytical solution, obtained by Equations (18) and (19), was computed for the following data (Singh & Kumari 2014):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: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.

Close modal
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.

Close modal
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.

Close modal
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.
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 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.

Close modal
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.

Close modal
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 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.

Close modal
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.

Close modal
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.
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.

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

DistanceAnalytical resultNumerical 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
DistanceAnalytical resultNumerical 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

DistanceAnalytical resultNumerical 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
DistanceAnalytical resultNumerical 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.

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

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.

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.

Aral
M. M.
Liao
B.
1996
.
J. Hydrol. Eng.
1
(
1
),
20
32
.
Ataie-Ashtiani
B.
Hosseini
S. A.
2005
.
Environ. Modell. Softw.
20
,
817
826
.
Ataie-Ashtiani
B.
Lockington
D. A.
Volker
R. E.
1996
.
J. Contam. Hydrol.
23
,
149
156
.
Ataie-Ashtiani
B.
Lockington
D. A.
Volker
R. E.
1999
.
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
,
,
USA
.
Carnahan
B.
Luther
H. A.
Wilkes
J. O.
1969
Applied Numerical Methods
.
Wiley
,
New York
,
USA
, pp.
429
530
.
Chai
T.
Draxler
R. R.
2014
.
Geosci. Model. Dev.
7
,
1247
1250
.
Chaudhari
N. M.
1971
.
Soc. Petrol. Eng. J.
11
,
277
284
.
Chen
J. S.
2007
.
30
(
3
),
430
438
.
Chrysikopolous
C. V.
Roberts
P. V.
Kitanidis
P. K.
1990
.
Water Resour. Res.
26
(
6
),
1189
1195
.
Crank
J.
1975
The Mathematics of Diffusion
, 2nd edn.
Oxford University Press
,
London
,
UK
.
Dudley
L. M.
McLean
J. E.
Furat
T. H.
Jurinak
J. J.
1991
.
Soil Sci.
151
,
121
135
.
Freeze
R. A.
Cherry
J. A.
1979
Groundwater
.
Prentice-Hall
,
,
USA
.
Fried
J. J.
Combarnous
M. A.
1971
.
7
,
169
281
.
Gao
G.
Zhan
H.
Feng
S.
Fu
B.
Huang
G.
2012
.
J. Hydrol.
424–425
,
172
183
.
Hantush
M. M.
Marino
M. A.
1998
.
J. Hydrol. Eng.
3
(
4
),
241
247
.
Hill
D. J.
Minsker
B. S.
Valocchi
A. J.
Babovic
V.
Keijzer
M.
2007
.
J. Hydroinform.
9
(
4
),
251
266
.
Hogarth
W. L.
Noye
B. J.
Stagnitti
J.
Parlange
J. Y.
Bolt
G.
1990
.
Comput. Math. Appl.
20
(
11
),
67
82
.
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
.
Kumar
S. G.
Sekhar
M.
Misra
D.
2008
.
J. Hydrol. Eng.
13
(
4
),
250
257
.
Kwok
Y. K.
1992
.
Comput. Math. Appl.
23
(
12
),
3
11
.
Lantz
R. B.
1971
.
Soc. Petrol. Eng. J.
11
,
315
320
.
Liu
C.
Ball
W. P.
Ellis
J. H.
1998
.
Trans. Porous Med.
30
,
25
43
.
Logan
J. D.
1996
.
J. Hydrol.
184
(
3–4
),
261
276
.
Lowry
T.
Li
S. G.
2005
.
28
,
117
133
.
Meerschaer
M. M.
C.
2004
.
J. Comput. Appl. Math.
172
,
65
77
.
Mickley
H. S.
Sherwood
T. K.
Reed
C. E.
1957
Applied Mathematics in Chemical Engineering.
McGraw-Hill
,
New York
,
USA
.
Mohebbi
A.
M.
2013
.
Numer. Algor.
63
,
431
452
.
Moldrup
P.
Kruse
C. W.
Yamaguchi
T.
Rolston
D. E.
1996
.
Soil Sci.
161
,
347
354
.
Molenkamp
C. R.
1968
.
J. Appl. Meteor.
7
,
160
167
.
Neville
C. J.
Ibaraki
M.
Sudicky
E. A.
2000
.
J. Cont. Hydrol.
44
(
2
),
141
159
.
Notodarmojo
S.
Ho
G. E.
Scott
W. D.
Davis
G. B.
1991
.
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
.
Rezaei
A.
Zhan
H.
Zare
M.
2013
.
J. Cont. Hydrol.
152
,
117
136
.
Richard
A. S.
Jr
1988
.
J. Hydraul. Eng.
114
(
3
),
329
336
.
Roberts
D. L.
Selim
M. S.
1984
.
Int. J. Numer. Methods Eng.
20
,
817
844
.
Singh
M. K.
Das
P.
2015
.
J. Hydrol.
520
,
289
299
.
Singh
M. K.
Kumari
P.
2014
.
Modell. Simul. Diff. Process.
XII
,
257
276
.
Singh
P.
Singh
V. P.
2001
Snow and Glacier Hydrology
.
,
Amsterdam
,
The Netherlands
, p.
78
.
Singh
M. K.
Singh
V. P.
Singh
P.
Shukla
D.
2009
.
J. Eng. Mech.
135
(
9
),
1015
1021
.
Smedt
F. D.
2006
.
J. Hydrol.
330
,
672
680
.
Smith
G. D.
1978
Numerical Solution of Partial Differential Equations: Finite Difference Methods
, 2nd edn.
Oxford University Press
,
Oxford
,
UK
.
Tkalich
P.
2006
.
J. Hydroinform.
8
(
3
),
149
164
.
Towler
B. F.
Yang
R. Y. K.
1979
.
Int. J. Numer. Methods Eng.
14
,
1021
1035
.
van Genuchten
M. T.
1981
.
J. Hydrol.
49
,
213
233
.
van Genuchten
M. T.
Gray
W. G.
1978
.
Int. J. Numer. Methods Eng.
12
,
387
404
.
Wang
P. P.
Zheng
C.
Gorelick
S. M.
2005
.
28
,
33
42
.
Yeh
H. D.
Yeh
G. T.
2007
.
J. Environ. Eng.
133
(
11
),
1032
1041
.
Zhan
H.
Zhang
W.
Guanhua
H.
Dongmin
S.
2009
.
J. Contam. Hydrol.
107
,
162
174
.
Zheng
C.
Bennett
G. D.
2002
Applied Contaminant Transport Modeling
, 2nd edn.
Wiley
,
New York
,
USA
.
Zoppou
C.
Knight
J. H.
1994
.
Water Resour. Res.
30
(
11
),
3233
3235
.