Exact analytical solutions for two-dimensional advection–dispersion equation in cylindrical coordinates subject to third-type inlet boundary condition

Advances in Water Resources 34 (2011) 365–374

Contents lists available at ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

Exact analytical solutions for two-dimensional advection–dispersion equation in cylindrical coordinates subject to third-type inlet boundary condition Jui-Sheng Chen a,⇑, Yiu-Hsuan Liu a, Ching-Ping Liang b, Chen-Wuing Liu c, Chien-Wen Lin a a

Graduate Institute of Applied Geology, National Central University, Jhongli City, Taoyuan County 32001, Taiwan Department of Environmental Engineering and Science, Fooyin University, Ta-Liao Dist., Kaohsiung City 83102, Taiwan c Department of Bioenvironmental Systems Engineering, National Taiwan University, Taipei 10617, Taiwan b

a r t i c l e

i n f o

Article history: Received 26 August 2010 Received in revised form 17 December 2010 Accepted 18 December 2010 Available online 25 December 2010 Keywords: Advection–dispersion equation Analytical solution Finite Hankel transform Laplace transform Dispersion coefficient Third-type inlet boundary condition

a b s t r a c t Exact analytical solutions for two-dimensional advection–dispersion equation (ADE) in cylindrical coordinates subject to the third-type inlet boundary condition are presented in this study. The finite Hankel transform technique in combination with the Laplace transform method is adopted to solve the twodimensional ADE in cylindrical coordinates. Solutions are derived for both continuous input and instantaneous slug input. The developed analytical solutions are compared with the solutions for first-type inlet boundary condition to illustrate the influence of the inlet condition on the two-dimensional solute transport in a porous medium system with a radial geometry. Results show significant discrepancies between the breakthrough curves obtained from analytical solutions for the first-type and third-type inlet boundary conditions for large longitudinal dispersion coefficients. The developed solutions conserve the solute mass and are efficient tools for simultaneous determination of the longitudinal and transverse dispersion coefficients from a laboratory-scale radial column experiment or an in situ infiltration test with a tracer. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction The advection–dispersion equation (ADE) has been widely used to describe the movement of contaminants in the subsurface. Analytical solutions for ADE play important roles in giving initial or approximate estimates of contaminant concentration distributions in soil or aquifer systems, for determining transport parameters, and for validating more comprehensive numerical models. A number of analytical solutions for one-, two- and three-dimensional ADEs have been developed for predicting the transport of various contaminants in the subsurface. For example, van Genuchten and Alves [1] formulated several analytical solutions for the onedimensional ADE subject to various initial and boundary conditions. Batu [2,3] presented analytical solutions to two-dimensional ADE with various source boundary conditions. Leij et al. [4] and Park and Zhan [5] derived analytical solutions for three-dimensional ADE. Although a number of two- or three-dimensional analytical models are available at present, these models are mostly limited to the case of the ADE in Cartesian coordinates with steady uniform flow [5]. Analytical solutions for two-dimensional ADE in cylindrical coordinates are particularly useful for analyzing problems of the two-dimensional solute transport in a porous medium system with

⇑ Corresponding author. Tel.: +886 3 2807427; fax: +886 3 4263127. E-mail address: [email protected] (J.-S. Chen). 0309-1708/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2010.12.008

a radial geometry associated with a circular source. For example, they can be applied to interpret the field infiltration test with a tracer [6] or a two-dimensional laboratory column transport experiment [7,8] for simultaneous determining the longitudinal and transverse dispersion coefficients. Leij et al. [4] derived an analytical solution for two-dimensional ADE in cylindrical coordinates subject to the first- and third-type conditions using Laplace and Hankel transform technique by considering the medium is infinite in flow direction and in radial direction. Massabó et al. [9] has derived some analytical solutions to the two-dimensional ADE in cylindrical coordinates subject to the first-type condition using successive applications of the Laplace transform and the Fourier– Bessel series expansion by considering the medium is infinite in flow direction and finite in radial direction. Similar to Leij et al. [4], Zhang et al. [6] developed an analytical solution to the twodimensional ADE in cylindrical coordinates but considered a third-type instantaneous slug input. It should be noted that the analytical solutions derived by Massabó et al. [9] are based on the first-type inlet boundary condition which specifies a given solute concentration at the surface of the column. Previously work has shown that analytical solution for the first-type and third-type boundary conditions can be associated with flux- and volumeaverage concentration, respectively [10–13]. Analytical solution for the first-type inlet condition gives rise to physically improper mass conservation and significant errors in predicting the solute concentration distribution especially for a porous medium system with a large longitudinal dispersion coefficient if it is interpreted to

366

J.-S. Chen et al. / Advances in Water Resources 34 (2011) 365–374

represent the usual volume-average concentration. Selection of the appropriate inlet boundary condition has been the subject of much research such as solute transport either in a uniform flow or a radial flow field [14–20]. The purpose of this study is to develop exact analytical solutions for two-dimensional ADE in cylindrical coordinates subject to the third-type inlet boundary condition. The finite Hankel transform technique coupled with the Laplace transform method is applied to obtain the exact analytical solutions to the twodimensional ADE in cylindrical coordinates. The finite Hankel transform technique provides a systematic, efficient and straightforward approach for obtaining analytical solutions for both transient and steady flow transport problems with a radial geometry. The developed solutions are compared with the analytical solutions for the first-type inlet boundary condition to illustrate the impact of the inlet boundary condition. 2. Mathematical model This study considers the problem of two-dimensional solute transport in a porous media system with a radial geometry and a circular source as illustrated schematically in Fig. 1. The flow is steady and uniform along the downward column axis. The solute mass is injected over the inner zone of the upper surface. The injected solute migrates in the z direction by advection and longitudinal dispersion, whereas it spreads in the r direction by transverse dispersion. Solute transport in such a radial system can be described by ADE in cylindrical coordinates as

  @2C @C DT @ @C @C þ r ¼ DL 2  V @z @z @r @t r @r

VCðr; z ¼ 0; tÞ  DL

@Cðr; z ¼ 0; tÞ ¼ @z



VC 0

06r6q

0

q6r6R

ð3aÞ

and the instantaneous slug input as

@Cðr; z ¼ 0; tÞ ¼ VCðr; z ¼ 0; tÞ  DL @z

(

r dðtÞ /

06r6q

0

q6r6R

ð3bÞ

where C0 is the concentration of the applied solute solution (assumed to be constant) for continuous input; q is the radius of the inner inlet zone; r ¼ pMq2 is the solute mass distribution over the inner inlet zone for the instantaneous input; M is the injected solute mass. Another boundary condition for variable z required for obtaining a unique solution to Eq. (1) is imposed at infinity as follows:

Cðr; z ! 1; tÞ ¼ 0

ð4Þ

The boundary condition (4) indicates that the radial system is assumed to be of semi-infinite length. The boundary condition at r = R assumes impermeable condition as

@Cðz; r ¼ R; tÞ ¼0 @r

ð5Þ

ð1Þ 3. Derivation of the analytical solution

where C(z, r, t) denotes the solute concentration; t is the time; V stands for the averaged steady-state pore water velocity; and DL and DT represent the longitudinal and transverse dispersion coefficients, respectively. The system is free of solute at the initial time:

Cðz; r; t ¼ 0Þ ¼ 0 0 6 z 6 1; 0 6 r 6 R

where R is the radius of the radial system. The solute mass is applied over the inner zone of the upper surface. Based on mass conservation across the upper surface of the column system, the continuous input can be formulated as

ð2Þ

The aforementioned two-dimensional ADE in cylindrical coordinates and its associated initial and boundary conditions are analytically solved by successive implementation of the finite Hankel transform and the Laplace transform. First, utilizing the second kind of finite Hankel transform with Eq. (1) with respect to r results in

Fig. 1. Schematic representation of the two-dimensional advective–dispersive transport in a porous medium system with cylindrical geometry subject to steady and uniform flow field. The injected solute mass is applied over the inner zone of the upper surface.

367

J.-S. Chen et al. / Advances in Water Resources 34 (2011) 365–374

DL

@2CH @C H @C H V  DT k2n C H ¼ @z2 @z @t

ð6Þ

where kn is the finite Hankel transform parameter as determined by ðkn RÞ the transcendental equation dJ0 dr ¼ 0; J 0 ðÞ is the zero-order Bessel function of the first kind, CH(kn, z, t) is the second kind of finite Hankel transform for C(r, z, t) as defined by the following conjugate equations [21].

C H ðkn ; z; tÞ ¼ H½Cðr; z; tÞ ¼

Z

where a1 and a2 are two unknown coefficients which need to be determined. a1 and a2 in Eq. (15) can be straightforwardly determined with the aid of the boundary conditions found with Eqs. (13a), (13b) and (14). Introducing a1 and a2 into Eq. (15), the particular solution for continuous input can be expressed as follows:

C HL ðkn ; z; sÞ ¼

R

rCðr; z; tÞJ 0 ðkn rÞdr

ð7aÞ

0

Cðr; z; tÞ ¼

2 R2

1 2 X

C H ðkn ¼ 0; z; tÞ þ

R2

C H ðkn ; z; tÞ

n¼1

J 0 ðkn rÞ

ð7bÞ

jJ 0 ðkn RÞj2

  VC 0 Fðkn Þ Vz 1 pffiffiffiffiffiffi exp 2DL s DL  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V2 exp  pzffiffiffiffi 4D þ DT k2n þ s L DL qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  V2 pV ffiffiffiffi þ 4D þ DT k2n þ s L 2

DL

and that for instantaneous slug input as follows:

Accordingly, the initial and boundary conditions after the finite Hankel transform are

C H ðkn ; z; t ¼ 0Þ ¼ 0

ð8Þ

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V2  exp  pzffiffiffiffi 4D þ DT k2n þ s L DL r Fðkn Þ Vz qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C HL ðkn ; z; sÞ ¼ pffiffiffiffiffiffi exp 2DL / DL V2 pV ffiffiffiffi þ 4D þ DT k2n þ s L 

2

@C H ðkn ; z ¼ 0; tÞ ¼ VC 0 Fðkn Þ VC H ðkn ; z ¼ 0; tÞ  DL @z for continuous input

ð9aÞ

@C H ðkn ; z ¼ 0; tÞ r ¼ Fðkn ÞdðtÞ @z / for instantaneous slug input

ð9bÞ

C H ðkn ; z ! 1; tÞ ¼ 0 where Fðkn Þ ¼

ð10Þ

q2 2

q

J ðk kn 1 n

qÞ

if kn ¼ 0 ; and J1() is the first-order if kn –0

Applying the Laplace transform with respect to t, one has

 @ 2 C HL @C HL  V  DT k2n þ s C HL ¼ 0 @z2 @z

The solution for original domain C(r, z, t) is the consecutive inversion of the Laplace transform and the finite Hankel transform of CHL(kn, z, s). For convenience, the Laplace inversion is performed first. Eqs. (16a) and (16b) are expressed in different functional forms of the Laplace parameter s, the Laplace inversion is performed respectively for Eqs. (16a) and (16b). 3.1. Case I. Continuous input The Laplace inverse transform for Eq. (16a) for continuous input is obtained using thes-shift theorem and the following Laplace inverse formulas [1]

Bessel function of the first kind.

DL

ð11Þ

" 1

L

where s denotes the Laplace transform parameter and CHL(kn, z, s) represents the Laplace transform of CH(kn, z, t), which is defined by the following pair of equations

C HL ðkn ; z; sÞ ¼ L½C H ðkn ; z; tÞ ¼ C H ðkn ; z; tÞ ¼

1 2pi

Z

Z

1

C H ðkn ; z; tÞest dt

ð12aÞ

0 kþi1

C HL ðkn ; z; sÞest ds

ð12bÞ

where k is a line in the complex domain to the right of all poles. As a consequence, Eqs. (9a), (9b), and (10) after the Laplace transform become

pffiffi #   pffiffi expða sÞ 1 a 2 p ffiffi ¼  A t expðA t  A a Þerfc p ffiffi 2ðA þ BÞ 2 t ðs  A2 ÞðB þ sÞ   pffiffi 1 a 2  expðA t þ AaÞerfc pffiffi þ A t 2ðA  BÞ 2 t   pffiffi B a 2 p ffiffi þ B t þ 2 expðB t þ B a Þerfc 2 t A  B2 ð17Þ

"

ki1

@C HL ðkn ; z ¼ 0; sÞ VC 0 Fðkn Þ ¼ @z s for continuous input

L1

pffiffi # expða sÞ 1 expðA2 t pffiffi ¼ At þ 4A ðs  A2 ÞðA þ sÞ   pffiffi a 1  AaÞerfc pffiffi  A t  ð1 4A 2 t þ 2Aa þ 4A2 tÞ expðA2 t   pffiffi a þ AaÞerfc pffiffi þ A t 2 t

VC HL ðkn ; z ¼ 0; sÞ  DL

@C HL ðkn ; z ¼ 0; sÞ r VC HL ðkn ; z ¼ 0; sÞ  DL ¼ Fðkn Þ @z / for instantaneous slug input C HL ðkn ; z ! 1; sÞ ¼ 0

ð13aÞ

ð13bÞ ð14Þ

Eq. (11) is a homogeneous ordinary equation and has the following general solution:

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 3 20 V þ V 2 þ 4DL ðDT k2n þ sÞ Az5 þ a2 C HL ðkn ; z; sÞ ¼ a1 exp 4@ 2DL qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 3 20 V  V 2 þ 4DL ðDT k2n þ sÞ Az5  exp 4@ 2DL

DL

ð16bÞ

VC H ðkn ; z ¼ 0; tÞ  DL

(

ð16aÞ

Performing the s-shift theorem on Eq. (16a) and taking the value V2 a ¼ 4D þ DT k2n , one obtains L

C H ðkn ; z; tÞ ¼ L1 ½C HL ðkn ; z; sÞ " !# VC 0 Fðkn Þ Vz V2 2 ¼ pffiffiffiffiffiffi exp þ þ DT kn t L1 2DL 4DL DL   8 9 pffiffi > > zffiffiffiffi p > > exp  s < = DL    h  2 i p ffiffi > > > V pV ffiffiffiffi þ s > : s  4D ; þ DT k2n L

ð15Þ

ð18Þ

2

ð19Þ

DL

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V2 Taking Eq. (17) on Eq. (19) and letting A ¼ 4D þ DT k2n ; B ¼ pV ffiffiffiffi L 2 DL z and a ¼ pffiffiffiffi, we obtain the following inversion DL

368

J.-S. Chen et al. / Advances in Water Resources 34 (2011) 365–374

  

   V V U z  Ut z erfc pffiffiffiffiffiffiffi C H ðkn ; z; tÞ ¼ C 0 Fðkn Þ exp V þU 2DL 2 DL t  

   V V þU z þ Ut z erfc pffiffiffiffiffiffiffi exp þ V U 2DL 2 DL t  )   V2 Vz z þ Vt 2 p ffiffiffiffiffiffiffi þ exp  D k t erfc T n DL 2 DL t 2DL DT k2n

Table 1 Descriptive simulation conditions and transport parameters. Parameter

ð20Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2ffi where U ¼ 1 þ 4DLVD2T kn .

It should be noted that Eq. (20) is only valid for kn – 0. For kn = 0, Eq. (16a) can be simplified into the following form

C H ðkn ; z; tÞ ¼ L1 ½C HL ðkn ¼ 0; z; sÞ " # VC 0 Fðkn ¼ 0Þ Vz V 2t pffiffiffiffiffiffi exp þ t L1 ¼ 2DL 4DL DL   9 8 pffiffi > > zffiffiffiffi p > > exp  s < = DL    h i p ffiffi > > 2 > V pV ffiffiffiffi þ s > : s  4D ; 2

L

V Taking Eq. (18) on Eq. (21) and letting A ¼ 2D and a ¼ pzffiffiffiffi, we obL DL tain the following inversion:

C H ðkn ¼ 0; z; tÞ

ð22Þ

The original domain solution can be straightforwardly obtained using the finite Hankel inversion of Eq. (7b) and is expressed as

8 sffiffiffiffiffiffiffiffiffi " #   C 0 q2 <1 z  Vt V 2t ðz  VtÞ2 erfc pffiffiffiffiffiffiffi þ Cðr; z; tÞ ¼ 2 exp  pDL 4DL t 2 DL t R :2 !  )   2 1 Vz V t Vz z þ Vt exp  erfc pffiffiffiffiffiffiffi 1þ þ 2 DL DL DL 2 DL t   

   1 2C 0 q X 1 V V U z  Ut z erfc pffiffiffiffiffiffiffi exp þ 2 2DL 2 DL t R n¼1 kn V þ U  

   V V þU z þ Ut z erfc pffiffiffiffiffiffiffi exp þ V U 2DL 2 DL t  )   V2 Vz z þ Vt 2 ffiffiffiffiffiffiffi p þ exp  DT kn t erfc DL 2 DL t 2DL DT k2n 

1 0.4 0.1 0.2 0.1

0.01

0.01

0.001

1 1

Table 2 Solutions for continuous input at z = 0.8 m and r = 0 m (N = number of terms summed). DL = 1 m2/h and DT = 0.1 m2/h t

N=2

N=3

N=4

N=5

N=6

0.08 0.88 1.68 2.48 3.28 4.08

3.0720E03 4.1080E02 5.4486E02 6.1876E02 6.6378E02 6.9284E02

3.1678E03 4.1190E02 5.4597E02 6.1987E02 6.6489E02 6.9395E02

3.1675E03 4.1190E02 5.4597E02 6.1986E02 6.6489E02 6.9395E02

3.1674E03 4.1190E02 5.4596E02 6.1986E02 6.6489E02 6.9395E02

3.1674E03 4.1190E02 5.4596E02 6.1986E02 6.6489E02 6.9395E02

DL = 0.1 m2/h and DT = 0.01 m2/h t

N=4

N=8

N = 12

N = 16

N = 20

0.08 0.88 1.68 2.48 3.28 4.08

2.0974E09 2.0182E01 2.7530E01 2.8157E01 2.8210E01 2.8215E01

1.9178E09 2.0177E01 2.7525E01 2.8151E01 2.8205E01 2.8210E01

1.9239E09 2.0177E01 2.7525E01 2.8151E01 2.8205E01 2.8210E01

1.9238E09 2.0177E01 2.7525E01 2.8151E01 2.8205E01 2.8210E01

1.9238E09 2.0177E01 2.7525E01 2.8151E01 2.8205E01 2.8210E01

DL = 0.01 m2/h and DT = 0.001 m2/h t

N = 20

N = 30

N = 40

N = 50

N = 60

0.08 0.88

3.2763E73 7.0107E01

3.2102E73 7.0107E01

3.2075E73 7.0107E01

3.2076E73 7.0107E01

3.2076E73 7.0107E01

Table 3 Solutions for instantaneous slug input at z = 0.8 m and r = 0 m (N = number of terms summed). DL = 1 m2/h and DT = 0.1 m2/h

J 1 ðkn qÞJ 0 ðkn rÞ

ð23Þ

jJ0 ðkn RÞj2

The steady state solution can be obtained by invoking the Tauberian theorem [21, p. 187, Eq. (3-13-3)] on Eq. (16a) and using the finite Hankel inverse transform as

Cðr; zÞ ¼

Average velocity V (m/h) Radius of cylindrical system R (m) Radius of inner injected zone q (m) Porosity / (m) Longitudinal dispersion coefficient DL (m2/h) 1 Transverse dispersion coefficient 2 DT (m /h) 0.1 Injected mass for instantaneous slug input M (kg) Injected concentration for continuous input C0 (kg/m3)

ð21Þ

DL

8 sffiffiffiffiffiffiffiffiffi " #   <1 z  Vt V 2t ðz  VtÞ2 exp  erfc pffiffiffiffiffiffiffi þ ¼ C 0 Fðkn Þ :2 pDL 4DL t 2 DL t !    ) 1 Vz V 2 t Vz z þ Vt  exp erfc pffiffiffiffiffiffiffi 1þ þ 2 DL DL DL 2 DL t

Values

1 2 X VC 0 Fðkn Þ pffiffiffiffiffiffi 2 DL R n¼0

t

N=2

N=3

N=4

N=5

N=6

0.08 0.16 0.24 0.32 0.40

8.2585E02 9.1084E02 6.8502E02 5.2715E02 4.3249E02

9.6844E02 9.3239E02 6.8697E02 5.2730E02 4.3250E02

9.7632E02 9.3246E02 6.8697E02 5.2730E02 4.3250E02

9.7640E02 9.3246E02 6.8697E02 5.2730E02 4.3250E02

9.7640E02 9.3246E02 6.8697E02 5.2730E02 4.3250E02

DL = 0.1 m2/h and DT = 0.01 m2/h

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V2   exp  pzffiffiffiffi 4D þ DT k2n L DL Vz J 0 ðkn rÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  exp 2 2 2DL V2 jJ pV ffiffiffiffi þ 4D 0 ðkn RÞj þ D k T n L 2

t

ð24Þ

DL

3.2. Case II. Instantaneous slug input The Laplace inverse transform for instantaneous slug input in Eq. (16b) can be achieved using the s-shift theorem and the following Laplace inverse formula [1]

N = 12

N = 16

N = 20

0.08 5.4828E07 5.0487E07 0.16 4.9904E03 4.9121E03 0.24 7.0205E02 6.9950E02 0.32 2.0569E01 2.0550E01 0.40 3.2936E01 3.2927E01 DL = 0.01 m2/h and DT = 0.001 m2/h

N=4

N=8

5.0612E07 4.9122E03 6.9950E02 2.0550E01 3.2927E01

5.0611E07 4.9122E03 6.9950E02 2.0550E01 3.2927E01

5.0611E07 4.9122E03 6.9950E02 2.0550E01 3.2927E01

t

N = 20

N = 30

N = 40

N = 50

N = 60

0.08 0.16 0.24 0.32

8.1652E70 1.8902E27 5.7684E14 1.0867E07

8.0021E70 1.8864E27 5.7674E14 1.0867E07

7.9957E70 1.8864E27 5.7674E14 1.0867E07

7.9960E70 1.8864E27 5.7674E14 1.0867E07

7.9960E70 1.8864E27 5.7674E14 1.0867E07

369

J.-S. Chen et al. / Advances in Water Resources 34 (2011) 365–374

pffiffi 

  expða sÞ 1 a2 pffiffi ¼ pffiffiffiffiffiffi exp  L1  A expðA2 t 4t Aþ s pt   pffiffi a þ AaÞerfc pffiffi þ A t 2 t

C H ðkn ; z; tÞ ¼ ð25Þ

Imposing s-shift theorem on Eq. (16b) and taking the value V2 a ¼ 4D þ DT k2n , we get L

L1 ½C HL ðkn ; z; sÞ ¼

"

rFðkn Þ

!#

2

Vz V pffiffiffiffiffiffi exp þ þ DT k2n t L1 2DL 4DL / DL   2 pffiffi 3 exp pzffiffiffiffi s 6 7 DL 7 6 4 pV ffiffiffiffi þ pffiffis 5 2

0

0.4

1.2

1.6

ð27Þ

( " # 1 X 1 ðVt  zÞ2 exp   DT k2n t k 4DL t n¼0 n    ) V Vz z þ Vt 2 exp  DT kn t erfc pffiffiffiffiffiffiffiffiffiffi  2DL DL 4DL t ! q2 2qJ1 ðkn qÞJ0 ðkn rÞ þ  R2 R2 jJ 0 ðkn RÞj2

r

ð26Þ

2

0

0.12

a

1 ðVt  zÞ2 Fðkn Þ pffiffiffiffiffiffiffiffiffiffiffi exp   DT k2n t / 4DL t DL pt  )   V Vz z þ Vt 2 exp  DT kn t erfc pffiffiffiffiffiffiffiffiffiffi  2DL DL 4DL t

r

#

Cðr; z; tÞ ¼ pffiffiffiffiffiffiffiffiffiffiffi / DL pt

DL

0.8

"

The solution in original domain can be obtained using Eq. (7b) and is expressed as

Making use of Eq. (25) on Eq. (26) and letting A ¼ pV ffiffiffiffi and B ¼ pzffiffiffiffi 2 DL DL gives the following inversion:

0.12

(

0.4

0.8

1.2

2 0.5

0.5

b

analytical solution numerical solution

analytical solution numerical solution

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0.08

C/CR

C/CR

0.08

1.6

ð28Þ

0.04

0.04

0

0 0

0.4

0.8

1.2

1.6

0

2

0 0

0.4

0.8

1.2

t [hr]

1.6

2

t [hr] 0

C/CR

4

0.4

0.8

1.2

1.6

2 4

c

analyitcal solution numerical solution

3

3

2

2

1

1

0

0 0

0.4

0.8

1.2

1.6

2

t [hr]

Fig. 2. The comparison of the breakthrough curves obtained from the derived analytical solution and those from the corresponding numerical solution at z = 0.8 m and r = 0 m for instantaneous slug input under various dispersion coefficient values: (a) DL = 1 m2/h and DT = 0.1 m2/h; (b) DL = 0.1 m2/h and DT = 0.01 m2/h; (c) DL = 0.01 m2/h and DT = 0.001 m2/h. The numerical solution is obtained using Laplace transform method [24].

370

J.-S. Chen et al. / Advances in Water Resources 34 (2011) 365–374

4. Results and discussion 4.1. Verification of the developed solution The analytical solutions for continuous input and instantaneous input are both in the form of the sum of the infinite series expansion and can be straightforwardly evaluated. It is interesting to look at how many terms are required to numerically compute the accurate solution. The parameter values for the numerical results for Eqs. (24) and (28) are summarized in Table 1. Tables 2 and 3 illustrate the convergence of the solution for continuous input and instantaneous slug input computed at z = 0.8 m and r = 0 m under different longitudinal dispersion coefficients. The values for parameter N in Tables 2 and 3 are defined as the number of terms summed in the truncated series expansion. Clearly, the term required to obtain an accurate numerical

0

2

4

6

0.05

8

0

2

4

6

0.05

a

current study Leij et al. (1991)

result is strongly dependent on the longitudinal dispersion coefficient. The analytical solution requires N = 6, 20 and 60 terms to achieve convergence to 4 decimal places for dimensionless Peclet  number Pe ¼ DVzL equal to 0.8, 8, and 80. For the case of large Peclet number, the required terms summed in the truncated series expansion is large. Yeh and Yang [22] suggested using Shanks [23] method to accelerate the convergence of the infinite series solution. After determining value of the parameter N, the accuracy and correctness of the solution are examined by comparing with the corresponding numerical solution. The numerical solution is generated using the Laplace transform finite difference (LTFD) technique proposed by Moridis and Reddel [24]. The LTFD technique has several advantages over the conventional timemarching finite difference method. The parameter values are the same as those listed in Table 1. The shapes of the breakthrough curves are strongly influenced by the Peclet number. In

b

current study Leij (1991) 0.04

0.03

0.03

0.02

0.02

0.3

0.3

0.2

0.2

0.1

0.1

C/C0

C/C0

0.04

8

0.01

0.01

0

0 0

2

4 t [hr]

6

0

0

8

2

0 0

4

2

6

6

8

8

c

current study Leij (1991)

1

C/C0

1

4 t [hr]

0.5

0.5

0

0 0

2

4 t [hr]

6

8

Fig. 3. The comparison of the breakthrough curves obtained from the derived analytical solution and those from Leij et al. [4] solution at z = 0.8 m and r = 0 m for continuous input under various dispersion coefficient values: (a) DL = 1 m2/h and DT = 0.1 m2/h; (b) DL = 0.1 m2/h and DT = 0.01 m2/h; (c) DL = 0.01 m2/h and DT = 0.001 m2/h.

371

J.-S. Chen et al. / Advances in Water Resources 34 (2011) 365–374

many instances, an arithmetic operation cannot be performed in the programming for the analytical solution when the Peclet number is large (advection-dominated). Such conditions generally involve breakthrough curves with a steep front. We therefore perform a series of comprehensive comparison tests for a wide range of Peclet numbers. Fig. 2 depicts the breakthrough curves observed at r = 0 m and z = 0.8 m obtained from the developed analytical solution and the numerical solution for an instantaneous slug input under different Peclet number. The developed analytical solutions agree well with the numerical solutions for all Peclet numbers. Additionally, comparing the developed solution with the Leij et al. [4] solution provides another way to verify the accuracy and correctness of the developed solution. Because the Leij et al. [4] solution considered a medium with an infinite radial distance, a large outer radius (R = 20 m) for

0 0.15

1

2

3

4

our solution is used for the comparison. The numerical result of Leij et al. [4] solution is produced using 3DADE program embedded in the windows-based STANMOD code [25]. Fig. 3 depicts the breakthrough curves observed at r = 0 m and z = 0.8 m obtained from the developed analytical solution and the Leij et al. [4] solution for a continuous input under different Peclet number. As expected, two solutions are nearly identical for various dispersion coefficients. 4.2. Effect of the inlet boundary condition The analytical solution is compared to the solution derived for the first-type inlet boundary condition to illustrate the effect of the inlet boundary condition on the two-dimensional solute transport in a radial porous medium system.

5

0

0.15

a

0.4

2

3

4

5 0.4

b

first-type BC third-type BC

first-type BC third-type BC

0.3

0.3

0.2

0.2

0.1

0.1

0.1

C/C0

C/C0

0.1

1

0.05

0.05

0

0

0 0

1

2

3

4

0 0

5

1

2

3

4

5

t [hr]

t [hr]

0

1

2

3

c

4

5

first-type BC third-type BC 1

C/C0

1

0.5

0.5

0

0 0

1

2

3

4

5

t [hr]

Fig. 4. The comparison of the breakthrough curves obtained from the analytical solution for the first-type and third-type boundary conditions at z = 0.8 m and r = 0 m for continuous input under various dispersion coefficient values: (a) DL = 1 m2/h and DT = 0.1 m2/h; (b) DL = 0.1 m2/h and DT = 0.01 m2/h; (c) DL = 0.01 m2/h and DT = 0.001 m2/h.

372

J.-S. Chen et al. / Advances in Water Resources 34 (2011) 365–374

   1 X Vz r J 0 Z k;1 Cðr; z; tÞ ¼ C 0 Ak exp 2D R 8 k¼0 2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi3  2 2 < V L DT 2 5 z  exp 4z þ Z k;1  erfc pffiffiffiffiffiffiffi 2 2 : 2 DL t 4D DL R 2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi3 s ! Z 2k;1 V 2t V 2 L2 DT 2 5 þ þ DT 2 t þ exp 4z þ Z k;1 erfc 4DL R 4D2 DL R2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi19 0 Z 2k;1 = z V 2t  @ pffiffiffiffiffiffiffi  þ DT 2 t A ð31Þ ; 4DL 2 DL t R

For the inlet boundary, the first-type inlet boundary condition is formulated as [8]

 Cðr; z ¼ 0; tÞ ¼

C0

06r6q

0

q6r6R

ð29aÞ

and for instantaneous slug input as

Cðr; z ¼ 0; tÞ ¼

(r dðtÞ 0 6 r 6 q /V

ð29bÞ

q6r6R

0

Massabó et al. [9] derived an analytical solution for the first-type boundary condition subject to continuous input;

0 0.6

0.4

0.8

1.2

1.6

where Zk,1 is the kth root of the Bessel function of first order; the Ak are coefficients defined as follows [9]:

2

0 0.6

a

0.6

0.2

0.2

1.2

1.6

2 0.6

first-type BC third-type BC

0.4

0.4

0.2

0.2

C/CR

C/CR

0.4

0.8

b

first-type BC third-type BC

0.4

0.4

0

0 0

0.4

0.8

1.2

1.6

0

2

0 0

0.4

0.8

t [hr]

1.6

2

t [hr]

0 4

C/CR

1.2

0.4

0.8

1.2

1.6

2 4

c

first-type BC third-type BC

3

3

2

2

1

1

0

0 0

0.4

0.8

1.2

1.6

2

t [hr]

Fig. 5. The comparison of the breakthrough curves obtained from the analytical solution for the first-type and third-type boundary conditions at z = 0.8 m and r = 0 m for instantaneous slug input under various dispersion coefficient values: (a) DL = 1 m2/h and DT = 0.1 m2/h; (b) DL = 0.1 m2/h and DT = 0.01 m2/h; (c) DL = 0.01 m2/h and DT = 0.001 m2/h.

J.-S. Chen et al. / Advances in Water Resources 34 (2011) 365–374

Ak ¼

r q 2 / R

Ak ¼ 2

rq

for k ¼ 0   J 1 qR Z k;1

/ R Z k;1 ½J 0 ðZ k;1 Þ2

ð32Þ for k – 0

ð33Þ

Massabò et al. [9] derived an exact analytical solution for instantaneous slug input subjected to first-type condition by application of the Bessel series expansion technique in combination with the Fourier transform method. In dealing with the instantaneous slug input, Massabò et al. [9] specified it as a Dirac delta function and considered the transport domain to be infinite (1 6 z 6 +1). However, it has been shown by Coronado et al. [26], that this condition gives rise to physically improper backward dispersion (i.e. the solute mass can move upwardly outside of the column system, 1 6 z 6 0). This is obviously contradictory to the actual transport process in the soil column. The actual transport process is that the tracer mass is injected at the upper surface of column (z = 0) and migrates downwardly in the flow direction (0 6 z 6 1). By using boundary conditions (29b) and following the mathematical procedures of the Massabò et al. [9] analytical solution for continuous input we can obtain the analytical solution for instantaneous input within a semi-infinite domain without any difficulty. The newly obtained solution is different from Massabò et al. [8] solution by only a factor of Vtz and expressed as follows:

Cðr; z; t ¼ 0Þ

" # þ1 X z Ak ðz  VtÞ2 pffiffiffiffiffiffiffiffiffiffiffiffiffi exp  ¼ Vt 4ptDL 4ptDL k¼0 " # 2  Z k;1 r  exp DT 2 t J0 Z k;1 R R

ð34Þ

Figs. 4 and 5 show the comparisons of the breakthrough curves at z = 0.8 m and r = 0 m for continuous input and instantaneous slug input, respectively, for the first- and third-type boundary conditions for different Peclet number. The used parameters are the same as Table 1. These graphs clearly show that there is a pronounced discrepancy in the predicted breakthrough curves for the solution obtained with the first-type boundary condition and the solution for the third- type boundary condition for a small Peclet number. The solution for the first-type boundary conditions gives a concentration that is more than 1.6 times higher for continuous input and 5.0 times higher concentration for instantaneous slug input than the solution for the third-type boundary condition. This will lead to a significant error in determining the longitudinal and transverse dispersion coefficients if the solutions for the first-type inlet boundary conditions are used to interpret the results of the radial column experiment. The discrepancies between solutions for first- and third-type condition decrease with increasing Peclet number. For Peclet number equal to 80, the inlet boundary condition has negligible effects on the BTCs. Accordingly, the analytical solution derived for first-type inlet boundary condition should be cautiously applied to interpret the resident concentration when a system having a large dispersion coefficient. 5. Conclusions Previous study has demonstrated that the third-type inlet boundary solution correctly conserves solute mass and accurately predicts solute movement in the problem of solute transport in a porous medium system. Analytical solutions for two-dimensional solute transport in a porous medium system with a radial geometry subject to the third-type boundary condition are presented in this study. The two-dimensional advection–dispersion equation in cylindrical coordinates is solved using a finite Hankel transform coupled with the Laplace transform. The developed analytical solu-

373

tions are compared with the analytical solutions for first-type boundary condition to clarify how the inlet boundary conditions influence the two-dimensional transport in a porous medium system with a radial geometry. Results show that the first-type boundary condition solution departs from the third-type boundary condition solution especially for a system with a small Peclet number. The analytical solutions derived for the third-type boundary condition conserve the mass and will thus be particularly useful for analyzing problems such as the two-dimensional transport in a radial soil column experiment or an infiltration test with a tracer. Acknowledgement The authors thank the National Science Council of the Republic of China for financially supporting this work under Contract Nos. NSC 97-2116-M-008-011 and NSC 98-2628-M-008-010. References [1] van Genuchten MTh, Alves WJ. Analytical solutions of the one-dimensional convective–dispersive solute transport equation. US Department of Agriculture, Technical Bulletin No. 1661; 1982. 151 pp. [2] Batu V. A generalized two-dimensional analytical solution for hydrodynamic dispersion in bounded media with the first-type boundary condition at the source. Water Resour Res 1989;25(6):1125–32. [3] Batu V. A generalized two-dimensional analytical solute transport model in bounded media for flux-type finite multiple sources. Water Resour Res 1993;29(8):2881–92. [4] Leij FJ, Skaggs TH, van Genuchten MTh. Analytical solution for solute transport in three-dimensional semi-infinite porous media. Water Resour Res 1991;27(10):2719–33. [5] Park E, Zhan H. Analytical solutions of contaminant transport from finite one-, two, three-dimensional sources in a finite-thickness aquifer. J Contam Hydrol 2001;53(1–2):41–61. [6] Zhang X, Qi X, Zhou X, Pang H. An in situ method to measure the longitudinal and transverse dispersion coefficients of solute transport in soil. J Hydrol 2006;328(3–4):614–9. [7] Massabò M, Catania F, Paladino O. A new method for laboratory estimation of the transverse dispersion coefficient. Ground Water 2007;45(3):339–47. [8] Frippiat CC, Perez PC, Holeyman AE. Estimation of laboratory-scale dispersivities using an annulus-and-core device. J Hydrol 2008;362(1– 2):57–68. [9] Massabò M, Cianci R, Paladino O. Some analytical solutions for twodimensional convection–dispersion equation in cylindrical geometry. Environ Modell Softw 2006;21(5):681–8. [10] van Genuchten MTh, Parker JC. Boundary conditions for displacement experiments through short laboratory soil columns. Soil Sci Soc Am J 1984;48(4):703–8. [11] Parker JC, van Genuchten MTh. Flux-averaged and volume-averaged concentrations in continuum approaches to solute transport. Water Resour Res 1984;20(7):866–72. [12] Kreft A, Zuber A. Comment on ‘‘Flux averaged and volume averaged concentrations in continuum approaches to solute transport’’. Water Resour Res 1986;22:1157–8. [13] Parlange JY, Barry DA, Starr JL. Comments on ‘‘Boundary conditions for displacement experiments through short laboratory soil columns’’. Soil Sci Soc Am J 1985;49(5):1325. [14] Barry DA, Sposito G. Application of the convection–dispersion model to solute transport in finite soil columns. Soil Sci Soc Am J 1988;52(1):3–9. [15] Parlange JY, Starr JL, van Genuchten MTh, Barry DA, Parker JC. Exit condition for miscible displacement experiments in finite columns. Soil Sci 1992;153(3):165–71. [16] Chen JS, Ni CF, Liang CP, Chiang CC. Analytical power series solution for contaminant transport with hyperbolic asymptotic distance-dependent dispersivity. J Hydrol 2008;362(1–2):142–9. [17] Chen JS, Ni CF, Liang CP. Analytical power series solutions to the twodimensional advection–dispersion equation with distance-dependent dispersivities. Hydrol Process 2008;22(24):4670–8. [18] Chen CS. Analytical solutions for radial dispersion with Cauchy boundary at injection well. Water Resour Res 1987;23(7):1217–24. [19] Chen JS, Chen CS, Chen CY. Analysis of solute transport in a divergent flow tracer test with scale-dependent dispersion. Hydrol Process 2007;21(18):2526–36. [20] Chen JS, Liu CW, Liao CM. Two-dimensional Laplace transformed power series solution for solute transport in a radially convergent flow field. Adv Water Resour 2003;26(10):1113–24. [21] Sneddon IH. The use of integral transforms. New York: McGraw-Hill; 1972. [22] Yeh HD, Yang SY. A new method for laboratory estimation of the transverse dispersion coefficient-discussion. Ground Water 2010;48(1):16–7.

374

J.-S. Chen et al. / Advances in Water Resources 34 (2011) 365–374

[23] Shanks D. Non-linear transformations of divergent and slowly convergent sequences. J Math Phys 1955;34:1–42. [24] Moridis GJ, Reddell DL. The Laplace transform finite difference method for simulation of flow through porous media. Water Resour Res 1991;27(8):1873–84. [25] Simunek J, van Genuchten MTh, Sejna M, Toride N, Leij FJ. The STANMOD computer software for evaluating solute transport in porous media using

analytical solutions of convection–dispersion equation. Riverside, CA: US Salinity Laboratory, Agriculture Research Service, US Department of Agriculture; 1999. [26] Coronado M, Ramírez J, Samaniego F. New considerations on analytical solutions employed in tracer flow modeling. Transport Porous Med 2004;54(2):221–37.