A Two-Dimensional Mathematical Model of Transportation of Contaminants in Unsaturated Porous Media with Uniform Flow

A Two-Dimensional Mathematical Model of Transportation of Contaminants in Unsaturated Porous Media with Uniform Flow

Shobhankumar D M

Associate Professor, Department of Mathematics, Maharanis Science College for Women, Palace Road, Bengaluru-560001, Karnataka, India

Machhindranath M Dhane

Associate Professor, Department of Mathematics, GFGC, Yelahanka, Bengaluru- 560064, Karnataka, India

Ramesh T

Professor, Department of Mathematics Presidency University, Bengaluru-560064, Karnataka, India

Abstract

The main objective of this paper is to develop a mathematical model for two-dimensional flow of pollutants in unsaturated porous media with uniform flow using isotopes. The advection dispersion equation (ADE) is used to develop the model with suitable initial and boundary conditions to find the solution and compare with the other solution. The displacement experiment gives the dispersion coefficient D. The magnitude of dispersion depends on the particle size distribution and flow parameters. The average pore velocity has been considered only for the vertically downward direction. Here we have considered DL is the longitudinal dispersion coefficient and DT is the transverse dispersion coefficients which are perpendicular to flow directions of the fluid. The value of the DT is more complicated to obtain than the value of DL, because of the concentration distribution are required to be measured in vertical to the flow direction. Some studies have been performed to get the transverse dispersion coefficient on velocity of unpolluted water and porous media with distribution of various particle sizes to find out the values of longitudinal and transverse dispersion coefficient.

Keywords: Contaminants, porous media, Diffusion, Dispersion, Advection Adsorption, Uniform flow, Integral transforms.

1 INTRODUCTION

Mathematical models are being used extensively in groundwater studies. Groundwater modeling can be classified in to two model, they are ground water flow and solute transport models. The Solute transport models are applied in connection with groundwater quality problems. The solute transport models are extended in few cases with chemical sub models for description of the fate of non-conservative polluting species while in some case may be sufficient only to study the flow of the pollutant in the sub surface. This can be analyzed by considering two and three-dimensional model to prevent the groundwater pollution.

Previous works closely related to the work are also conducted by Yeh (1981), Domenico and Robbins (1984), Domenico (1987), Batu (1989, 1993), Leij et al (2000), and Park and Zhan (2001). Yeh (1981) provided a general framework for using Green's functions to solve transport equations without giving too many details. Domenico and Robbins (1984) and Demenico (1987) considered finite sources as boundary conditions when solving the advection-dispersion equation. They have not included the effect from the upper and lower boundaries of an aquifer. Batu (1989, 1993) provided a two-dimensional analytical solute transport model in a bounded aquifer by using the same source as the aquifer thickness along the z-axis and include the contaminant source as a boundary condition. The general solutions were derived with the help of Fourier analysis and Laplace domain. The contaminant concentration in the near field is found to be sensitive to the source geometry and anisotropy of the dispersion coefficients. The contaminant concentration in the far field is found to be less sensitive to the source geometry. All the above solutions have some restrictions. Therefore, our objective is to provide an analytical solution for two-dimensional transport problem which is applicable for both longitudinal and transverse dispersion.

2 MATHEMATICAL MODEL

Hence, a comprehensive groundwater modeling must get the models for saturated and unsaturated porous media is integrated in the groundwater/subsurface water models. The affects of hydrodynamic processes and physical aspects that influences the mass balance on groundwater includes, advection, dispersion, diffusion, transportation, adsorption, and transformation in chemical reactions.

The infiltration is assumed that the solute transport is obtained with ADE. Let us consider the two-dimensional dispersion with one-dimensional steady state flow is of the form

C

t DL

2C

z2

• w C

  z

  DT

  2C

  x2

• C

  (1)

  where V is the pure water velocity [LT-1], t is time dependent [T], C is the concentration solution [ML-3], DL and DT are longitudinal and transversion [L2T-1] coefficient, z and x are the particular points along the Cartesian coordinate axes which is perpendicular and parallel to the direction of ground water flow [L].

  For groundwater flow, a linear relation b/w the dispersion coefficient and seepage velocity is generally adopted. If the aquifer is isotropic dispersion coefficient can be characterized by a transverse and longitudinal coefficient (Fig.1).

  Figure 1. Physical layout of two-dimensional transport of contaminants in unsaturated porous media.

  For this problem the initial and boundary conditions are as follows:

  C x, z, 0 f x, 0 z , x (2a)

  C

  z

  z

  0,

• x , t 0

  (2b)

  C 0, x, t g x CL , x 0, t 0

  C C

  C 0, x, t g x L R , x 0, t 0

  (2c)

  2

  C 0, x, t g x CR , x 0, t 0

  C

  x

  0, 0 z , t 0

  (2d)

  x

  The initial condition (2a) is subjected to find the solution of equation (1) which can be obtained with integral transformation as, Fourier transformation for x variable and Laplace transformation for t variable. The solution is the combination of semi-infinite plane for steady state conditions considering one-dimensional advection- dispersion model and one-dimensional advection and two-dimensional dispersion model.

  The solution for the problem is solved first by using the well-known Laplace transforms applications for C and

  C w.r.t to t the utilization of the initial condition (2a), for the transform equation (1) is transformed in the form

  t

  LC x, z,t Exp pt C x, z,t dt Cx, z, p

  0

  L C pC(x, z, p) C(x, z, 0)

  t

  Applying the Laplace transforms to Eqns. (1) and (2), it transforms in the form

  (6.3)

  (3)

  pC(x, z, p) C(x, z, 0) DL

  2C

  z2

  v C

  z

  DT

  2C

  x2

• C

  (4)

  pC(x, z, p)

• f DL

2C

z2 DT

2C

x2 u1

C v

x 1

C C

z

(5)

C

z

and boundary conditions reduces to

0

z

(6a)

C0, x, p g

C

x

p

0

(6b)

(6c)

x

For the infinite x-domain the Fourier transform is applied and for C

transform is given by

and

2C

t 2

the Fourier

F C z, x, p 1

2

Exp i xC z, x, p dx C z, , p

(7)

and

C 2

2

F

x C

(8)

The Fourier transforms of Eqns. (5) and (6) becomes

D

D/p>

d 2C 2 dC dC

L dz2

T C u1 dx

v

1 dz

C pC (x, z, p) f 0

d 2C

D v

dC u

dC (D 2

(9)

L dz2

1 dz

1 dx T

p )C f 0

and the reduced boundary conditions are

dC 0

dz z

p

C0, x. p g

dC

dx

(10)

0

x

The Fourier transformation of f and g are and .

The transport equation (9) which is an ordinary differential equation is subjected to the boundary conditions, then equation (10) is of the form

where

C z, , p AExp(R1z) B Exp(R2 z)

f

T

D 2 p

(11)

1

L

v1		v 2	4D	D 2 p T
2D v1		v 2	4D	2D D 2	p

R1

L L

(12a)

1

L T

R2 2D 2D

(12b)

L L

It follows from the Eqn. (10) that A 0.Therefore we will use R instead of R2 . Applying the

inner boundary conditions of Eqn.(10), we get

B g f

(13)

T

p D 2 p

Substituting Eqn (13) in (11), we get

C z,, p g f Exp(Rz) g f

(14)

p D 2 p

p D 2 p

T T

Inverse Laplace transform has been applied for the given domain z, x,t in order to get the solution. This transport

problem can also be solved using numerical technique but to solve using analytical method. The RHS of Eqn.(14) has divided into sum of three terms, the first integral term reduces to

C z, , p L1 g Exp(Rz)

p

Exp

v z 1

z v 2

12

C z, , p g Exp 1

L1

1 D 2 p

(15)

2D p

D 4D T

L

L L

From inverse Laplace transformation table [van Genuchten and Alves (1982)], by the application of convolution theorem it is preferable to solve the inverse Laplace transform than inverse Fourier transform. Let f and g be

two continuous functions in the interval 0, , then convolution of f and g is denoted f g and defined by the integrals (G Zill and R cullin (2000))

The convolution

t

f g f g t d

0

f g is a f(t). The limitation of convolution is to find the Laplace transformation of two

functions without evaluate the integral.

Let h(t) and k(t) be two functions and its Laplace transformations are these integrals are

h( p) & k( p) then the convolution of

L1 h p.k p h k

t t

L1 h pk p h k t d h t k d

0 0

Here is a variable. The difference in two functions in the Eqn.(15) is of the form

(16)

h p 1

p

z v 2

12

(17)

k t Exp 1 D 2 p

T

DL

4D

L

The values of h(t) and k(t) are calculate by using shift property and Laplace transformation is given by

h t 1

z

v 2

2 z2

(18)

k t Exp 1 DT t

t 4 DL

4DL

4DLt

Substitute (18) in (16) and subsequently into Eqn.(15) reduces to

gz t 3

(z v )2

C z,, p

2 Exp 1 Exp (D 2 ) d

(19)

0

T

4 DL

4DL

By using Laplace transforms table, the Eqn.(19) may be transformed in the form of complimentary error function. For the integral of the error function is given by

erf z

2 e2 d

z

0

The Complementary error function of

f (z)

is given by

erfc z 1 erf z

2

e2 d

z

Now applying the complimentary error function to the expression (19), becomes

g v z

z

v 2

12

C z, , t Exp 1 Exp 1 D 2

1 2 2D D 4D T

L L L

z v 2t z v 2

12

erfc 1 D 2t t Exp 1 D 2

4D T 4D T

D

4D T

L L

z

L L

v 2t

2

erfc

1 DT t t

4DLT 4DL (20)

Inverse Fourier transform of the Eqn.(20) is difficult than for Eqn.(14). The Fourier transforms variable displays in the argument of both error and exponential functions. Applying the inverse Fourier transform to the second term C2(z,, p) of RHS of the Eqn. (14), we have

C (z,,t) L1 f Exp Rz

(21)

T

2 D 2 p

C (z,, t) f Exp D 2 t

z v1t

v z

Exp 1 erfc

z v1t

(22)

2 2 T

erfc

D

L

4DLt

4DLt

Now applying inverse Laplace transform of the 3rdterm of the RHS of Eqn.(14) may be written the form [Oberhettinger and Baddi(1973)]

C (z,,t) L1 f f Exp D 2 t

(23)

3 D 2 p T

T

The required Equation for C (z,,t) can be written as

C(z,,t) C1(z,,t) C2 (z,,t) C3 (z,,t)

(24)

The final step of the solution is in the other form, then the application of the Fourier inversion

in terms of the Eqn. (24). The inverse Fourier transform of C (z,,t) is

C(z, x, t) F

1 C(z, x, t)

2

1

Exp i x C z, , t d

(25)

It is left with the Fourier inversion of the 1stterm of the Eqn. (24) may be written as

1

gz t 3

(z v )2

F 1

C1 z, x, t

Exp i x

2

4 D

2 Exp

1

4D

Exp

D 2

d d

T

L 0 L

t

z 3

(z v )2

F 1 C z, x, t 2 Exp 1

0

1 4 D

L

4DL

(26)

1 Exp i xgExp D 2 d d

2

T

The convolution h k is to determine convolution of product of two functions h and k

is

1

1

2

2

h k Exp i x d h k

h x k d

(27)

1

Where is the fake variable, after verifying the last term for C (z, x,t) in equation (26), using

the convolution integral for the other two functions are

h g

k Exp D 2

(28)

T

The Fourier transform of the single step function we have h(x) is equal to g(x) , but not

determining g and to find the inverse transform of k () , we have

2

F 1 Exp DT

1 Exp

x2

(29)

2 DT

The resultant h(x) and

k(x)

4 DT

are

h(x) g(x) CL , : x 0

(C C )

h(x) g(x) L R ,: x 0

(30)

and

h(x) g(x) CR

2

p>

, : x 0

1

1 x2

(31)

2 DT

k (x)

k Exp i x d

Exp 4 D

2

T

From Eqn. (27) is of the form

1 g (x )2

2(DT )

h k

Exp 4(D

) d

(32)

2

The above expression can be evaluated by

4(DT )

T

x

The conditions of Eqn.(30) and the complimentary error function [Abramowiz and Stegan(1970 and Crank(1975))] converts to

h k CL erfc x CR erfc x

(33)

2 4(D ) 2 4(D )

T T

The result is then substituted in (30) to get expression for C1(z, x, t) , is

t z t 3 (z v )2

C1 z, x, t h k 2 Exp 1

0 4 DL 0 4DL

CL erfc x CR erfc x

d

(34)

2 4(D ) 2

4(D

)

T

T

The inverse of the second term of the Eqn. (26), C2 (z, x,t) is

erfc

C (z,, t) 1

z t

Exp z erfc

z t

F 1 f Exp D

t

(35)

4DLt

2 2

D

T

L

4DLt

Applying the convolution theorem, the inverse Fourier transform in Eqn. (35),

h f

k Exp D 2

(36)

T

Let us assume that the initial concentration Ci is constant and Eqn. (29) is utilized to find k(x)

may be written as

h x Ci

1

x2

k x

Exp 4(D

)t

2(DT )t

T (37)

The inverse transformation of equation is carried out by using the properties or error function

F 1

f Exp

D 2 t

h k Ci

4 DT t

4 DT t

Exp

(x )2

d

T

(38)

Substitution of Eqn.(31) into Eqn.(28), then C2 (z, x, t) reduces to

C (z, x,t) Ci

z t

z

z t

(39)

4DLt

4DLt

2 2 erfc Exp D erfc

L

The inverse Fourier transform of the last term of Eqn. (24) then C3(z, x,t) is evaluated in Eqn.

(38) and we get

C (z, x, t) F 1 f Exp D 2 t C

3 T

C

i

x

4 (DT )t

4 (DT )t

C3 z, x, t i Exp

d

(40)

Substitution of C1(z, x, t) , C2 (z, x, t) and C3(z, x,t) in equation (24) and it reduces to

t 3

C z, x, t

z 2 CL erfc

x CR erfc x

4 DL

4(D 1 n )

T d

n

4(D 1 n )

T d

n

2 2

0

(z )2 C z t z z t

4DLt

4DLt

Exp d i erfc Exp erfc C

(41)

4DL

2

DL

Various numerical methods can be used to evaluate Eqn.(41), among the suitable method of solution Gauss Chebyshev quadrature. The transport of solute is linearly adsorbed using the rough layer that can be solved by splitting variables DL, DT and v with the factor of retardation. The solution may be different for initial and inlet conditions, like f and g.

Figure 2. BTC curve with sand for Q = 100mL/d

Figure 3. BTC curve with sand for Q = 50mL/d

Figure 4. BTC curve with sand for Q = 20mL/d

Figure 5. BTC curve with sand for Q = 2mL/d

Figure 6. BTC curve with sand and clay (3:1) for Q = 20mL/d

Figure 7. BTC curve with sand and clay (1:1) for Q = 2mL/d

Figure 8. BTC curve sand and clay (1:3) for Q = 2mL/d

Figure 9. BTC curve with clay for Q = 2mL/d

3 RESULTS AND CONCLUSIONS

The solution provided in Eqn. (41), wherein the mathematical integration of the first term is verified through comparing with various mathematical solutions for particular values of CL, CR, and Ci. First, the solution of stable condition Eqn.(1) is used below the idea that longitudinal dispersion DL may be neglected. This solution could be obtained through using the Fourier transformation of x may be written as Harleman and Rumer (1963).

For the affirmation we assume that a permeable medium with the arbitrary transportation parameters are DL=25cm2/d, DT=5cm2/d, and v=50cm/d. These concentration parameters are expressed as dimensionless quantity C/Co with Co =1.

The primary concentrations are CL=1, CR=0, and C1=0. Figures 2 to 9 shows that the solution for several times according to equation (41). The incursion of the solute into the medium at x > 0 and the consequent destruction of the solute in both directions may be observed clearly.

Figures 2 and 9 represent the Break-Through-Curves for C/C0 verses time and intensity for specific depths z and x and t and x. It is appearing to be awareness location will boom within the starting and reaches constant nation price for consistent z and x but decreases with a growth within the radioactive decay coefficient 1. A growth in 1 will make the solute awareness lower as is evident from the bodily floor.

It should be noted that equation (41) can suitably use to determine DT and DL simultaneously. The awareness of these transport parameters is necessary for the prediction of two-dimensional solute transport.

In this paper we are capable of see that an analytical solution is received for the two-dimensional ADE for a semi- infinite region (half of plane) with one-dimensional waft using integral transforms. Even though the solution includes a fundamental expression, this essential can also want to effectively be evaluated using Gauss-Chebyshev quadrature.

4 REFERENCES

1. Sudheendra S.R., 2011. A solution of the differential equation of dependent dispersion along uniform and non-uniform flow with variable coefficients in a finite domain, Int. J. of Mathematical Analysis, Vol.3, No. 2, 89-105.

2. Sudheendra S.R. 2012. An analytical solution of one-dimensional advection-diffusion equation in a porous media in presence of radioactive decay, Global Journal of Pure and Applied Mathematics, Vol.8, No. 2, 113-124.

3. A. Verruijt (1990): Modelling Groundwater Flow and Pollution. D Reidel Publishing Co., Tokyo.

4. Abramowitz, M., and I. A. Stegun (1970): Handbook of Mathematical Functions, Appl. Math. Ser. Vol. 55, National Bureau of Standards, Washington, D.C.

5. Al-Niami A N S and Rushton K R 1977 Analysis of flow against dispersion in porous media; J. Hydrol. 33 8797.

6. Aral M M and Liao B 1996 Analytical solutions for two-dimensional transport equation with time-dependent dispersion coefficients; J. Hydrol. Engg. 1(1) 2032.

7. Banks R B and Ali J 1964 Dispersion and adsorption in porous media flow; J. Hydraul. Div. 90 13-31.

8. Batu,V.,1993. A genaralized two-dimensional analytical solute transport model in bounded media for flux-type multiple sources. Water Resour.Res.,29,2881-2892.

9. Bear J 1972 Dynamics of fluids in porous media (New York: Amr. Elsev. Co.).

10. Bruce, J.C.; and Street, R.L. (1966). Studies of free surface flow and twodimensional dispersion in porous media. Report No. 63, Civil Eng. Dept. Stanford Univivrsiy.

11. Brutsaert F. W.(1971): A functional iteration technique for solving the Richards equation applied to two dimensional infiltration problems. Water Resour. Res., 7(6): 1583-1596.

12. Guvanasen V and Volker R E 1983 Experimental investigations of unconfined aquifer pollution from recharge basins; Water Resour. Res. 19(3) 707717.

13. Oberhettinger and Baddi(1973) tables of Laplace transforms Springer- verlag, Newyork

14. Van Genuchten, M. Th. (1981). Analytical solutions for chemical transport with simultaneous adsorption, zero-order production and first order decay. Journal of Hydrology, 49(3-4), 213-233. 468

15. Zyvolski, G., Bruch, J.C., and James M. Sloss (1976): Solution of equation for two- dimensional infiltration problems. Soil Science 122 (2): 65-70.