Two-dimensional power series solution for non-axisymmetrical transport in a radially convergent tracer test with scale-dependent dispersion

Advances in Water Resources 30 (2007) 430–438 www.elsevier.com/locate/advwatres

Two-dimensional power series solution for non-axisymmetrical transport in a radially convergent tracer test with scale-dependent dispersion Jui-Sheng Chen

*

Institute of Applied Geology, National Central University, No. 300, Jungda Rd., Jhongli, Taoyuan 320, Taiwan, ROC Received 31 October 2005; received in revised form 8 May 2006; accepted 11 May 2006 Available online 7 July 2006

Abstract It has been known for many years that dispersivity increases with solute travel distance in a subsurface environment. The increase of dispersivity with solute travel distance results from the significant variation of hydraulic properties of heterogeneous media and was identified in the literature as scale-dependent dispersion. This study presents an analytical solution for describing two-dimensional non-axisymmetrical solute transport in a radially convergent flow tracer test with scale-dependent dispersion. The power series technique coupling with the Laplace and finite Fourier cosine transform has been applied to yield the analytical solution to the two-dimensional, scale-dependent advection–dispersion equation in cylindrical coordinates with variable-dependent coefficients. Comparison between the breakthrough curves of the power series solution and the numerical solutions shows excellent agreement at different observation points and for various ranges of scale-related transport parameters of interest. The developed power series solution facilitates fast prediction of the breakthrough curves at any observation point. Ó 2006 Elsevier Ltd. All rights reserved. Keywords: Scale-dependent dispersion; Convergent flow field; Power series solution

1. Introduction To describe solute transport in a subsurface porous medium, the advection–dispersion equation (ADE) is commonly used to describe the physical process governing advective, molecular diffusive and hydrodynamic dispersive transport. In the theory of dispersive transport, the dispersivity is a controlling parameter for the dilution of concentrations in a solute plume that is displaced by groundwater flow in an aquifer. Research works in groundwater contamination have shown that the solutions of the ADE give a satisfactory agreement with experimental data only if the value of dispersivity needed to fit the data to an ADE are increased with solute displacement distance [14]. The increase of dispersivity with solute travel distance was iden*

Tel.: +886 3 2807427; fax: +886 3 4263127. E-mail address: [email protected]

0309-1708/$ - see front matter Ó 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2006.05.004

tified in the literature as scale-dependent dispersion. Notably, a few researchers have presented analytical solutions to ADE in Cartesian coordinates to describe scale-dependent dispersive solute transport in a unidirectional flow field with distance-dependent dispersivity [8,13,18,23,24]. Su et al. [22] presented analytical solutions to advection–dispersion equation with a spatially and temporally varying dispersion coefficient for modelling solute transport in steady, saturated subsurface flow through heterogeneous porous media. Few researchers have developed solutions to the two-dimensional advection–dispersion equation with scale-dependent dispersion in a uniform flow field. Hunt [9] developed one-, two- and three-dimensional analytical solutions for a scale-dependent dispersion equation. Hunt [10] used a singular perturbation method to obtain a twodimensional approximate analytical solution for scaledependent dispersive transport in a uniform flow field. Alternatively, Aral and Liao [1] developed analytical

J.-S. Chen / Advances in Water Resources 30 (2007) 430–438

solutions for two-dimensional advection–dispersion equation with a time-dependent dispersion coefficient. The ADE in cylindrical coordinates refers to the problem of analyzing the dispersive transport of a contaminant in the radial flow field generated by the extraction or injection well that fully penetrates a homogeneous confined aquifer. Current contaminant transport models in a radial flow field are based on ADE with constant dispersivity [2– 4,15]. Relatively few studies have developed solutions to ADE in cylindrical coordinates to investigate how scaledependent dispersion affects solute transport in a radial flow field. Kocabas and Islam [11] developed analytical solutions to the radial advection–dispersion equation, assuming spatial velocity and scale dependence of dispersivity for a radially divergent flow field. Chen et al. [6] presented a semi-analytical solution for describing the transport of dissolved substances in a radially convergent flow field with linearly scale-dependent dispersivity. Their solution was applied to a set of reported field data to examine its applicability. The application results reveal that the solute transport process at the test site obeys the linearly scale-dependent dispersion assumption. To my knowledge, the analytical solution to the twodimensional, scale dependent advection–dispersion equation in cylindrical coordinates is not currently available. The lack of progress in this area is attributed to difficulty in deriving the analytical solution to ADE with variable coefficients that results from the dependency of both the longitudinal and transverse dispersivities on solute travel distance. Chen et al. [5] presented a two-dimensional mathematical model in cylindrical coordinates that was solved by the Laplace transform finite difference (LTFD) method to illustrate the salient feature of scale-dependent dispersion for non-axisymmetrical transport in a radially convergent flow field. In addition, a curve-fitting method involving a theoretical breakthrough curves at the extraction well and one observation well was proposed to simultaneously evaluate the dispersivity/distance ratios in longitudinal and transverse directions, respectively. However, the numerical solution has the propensity to be computationally cumbersome and expensive and, in absence of interpolation, can be used to yield breakthrough curves only at discrete nodal points. Conversely, analytical solutions can inherently provide a temporally and spatially continuous concentration distribution and are extremely useful in generating theoretical breakthrough curves at any observation point with small computational efforts. Analytical solutions also make them valuable as screening tools or benchmarks for the more complex numerical codes. In this paper the problem is approached analytically to yield a temporally and spatially continuous solution. A power series technique coupling with Laplace and finite Fourier cosine transform was applied to obtain the analytical solution to the two-dimensional, scale dependent advection–dispersion equation in cylindrical coordinates. The developed analytical solution will be compared to the numerical solution to examine its accuracy and robustness.

431

2. Problem formulation An ideal convergent tracer test system considered herein is depicted in Fig. 1. An extraction well with a radius of rW located along vertical axis at (r = 0, h = 0), screened entirely over the statistically homogeneous and isotropic aquifer of constant thickness, withdraws water from the aquifer at a constant rate of Q. The initial transient phase is neglected and the steady and axisymmetrical flow condition is established prior to start of tracer injection. The extraction-induced mean pore water velocity (V) is given by V ¼

A r

ð1Þ

Q where A ¼ 2pb/ , and b and / represent the aquifer thickness and effective porosity. When a field test is initiated, once water levels are stabilized, a tracer is introduced into the injection borehole with a radius of rI located at (r = rL, h = p), from which it flows out of the injection borehole. Unlike a diverging flow tracer test, the plume is non-axisymmetrical during tracer transport into extraction well in a converging flow tracer test (Fig. 1). The two-dimensional advection–dispersion equation in cylindrical coordinates for non-axisymmetrical solute transport in a tracer test system as shown in Fig. 1 can be described as     oC 1 o oC oC 1 o oC R ¼ rDL þ DT V ð2Þ ot r or or or r2 oh oh

where C(r, h, t) is the concentration of a tracer, r and h are the cylindrical coordinates, t is time, DL and DT denote hydrodynamic dispersion coefficients in the longitudinal

Fig. 1. Schematic diagram of radically convergent tracer test: (A) side view, (B) plane view.

432

J.-S. Chen / Advances in Water Resources 30 (2007) 430–438

and transverse directions, respectively, and R is retardation factor. Hydrodynamic dispersion includes mechanical dispersion and molecular diffusion. Molecular diffusion is generally much smaller than mechanical dispersion and can be neglected in a subsurface groundwater system. The longitudinal and transverse dispersion coefficients, DL and DT, are generally considered to be a linear function of the pore-water velocity and can be expressed in the following form: DL ¼ aL jV j DT ¼ aT jV j

ð3aÞ ð3bÞ

where aL and aT are longitudinal and transverse dispersivities, respectively. Field investigations have indicated that longitudinal dispersivity in fact is not constant but generally increases with the travel distance from the contaminant source. Although a detailed dispersivity–distance relationship is not available for each site, it is reasonable to postulate the longitudinal dispersivity increase linearly with solute distance from the injection well located at r = rL. The growth of the longitudinal dispersivity can be adequately described using an equation of the form (see [20]) aL ðrÞ ¼ eL ðrL  rÞ

ð4aÞ

eL is longitudinal dispersivity/distance ratio (i.e. the linearly scale-proportional factor), which characterizes the linearly scale-dependent dispersion process. The relationship between the longitudinal and transverse scale-dependent dispersivities may be assumed to be constant as that in a constant dispersivity model. Several researchers including Hunt [9,10] and Aral and Liao [1] extended this constant relationship of dimensionless dispersivity ratio to describe the transverse scale-dependent dispersivity in the development of their analytical solution. Following Hunt [9,10] the transverse scaledependent dispersivity, aT(r), can be expressed in the following form: aT ðrÞ ¼ X s aL ðrÞ ¼ X s eL ðrL  rÞ ¼ eT ðrL  rÞ

ð4bÞ

where X s ¼ aaTL is the dimensionless dispersivity ratio, eT is transverse dispersivity/distance ratio. Substituting Eqs. (1), (3a), (3b), (4a) and (4b) into Eq. (2) yields   oC 1 o A oC A oC ¼ reL ðrL  rÞ R þ ot r or r or r or   1 o A oC eT ðrL  rÞ þ 2 ð5Þ r oh r oh Rewriting Eq. (5) in more clearly visible form as R

oC eL A o2 C Að1  eL Þ oC eT AðrL  rÞ o2 C ¼ ðrL  rÞ 2 þ þ ot r r or r3 or oh2 ð6Þ

Fig. 2. Plane view of the injection well boundary condition (after Zlotnik and Logan [25]).

The aquifer’s initial tracer concentration is assumed to be zero before starting the test: Cðr; h; t ¼ 0Þ ¼ 0

ð7Þ

The boundary, which is used to describe solute transport at the interface between the extraction well and aquifer while neglecting the finite volume in the well bore is as oCðr; h; tÞ ¼ 0 at r ¼ rw or

ð8Þ

The solute transport at the injection well boundary is considered as the other boundary condition. It should be noted that the velocity field near injection field will be distorted by the mere presence of the injection well (Fig. 2). The maximal width of the discharge near the injection well is w = 2arI, where a is a factor that defines the distortion of distance between the two most separated streaming line entering (or leaving) the injection well (this parameter can also depend on skin effect for an injection well). This small zone with advection-dominated flow has an aperture angle of 2d (Fig. 2). The aperture angle of this narrow zone at the distance r*  rL  l (l  rL) from the center of extraction is (see [19]; Eq. (3)) 2d ¼

2arI rL

ð9Þ

Because the concentration is symmetrical across h = 0 and h = p, only upper half domain is considered for the sake of mathematical simplicity. Therefore the boundary condition at injection well can be formulated as [19]  Cðr; h; tÞ ¼

C I ðtÞ

pd
0

0
ð10Þ

J.-S. Chen / Advances in Water Resources 30 (2007) 430–438

where C ID ¼ CCRI

Table 1 Dimensionless parameters used in this study

 C ID ðtD Þ ¼ lI

Dimensionless quantity

Expression

Time

tD ¼ pb/ðrQt2 r2 Þ

Distance

rD ¼ rrL

Extraction well radius

rwD ¼ rrwL

Concentrationa

C D ¼ CCR

Dispersivity ratio

X s ¼ aaTL hI lI ¼ /h rðrI r2Lr 2Þ I L w

Injection well mixing factor a

433

L

w

dC ID ðtD Þ dtD

ð19Þ

C ID ðtD ¼ 0Þ ¼ C 0 =C R

ð20Þ

oC D ðr; h; tÞ ¼0 oh

ð21Þ

at h ¼ 0

oC D ðrD ; h; tD Þ ¼0 oh

at h ¼ p

ð22Þ

Here C R ¼ /pbðrM2 0r2 Þ. L

W

3. Power series solution where CI(t) is the concentration generated in the injection well and transported downstream the narrow and short (a few well diameters) discharge zone by advection. The effluent concentration from the injection well in an ambient horizontal flow is yet to be determined. The effluent concentration from the well with initial dissolved tracer mass, M0, satisfies a mass balance equation for the tracer in a borehole:  2arI /bjV ðrL ÞjC I ðtÞ ¼ pr2I hI C I ðt ¼ 0Þ ¼

dC I ðtÞ dt

M0 ¼ C0 pr2I hI

ð11Þ

It is difficult to derive the analytical solution to Eq. (15) with variable coefficients using a special function. However, the initial and boundary value problem (15)–(22) can be approached analytically by using the power series method after successively applying the Laplace transform with respect to tD and finite Fourier cosine transform with respect to h. First, by performing the Laplace transform in Eqs. (15), (16) and its associated boundary conditions (17)– (22) with respect to tD yields

ð12Þ

eL ð1  rD Þ

where hI is the mixing length of the injection well ([19]; Eq. (6)). After integration, the known effluent concentration CI(t) can be substituted into the boundary condition (10). The physics of the problem stipulates that C is a singlevalued function in r and h coordinates. In addition, C is a continuous and symmetrical function across h = 0 and h = p. Thus, boundary conditions in the transverse directions are as follows: oCðr; 0; tÞ ¼0 oh oCðr; p; tÞ ¼0 oh

ð13Þ ð14Þ

Substituting the dimensionless variables given in Table 1 into Eq. (7), gives a dimensionless, scale-dependent advection–dispersion equation in the following form: o2 C D oC D eT ð1  rD Þ o2 C D þ ð1  e Þ þ L r2D orD o2 r D oh2 2rD R oC D ¼ ð1  r2wD Þ otD

eL ð1  rD Þ

ð15Þ

Consequently, the initial and boundary conditions (7)–(14) become C D ðrD ; h; 0Þ ¼ 0 oC D ðrD ; h; tD Þ ¼ 0 at rD ¼ rwD orD  C ID ðtD Þ p  d < h < p C D ðrD ; h; tD Þ ¼ 0 0
ð16Þ ð17Þ ð18Þ

¼

o2 G oG eT ð1  rD Þ o2 G þ ð1  eL Þ þ 2 orD r2D o rD oh2

2rD R sG ð1  r2wD Þ

ð23Þ

oGðrD ; h; sÞ ¼ 0 at rD ¼ rwD orD ( GID ðsÞ p  d < h < p GðrD ; h; sÞ ¼ 0 0
ð24Þ ð25Þ

 GID ðsÞ ¼ lI ½sGID ðsÞ  C 0 

ð26Þ

oGðrD ; h; sÞ ¼0 oh

at h ¼ 0

ð27Þ

oGðrD ; h; sÞ ¼0 oh

at h ¼ p

ð28Þ

where s denotes the Laplace transform parameters and G represents the Laplace transform of CD, as defined by Z 1 GðrD ; h; sÞ ¼ C D ðrD ; h; tD ÞestD dtD Z 1 0 GID ðsÞ ¼ C ID ðtD ÞestD dtD

ð29Þ ð30Þ

0

The application of the finite Fourier cosine transform with respect to h of Eqs. (23)–(28) gives: eL ð1  rD Þ 

d2 W dW þ ð1  eL Þ 2 dr d rD D

 eT ð1  rD Þ 2 2rD R s W ¼0  n þ r2D ð1  r2wD Þ

ð31Þ

434

J.-S. Chen / Advances in Water Resources 30 (2007) 430–438

dW ðrD ; n; sÞ ¼ 0 at rD ¼ rwD drD W ðrD ; n; sÞ ¼ F ðnÞ

ð32Þ ð33Þ

into Eq. (36) yields the following equation 1 1 X X m  ðm  1Þ  am rm2 þ ð1  eL Þ m  am  rm1 eL ð1  rD Þ D D m¼2

where n denotes the finite Fourier cosine transform parameter and W represents the finite Fourier cosine transform of G, as defined by W ðrD ; n; sÞ ¼ ( F ðnÞ ¼

Z

p

GðrD ; h; sÞ cosðnhÞ  dh

ð34Þ

0

GID ðsÞd n¼0 l C0 h n i : and GID ðsÞ ¼ I ð1Þ sinðndÞ l GID ðsÞ n ¼ 1; 2; 3 Is þ 1 n



m¼1 1 X

2Rs  rD  am rmD ¼ 0 ð1  r2wD Þ m¼0

ð42Þ

Rearranging the terms and shifting the summation indices, we obtain 1 1 X X ðm þ 2Þ  ðm þ 1Þ  amþ2  rmD  eL  ðm þ 1Þ  m eL  m¼0

m¼1

 amþ1 ðsÞ 

rmD

1 X þ ð1  eL Þ ðm þ 1Þ  amþ1 ðsÞ  rmD m¼0

Such a transform is advantageous in that the inversion is directly given by the following formula [21]: 1 1 2X GðrD ; h; sÞ ¼ W ðrD ; 0; sÞ þ W ðrD ; n; sÞ cosðnhÞ p p n¼1

ð35Þ

Governing equation (31) is an ordinary differential equation with variable coefficients and can be solved by directly applying the power series method. There are two different types of power series solutions for Eq. (31): Case 1, in which the coefficients are analytical for n = 0; Case 2, in which the coefficients are analytical for n 5 0. Therefore the analytical solution to the governing equation (31) of Case 1 is obtained by the application of the direct power series method. The extended power series method, the Frobenius method, is applied to yield the analytical solution to Eq. (31) of Case 2 [12]. The detailed procedures of the two approaches are described below: Case 1, n = 0, the coefficients are analytical. Rewriting Eq. (31) for n = 0 gives

1 X 2Rs   am1 ðsÞ  rmD ¼ 0 ð1  r2wD Þ m¼1

ð43Þ

Since Eq. (43) is an identity in rD, the sum of the coefficients of each power of rD in left hand side of Eq. (43) is set to zero as the identity property is invoked. For m = 0, this yields a2 ¼ 

ð1  eL Þa1 2eL

ð44Þ

and in general when m = 1, 2, 3, . . . amþ2 ¼

1 eL ðm þ 2Þðm þ 1Þ   ðm þ 1Þð1  eL  eL mÞamþ1 þ

 2Rs am1 : 1  r2wD ð45Þ

Substituting the coefficients correlation obtained from of Eqs. (44) and (45) into Eq. (39) gives the following general solution W ðrD ; 0; sÞ ¼ b1 Z 1 ðrD ; 0; sÞ þ b2 Z 2 ðrD ; 0; sÞ

d2 W dW 2Rs eL ð1  rD Þ 2 þ ð1  eL Þ rD W ¼ 0  drD drD ð1  r2wD Þ dW ðrD ; 0; sÞ ¼ 0 at rD ¼ rwD drD W ðrD ; 0; sÞ ¼ GID ðsÞd at rD ¼ 1

ð36Þ ð37Þ ð38Þ

The solution of Eq. (36) is assumed to be in the form of power series with unknown coefficients as W ðrD ; 0; sÞ ¼

1 X

am  rmD :

ð39Þ

m¼0

Substituting the power series of Eq. (39) along with its term-wise differentiation dW ðrD ; 0; sÞ ¼ drD

1 X

m1 m  am  r D

ð40Þ

m¼1

1 d2 W ðrD ; 0; sÞ X ¼ m  ðm  1Þ  am  rm2 D dr2D m¼2

ð41Þ

ð46Þ

where W1(rD, 0, s) and W2(rD, 0, s) are two linearly independent functions which are the power series form of (39) with coefficients determined by Eqs. (44) and (45). Case 2, n 5 0, the coefficients are not analytical. The coefficients of the transformed equation (36), except for n = 0 are not analytical. Therefore, the Frobenius power series method is applied to derive the analytical solution to Eq. (36). The solution of Eq. (36) is assumed to be in the form of the extended power series with undetermined coefficients as 1 X W ðrD ; n; sÞ ¼ am  rmþr ð47Þ D : m¼0

Insert this power series and its term-wise differentiation, 1 dW ðrD ; n; sÞ X mþr1 ¼ ðm þ rÞ  am  rD ð48Þ drD m¼0 1 d2 W ðrD ; n; sÞ X mþr2 ¼ ðm þ rÞ  ðm þ r  1Þ  am  rD dr2D m¼0

ð49Þ

J.-S. Chen / Advances in Water Resources 30 (2007) 430–438

W ðrD ; n; sÞ ¼ b1 Z 1 ðrD ; n; sÞ þ b2 Z 2 ðrD ; n; sÞ

into Eq. (36). This yields 1 X mþr2 ðm þ rÞðm þ r  1Þam rD eL ð1  rD Þ m¼0

þ ð1  eL Þ þ

1 X

mþr1 ðm þ rÞam rD 

m¼0 1 X

eT ð1  rD Þ 2 n r2D

2Rs rD am rmþr ¼0 D 1  r2wD m¼0

ð50Þ

m¼0 mþr  ðm þ r  2Þam1 rD þ ð1  eL Þ

P 11 ¼ W 1 ðrwD ; n; sÞ

mþr ðm þ r  1Þam1 rD

m¼1

 eT n2

1 X

2 am rmþr D þ eT n

m¼0

1 X

mþr am1 rD

m¼1

1 2Rs X  am3 rmþr ¼0 D 1  r2wD m¼3

ð51Þ

Equating the sum of the coefficient of each power of rD in the left hand side of Eq. (51) to zero yields a general formula for all m. eL rðr  1Þa0  eT n2 a0 ¼ 0

ðm ¼ 0Þ

ð52Þ

eL ðr þ 1Þra1  eL rðr  1Þa0 þ ð1  eL Þra0  eT n2 a1 þ eT n2 a0 ¼ 0 ðm ¼ 1Þ eL ðr þ 2Þðr þ 1Þa2  eL ðr þ 1Þra1

ð53Þ

þ ð1  eL Þðr þ 1Þa1  eT n2 a2 þ eT n2 a1 ¼ 0 ðm ¼ 2Þ eL ðm þ rÞðm þ r  1Þam  eL ðm þ r  1Þ

ð54Þ

 ðm þ r  2Þam1 þ ð1  eL Þ  ðm þ r  1Þam1  eT n2 am þ eT n2 am1 2Rs  am3 ¼ 0 ðm ¼ 3; 4; 5; . . .Þ 1  r2wD

F ðnÞ½P 22 W 1 ðrD ; n; sÞ þ P 12 W 2 ðrD ; n; sÞ P 12 Q21  P 22 Q11 ð61Þ

m¼1 1 X

ð60Þ

where Z1(rD, n, s) and Z2(rD, n, s) are two linearly independent functions that have the form of Eq. (47) with coefficients determined by (57)–(59) and set r = r1 and r = r2. The particular solution is obtained straightforwardly by imposing boundary conditions (32) and (33) on the general solution specified as (17). Therefore, this particular solution may be expressed as W ðrD ; n; sÞ ¼

By shifting summation indices, we obtain 1 1 X X mþr ðm þ rÞðm þ r  1Þam rD  eL ðm þ r  1Þ eL

435

P 21 ¼ W 2 ðrwD ; n; sÞ dW 1 ðrwD ; n; sÞ P 12 ¼ drD dW 2 ðrwD ; n; sÞ P 22 ¼ drD Q11 ¼ W 1 ð1; n; sÞ Q21 ¼ W 2 ð1; n; sÞ dW 1 ð1; n; sÞ drD dW 2 ð1; n; sÞ Q22 ¼ drD Q12 ¼

Solutions in the original domain, C(rD, h, tD), are the Laplace and finite Fourier cosine inverse transform of W(rD, n, s). The finite Fourier cosine inversion is performed using Eq. (35). An algorithm for the numerical inversion of the Laplace transform, such as De Hoog et al. [7], is required to execute the Laplace inverse transform; the robust nature of the De Hoog et al. [7] algorithm was demonstrated by Moench [17]. 4. Results and discussion

ð55Þ

From Eq. (52), the indicial equation is eL rðr  1Þ  eT n2 ¼ 0 pffiffiffiffiffiffiffiffiffiffiffiffi ð56Þ pffiffiffiffiffiffiffiffiffiffiffi ffi e e 1þ 1þ4n2 eT 1 1þ4n2 eT L L The roots are r1 ¼ and, r2 ¼ , respec2 2 tively. Eqs. (53)–(55) yield

It is necessary to verify the developed power series solution prior to its application. To date, no analytical or semi-analytical solutions expressed as special functions are available for verification. Therefore, a comparison with a numerical solution must be made to demonstrate the accuracy of the developed power series solution. In this study,

½eL rðr  1Þ  ð1  eL Þr  eT n2  a0 ½eT ðr þ 1Þr  eT n2  ½eL ðr þ 1Þr  ð1  eL Þðr þ 1Þ  eT n2  a1 a2 ¼ ½er ðr þ 2Þðr þ 1Þ  eT n2  n o 2Rs ½eL ðm þ r  1Þðm þ r  2Þ  ð1  eL Þðm þ r  1Þ  eT n2 am1 þ 1r 2 am3 wD am ¼ ½eL ðm þ rÞðm þ r  1Þ  eT n2 

a1 ¼

By inserting the coefficients correlation of Eqs. (57)–(59) into Eq. (47) the general solution is obtained as

ð57Þ ð58Þ

ð59Þ

the developed power series solution is compared with the numerical solution derived using the Laplace transform

J.-S. Chen / Advances in Water Resources 30 (2007) 430–438

Table 2 Descriptive simulation conditions and transport parameters in the model verification Parameter

Tracer test

Pumping rate (Q), m3 min1 Aquifer thickness (h), m Effective porosity (/), dimensionless Radius of extraction well (rw), m Extraction well mixing length (hw), m Radius of injection well (rI), m Injection well mixing length (hI), m Distance to the injection well (rL), m Injected mass (M), kg Longitudinal dispersivity/distance ratio (eL), dimensionless Transverse dispersivity/distance ratio (eT), dimensionless Dimensionless dispersivity ratio (X), dimensionless

2 10 0.2 0.05 10 0.05 10 5 10 1, 0.2, 0.05 0.2, 0.04, 0.01 0.2

0

50

100

150

2

3 power series solution numerical solution eL=0.05

1.6

Concentration [kg/m3]

finite difference (LTFD) technique [5]. The LTFD method was first developed by Moridis and Reddell [16] to solve partial differential equation of transient flow through porous media. In LTFD approach, the Laplace transform is performed first on the partial differential equation to eliminate the temporal derivative. The transformed equation is then approximated using a finite difference (FD) method to discretize the space variable. Unlike the LTPS solution that is analytical and continuous in both time and space, the LTFD solution is only exact in time and can be yield solution only at discrete nodal points. Additionally, the computation of LTFD solution is time-consuming. A hypothetical tracer test is considered for validation. The shapes of breakthrough curves are affected by physical parameters and they generally do not permit evaluation of the analytical solution for the breakthrough curve with a steep front that results from advective-dominated transport or the observation point located near the injection well. In this study the developed power series solution are verified with LTFD solution over a wide range of physical parameters of interest. Table 2 summarizes the simulation conditions and transport parameters for verification. Figs. 2–5 compare the breakthrough curves at the observation point with different radial distance values (r ¼ 15 rL , 25 rL , 35 rL and 45 rL ) and fixed transverse angle (h = p) obtained using the proposed power series solution and the LTFD solution for various dispersivity/distance ratios, eL = 1, 0.2 and 0.05. As expected, the power series solution agrees closely with the LTFD solution, thus validating the power series solution. The comparisons of the two solutions reveal that the power series technique can effectively and accurately solve the twodimensional, scale-dependent advection–dispersion equation with the variable-dependent coefficient. The developed power series solution is useful to generate the theoretical breakthrough curves at any observation point with small computational effort and helps to fast evaluate the longitudinal and transverse dispersivity/distance ratios by curvingmatching of theoretical breakthrough curves against breakthrough curves from a field tracer test.

2 1.2

0.8 1 eL=0.2 0.4 eL=1 0

0 0

50

100

150

Time [min]

Fig. 3. Comparison of the breakthrough curves obtained using the proposed power series solution and the numerical solution at the observation point ðr ¼ 15 rL ; h ¼ pÞ for eL = 1, 0.2 and 0.05.

0

50

4

100 eL=0.05

150 3

power series solution numerical solution

3 Concentration [kg/m3]

436

2

2

1 eL=0.2 1

eL=1

0

0 0

50

Time [min]

100

150

Fig. 4. Comparison of the breakthrough curves obtained using the proposed power series solution and the numerical solution at the observation point ðr ¼ 25 rL ; h ¼ pÞ for eL = 1, 0.2 and 0.05.

5. Summary and discussion In this paper, analytical solutions have been developed for a model investigating non-axisymmetrical contaminant

J.-S. Chen / Advances in Water Resources 30 (2007) 430–438 0

50

100

150 3

10 power series solution numerical solution

Concentration [kg/m3]

eL=0.05

2

437

excellent agreement between the power series solution and the numerical solution. The solutions obtained may be useful to estimate the scale-related transport parameters in longitudinal and transverse directions on the basis of the curve-fitting method (Fig. 6). Acknowledgements The author would like to thank the National Science Council of the Republic of China for financially supporting this work under Contract Nos. NSC 92-2313-B-008-003 and NSC 93-2313-B-008-002.

5

1 eL=0.2

References eL=1

0

0 0

50

Time [min]

100

150

Fig. 5. Comparison of the breakthrough curves obtained using the proposed power series solution and the numerical solution at the observation point ðr ¼ 35 rL ; h ¼ pÞ for eL = 1, 0.2 and 0.05.

transport in a radially convergent flow field with scaledependent dispersion. The analytical solution is derived by application of the power series method coupling with the Laplace and finite Fourier cosine transform. The power series method appears to be capable of effectively solving the two-dimensional, scale-dependent advection–dispersion equation with the variable-dependent coefficient due to

0

50

100

150

30

3 power series solution numerical solution

eL=0.05

Concentration [kg/m3]

20

2

eL=0.2

10

1

eL=1

0

0 0

50

100

150

Time [min]

Fig. 6. Comparison of the breakthrough curves obtained using the proposed power series solution and the numerical solution at the observation point ðr ¼ 45 rL ; h ¼ pÞ for eL = 1, 0.2 and 0.05.

[1] Aral MM, Liao B. Analytical solutions for two-dimensional transport equation with time-dependent dispersion coefficients. J Hydrol Eng 1996;1:20–32. [2] Becker MW, Charbeneau RJ. First-passage-time transfer functions for groundwater tracer tests conducted in radially convergent flow. J Contam Hydrol 2000;40:299–310. [3] Chen CS. Analytical and approximate solutions to radial dispersion from an injection well to a geological unit with simultaneous diffusion into adjacent strata. Water Resour Res 1985;21(8):1069–76. [4] Chen CS. Analytical solution for radial dispersion problem with Cauchy boundary at injection well. Water Resour Res 1987;23(7):1217–24. [5] Chen JS, Liu CW, Liang CP. Evaluation of longitudinal and transverse dispersivities/distance ratios for tracer test in a radially convergent flow field with scale-dependent dispersion. Adv Water Resour 2005;29(3):887–98. [6] Chen JS, Liu CW, Hsu HT, Liao CM. A Laplace transformed power series solution for solute transport in a convergent flow field with scale-dependent dispersion. Water Resour Res 2003;39(8). doi:10.1029/2003WR002299. [7] De Hoog FR, Knight JH, Stokes AN. An improved method of Laplace transforms using a Fourier series approximation. SIAM J Sci Stat Comput 1982;3(3):357–66. [8] Huang K, van Genuchten MT, Zhang R. Exact solutions for onedimensional transport with asymptotic scale-dependent dispersion. Appl Math Modell 1996;20(4):297–308. [9] Hunt B. Contaminant source solutions with scale-dependent dispersivity. J Hydrol Eng 1998;3(4):268–75. [10] Hunt B. Scale-dependent dispersion from a pit. J Hydrol Eng 2002;7(4):168–74. [11] Kocabas I, Islam MR. Concentration and temperature transients in heterogeneous porous media. Part II: Radial transport. J Petrol Sci Eng 2000;26:221–33. [12] Kreyszig E. Advanced engineering mathematics. New York: John Wiley and Sons; 1968. [13] Logan JD. Solute transport in porous media with scale-dependent dispersion and periodic boundary conditions. J Hydrol 1996;184:261–76. [14] Mishra S, Parker JC. Analysis of solute transport with a hyperbolic scale dependent dispersion model. Hydrol Proc 1989;4(1):45–57. [15] Moench A. Convergent radial dispersion: a Laplace transform solution for aquifer testing. Water Resour Res 1989;25(3):439–47. [16] Moridis GJ, Reddell DL. The Laplace transform finite difference method for simulation of flow through porous media. Water Resour Res 1991;27(80):1873–84. [17] Moench AF. Convergent radial dispersion: a note on evaluation of transform solution. Water Resour Res 1991;27(12):3261–4. [18] Pang L, Hunt B. Solutions and verification of a scale-dependent dispersion model. J Contam Hydrol 2001;53:21–39.

438

J.-S. Chen / Advances in Water Resources 30 (2007) 430–438

[19] Pickens JF, Grisak GE. Scale-dependent dispersion in a stratified granular aquifer. Water Resour Res 1981;17(4):1191–211. [20] Pickens JF, Grisak GE. Modelling of scale-dependent dispersion in a hydrogeological system. Water Resour Res 1981;17(6):1701–11. [21] Sneddon IN. The use of integral transform. New York: McGrawHill; 1972. [22] Su N, Sander GC, Liu F, Anh V, Barry DA. Similarity solutions for solute transport in fractal porous media using a time- and scaledependent dispersivity. Appl Math Modell 2005;29(9):852–70.

[23] Yates SR. An analytical solution for one-dimensional transport in heterogeneous porous media. Water Resour Res 1990;26(10): 2331–8. [24] Yates SR. An analytical solution for one-dimensional transport in heterogeneous porous media with an exponential dispersion function. Water Resour Res 1992;28(18):2149–54. [25] Zlotnik VA, Logan JD. Boundary conditions for convergent radial tracer tests and the effect of well bore mixing volume. Water Resour Res 1996;32(7):2323–8.