A mathematical solution and analysis of contaminant transport in a radial two-zone confined aquifer

(1)

A mathematical solution and analysis of contaminant transport in a radial

two-zone con

ﬁned aquifer

Yen-Ju Chen, Hund-Der Yeh

⁎

, Kai-Ju Chang

Institute of Environmental Engineering, National Chiao Tung University, 1001 University Road, Hsinchu 300, Taiwan

a r t i c l e i n f o

a b s t r a c t

Article history: Received 19 July 2011

Received in revised form 17 June 2012 Accepted 22 June 2012

Available online 29 June 2012

An aquifer with a wellbore surrounded by a finite-thickness skin, such as a gravel pack, can be regarded as a radial two-zone system. In this study, a mathematical model is developed to describe contaminant transport in a radial two-zone confined aquifer system. The solution of the model equations can be used to delineate contaminant transport in such a two-zone aquifer system and to investigate the skin effect on the temporal and spatial concentration distribution. The present solution is shown to reduce to an existing solution for transport in a homogeneous aquifer in the case when there is no wellbore skin. The results predicted from the solution indicate that the skin effect has a significant impact on the concentration distribution at early time. But an abrupt change in the spatial concentration distribution may occur at the interface of the skin zone and aquifer formation zone. The results from a sensitivity analysis reveal that dispersivity in the formation zone has a more significant effect on the concentration distribution than do the effects of skin thickness and dispersivity in the skin zone.

Keywords: Radial dispersion Contaminant transport Skin effect Laplace transform Groundwater 1. Introduction

A skin zone that occurs around a well due to drilling practices or well completion may have properties that are different from the aquifer formation. This situation can be considered as a radial two-zone aquifer system with a finite-thickness skin adjacent to the aquifer region. Numer-ous studies have considered the influence of a skin zone on

groundwater flow induced by various aquifer tests (e.g.,Wen

et al., 2011; Yang and Yeh, 2005, 2006), as well as on the

estimation of hydrogeological properties (e.g.,Chen and Yeh,

2009; Kabala, 2001; Yeh and Chen, 2007). It has been widely acknowledged that the existence of a skin zone affects the early-time head data near the well. However, little is known about the influence of a skin zone on contaminant transport in a radial two-zone aquifer system.

The radial contaminant transport problem represents the injection of a dissolved contaminant (e.g., waste water or treated

water) into an aquifer, and several articles have examined this

radial dispersion problem. Ogata (1958) first presented an

analytical solution for the radial advection–dispersion equation

(RADE) based on the Laplace transform and complex integral methods, though that study did not evaluate the integral

solution.Hoopes and Harleman (1967)used a finite difference

method to simulate solute transport in radial flow induced by a

line source injection. Subsequently, both Dagan (1971) and

Gelhar and Collins (1971)used perturbation methods to obtain approximate solutions for radial dispersion due to injection in a

finite-radius well. For the radial dispersion problem,Tang and

Babu (1979)developed an analytical solution expressed in terms of three integrals, including the combination of Bessel functions and/or modified Bessel functions, and an approximate solution in terms of the complementary error function. In addition, they also presented small-time and large-time solutions for the radial

dispersion problem.Moench and Ogata (1981)solved the RADE

by the Laplace transform technique. Their Laplace-domain solution was expressed in terms of Airy functions and

in-verted numerically to the time domain by theStehfest (1970).

Hsieh (1986)presented an analytical solution in terms of Airy ⁎ Corresponding author. Tel.: +886 3 5731910; fax: +886 3 5725958.

E-mail address:[email protected](H.-D. Yeh).

Contents lists available atSciVerse ScienceDirect

Journal of Contaminant Hydrology

j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / j c o n h y d

(2)

functions for the radial dispersion problem, which was simpler

than that ofTang and Babu (1979). These analytical solutions

and approximate solutions serve as useful and efficient tools to predict solute concentration distributions in homogeneous

aquifer systems. Recently, Zhan et al. (2009a, 2009b)dealt

with solute transport in a layered system with two-dimensional transport in the aquifer and one-dimensional transport in the aquitard. The concentration distributions in the aquifer and aquitard were solved based on the Laplace transform technique

along with the continuity requirements at the aquifer–aquitard

interface. The Stehfest numerical inversion algorithm (1970) was used to obtain time-domain results. Though the presence of a wellbore skin is expected to affect the concentration distribu-tions in the aquifer system, the dispersion problem in a radial two-zone aquifer, which accounts for the finite-thickness well-bore skin, has rarely been considered in the field of groundwater transport.

The purpose of this study is to develop a new mathematical solution for solute transport in a radial two-zone confined aquifer system. The case considered is of an injection well of finite radius that fully penetrates the aquifer. The Laplace-domain solution of the model is developed by the Laplace transform technique. Due to its complicated form, the time-domain solution is obtained by applying a numerical inversion. This solution is then used to investigate the effect of wellbore skin properties on the temporal and spatial concentration distribution in the aquifer system. In addition, the influence of dispersivities in the skin and formation zones as well as the skin thickness on predicted concentration values is assessed through a sensitivity analysis.

2. Methodology

2.1. Mathematical model and Laplace-domain solution

Fig. 1shows the configuration considered; a well in a radial two-zone aquifer system. The injection well is located at the origin of the cylindrical coordinate system, which has a finite radius, rw, and a screen that fully penetrates the aquifer. The contaminant is considered to be nonreactive and continuously released from the injection well with a constant concentration

C0. The confined aquifer has a finite thickness b and consists of

two zones in the radial direction. The first zone, known as the wellbore-skin zone, is located around the injection well and

has a thickness of r1–rw. The longitudinal dispersivity in the

first zone is denoted byα1. As indicated inFig. 1, r1is the radial

distance from the central line of the well to the outer radius of the wellbore skin. The second zone is known as the aquifer formation zone, which is homogeneous, isotropic, and of infinite lateral extent with a longitudinal dispersivity denoted

byα2. The contaminant concentrations in the first and second

zones are denoted as C1and C2, respectively.

For convenience, the mathematical model is expressed in dimensionless form. The contaminant transport is dominated by radial advection and longitudinal dispersion in the two-zone system. The molecular diffusion coefficient is assumed to be negligible. The dimensionless concentration distribution in the skin zone is written as

κ∂ 2 G1 ∂ρ2 − ∂G1 ∂ρ ¼ ρ∂G 1 ∂τ ; ρwb ρ ≤ ρ1; τ > 0 ð1Þ

and in the aquifer formation zone as

∂2 G2 ∂ρ2 − ∂G2 ∂ρ ¼ ρ∂G 2 ∂τ ; ρ1< ρ < ∞; τ > 0 ð2Þ

where G1=C1/C0and G2= C2/C0are the dimensionless

con-centrations in the skin and aquifer formation zones, respectively;

κ=α1/α2represents the ratio of longitudinal dispersivities in

the skin and formation zones;ρ=r/α2represents the

dimen-sionless radial distance from the well;ρw=rw/α2andρ1=r1/α2

are the dimensionless well radius and outer radius of the skin

zone, respectively;τ=Qt/(2πbnα22_{) is the dimensionless time}

since injection of contaminant, and Q is the constant well injection rate, and n is the aquifer porosity. Both the skin and aquifer formation zones are initially free from contamination, expressed as

G1ðρ; 0Þ ¼ G2ðρ; 0Þ ¼ 0; ρ ≥ ρw: ð3Þ

The boundary condition at the injection well is a constant concentration denoted as

G1ðρw; τÞ ¼ 1; τ > 0: ð4Þ

The aquifer formation is assumed to be free from contamination far from the well, where the concentration is written as

G2ð∞; τÞ ¼ 0; τ > 0 ð5Þ

The concentration and the flux at the interface between the skin and aquifer formation zones are continuous, expressed as

G1ðρ1; τÞ ¼ G2ðρ1; τÞ; τ > 0 ð6Þ and κ∂G1ðρ1; τÞ ∂ρ ¼∂G 2ðρ1; τÞ ∂ρ ; τ > 0: ð7Þ

Taking the Laplace transform with respect to the

dimen-sionless time_{τ of Eqs.}(1) to (7)gives

κd 2_G 1 dρ2 − d G1 dρ −ρsG1¼ 0; ρw< ρ ≤ ρ1 ð8Þ 1 α α2 w r 1 r ) , ( 1rt C C2(r,t) Dimensional form

Well Skin zone Formation zone

Injection rate, Q Input concentration, C0

(First zone) (Second zone) Concentrations Dispersivities Aquifer thickness, b 0 = r

Fig. 1. Schematic diagram of an injection well in a radial two-zone aquifer system.

(3)

d2G2 d_ρ2 − d G2 dρ −ρsG2¼ 0; ρ1< ρ < ∞ ð9Þ G1ðρw; sÞ ¼ 1 s ð10Þ G2ð∞; sÞ ¼ 0 ð11Þ G1ðρ1; sÞ ¼ G2ðρ1; sÞ ð12Þ and κd G1ðρ1; sÞ dρ ¼ d G2ðρ1; sÞ dρ ð13Þ

where G1 and G2 are the dimensionless Laplace-domain

concentrations in the skin and formation zones, respectively,

and s is the Laplace variable. The development of Eqs.(8) to

(13)is given inAppendix A. The results for G1and G2are,

respectively, given by G1ðρ; sÞ ¼ 1 sexp ρ−ρw 2κ ( Ai Z1;ρ Ai Z1;ρw 1−Bi Z1;ρw _∇ 1 ∇ þBi Z1;ρ _∇ 1 ∇ ) ð14Þ and G2ðρ; sÞ ¼ 1 sexp ρ1−ρw 2κ þρ−ρ2 1 Ai Z2;ρ Ai Z2;ρ1 × Ai Z1;ρ1 Ai Z_1;ρ w 1−Bi Z1;ρw _∇ 1 ∇ þ Bi Z1;ρ1 _∇ 1 ∇ 8 < : 9 = ; ð15Þ

where Ai(_{∙) and Bi(∙) are Airy functions. Other associated}

functions and variables are defined in Appendix A. The

argument Z1,ρis the abbreviation for the function Z1(ρ, s), as

defined in Eq.(A.5). In addition, Z1,ρ1and Z1,ρware referred to

as the function Z1by replacingρ with their second indexes ρ1

andρw, respectively. Similarly, Z2,ρand Z2,ρ1are the

abbrevi-ations for the functions Z2(ρ, s) and Z2(ρ1, s), respectively,

defined in Eq.(A.6).

2.2. Numerical evaluation

Numerical inversion and evaluation were carried out using a FORTRAN code with double precision accuracy. The Airy functions with argument x can be expressed in terms of

the modified Bessel functions as (Abramovitz and Stegun,

1972) Ai xð Þ ¼1 π ffiffiffi x 3 r K1=3ð Þξ ð16Þ and Bi x_{ð Þ ¼} ffiffiffi x 3 r I₋₁₌₃_{ð Þ þ I}_ξ 1=3ð Þξ h i ð17Þ

where I(_{∙) and K(∙) are the modified Bessel functions of}

the first and second kinds, respectively, with the argument

ξ=2/3⋅x3/2_{. The derivatives of the Airy functions, which}

appear in Eqs.(A.13) and (A.14), are expressed by

Ai′ð Þ ¼ −x 1_π xffiffiffi 3 p K2=3ð Þξ ð18Þ and Bi′ð Þ ¼x xffiffiffi 3 p I−2=3ð Þ þ Iξ 2=3ð Þξ h i : ð19Þ

The functions I(∙) and K(∙) are calculated by the

sub-routines DCBIS and DCBKS ofIMSL (2003a), respectively. The

numerical inversion of Eqs.(14) and (15)to the time domain

is achieved by the subroutine DINLAP ofIMSL (2003b)which

was developed based on the Crump algorithm (Chen et al.,

1996).

2.3. Reduction to Moench and Ogata's solution

Moench and Ogata (1981) provided a Laplace-domain solution for the problem of radial dispersion with the contaminant injected into a well in an aquifer without the presence of wellbore skin. Their Laplace-domain solution is expressed as G ρ; sð Þ ¼1_sexp ρ−ρw 2 Ai Yρ Ai Yð Þw ð20Þ where Y_ρ¼4ρs þ 1 4s2=3 ð21Þ and Yw¼ 4_ρwsþ 1 4s2=3 : ð22Þ

Note that the present solution, i.e., Eqs.(14) and (15), can

be reduced to Eq.(20)by settingρ1=ρwandκ=1, i.e., Z1=

Z2= Yρ, Z1,ρ1= Z2,ρ1= Z1,ρw= Yw, and∇1= 0.

2.4. Sensitivity analysis

Sensitivity analysis can be performed to examine the effect of changing the input parameters for model output, i.e., the concentration. The normalized sensitivity coefficients of

Uα1, Uα2and Ur1with respect to each of the parametersα1,α2

and r1indicate the relative importance of model parameters,

which are expressed as (Liou and Yeh, 1997)

U_α 1 ¼ α1 ∂C ∂α1 ; Uα2¼ α2 ∂C ∂α2 ; and Ur1¼ r1 ∂C ∂r1 : ð23Þ

Note that these three coefficients have the same di-mensions as concentration. Three new sensitivity coefficients are further introduced as

u_α 1¼ κ ∂G ∂κ¼ Uα1=C0 ð24Þ u_α 2¼ 1 ρw⋅ ∂G ∂ 1 ρw ¼ U_α₂=C0 ð25Þ

(4)

and

ur1¼ ρ1

∂G

∂ρ1

¼ Ur1=C0 ð26Þ

where uα1, uα2and ur1represent the dimensionless

normal-ized sensitivity coefficients of the dimensionless

concentra-tion with respect to the dimensionless parametersκ, 1/ρw

andρ1, respectively. They are approximated by the

finite-difference formulas as u_α 1≈κ Gðκ þ ΔκÞ−G κð Þ Δκ ð27Þ u_α 2≈α2 G½ρwðα2þΔα2Þ; ρ1ðα2þΔα2Þ; κ αð 2þΔα2Þ; τ αð 2þΔα2Þ−G ρ½ w; ρ1; κ; τ Δα2 ≈_ρ1 w G½ρwðα2þΔα2Þ; ρ1ðα2þΔα2Þ; κ αð2þΔα2Þ; τ αð 2þΔα2Þ−G ρ½ w; ρ1; κ; τ Δ _ρ1 w ð28Þ and ur1≈ρ1 Gðρ1þ Δρ1Þ−G ρð Þ1 Δρ1 ð29Þ

whereΔκ, Δα2andΔρ1denote small changes in the

param-eters, which are usually chosen as the parameter values

multiplied by a factor of 10−3. Thus, the following formulas

are applied to calculate Eq.(28); i.e.,ρw(α2+Δα2) =ρw/1.001,

ρ1(α2+Δα2) =ρ1/1.001, κ(α2+Δα2) =κ/1.001, and τ(α2+

Δα2) =τ/1.001.

3. Results and discussion

3.1. Effect of dispersivities in the formation zone

Fig. 2displays the curves of dimensionless concentration

versus dimensionless distance forρw= 1,ρ1= 2, andκ=0.5, 1

and 2 when_{τ=0.5, 4.5 and 18. The case of κ=1 represents the}

absence of the skin zone. Furthermore, the cases ofκ=0.5 and

2 represent the situations forα1bα2andα1>α2, respectively.

At the same distance, the highest dimensionless concentration

is produced in the case of greatestκ at a specific time. Overall,

the differences in dimensionless concentrations forκ=0.5, 1

and 2 decrease with time. The slopes of the dimensionless concentration curve are markedly different at the interface

between the skin zone and formation zone (ρ=2) in the cases

of κ=0.5 and 2 when τ=0.5 and 4.5. The variations in the

slope of the concentration curve may provide a diagnostic indication for the presence of the skin zone in the aquifer system, although it is recognized that such detailed data on concentration versus distance may rarely be available in the field. 0 0.2 0.4 0.6 0.8 1 τ = 0.5 1 2 3 4 5 6 7 8 9 10 τ = 4.5 τ = 18 κ = 0.5 κ = 1 κ = 2 τ = 0.5 τ = 4.5 τ = 18

Dimensionless distance,

ρ

Dimensionless concentration,

G = C

/C

0 ρw= 1, ρ1 = 2

Fig. 2. Dimensionless concentration distributions forρw= 1,ρ1= 2, and

κ=0.5, 1 and 2 when τ=0.5, 4.5 and 18.

1 2 3 4 5 6 7 8 9 10 0 0.2 0.4 0.6 0.8 1 τ = 0.5 τ = 4.5 τ = 18 κ = 0.5 κ = 1 κ = 2 τ = 0.5 τ = 4.5 τ = 18

Dimensionless distance,

ρ

Dimensionless concentration,

G = C

/C

0 ρ_w= 1, ρ₁ = 4

κ=0.5, 1 and 2 when τ=0.5, 4.5 and 18.

0 0.2 0.4 0.6 0.8 1 τ = 50 10 τ = 450 τ = 1800 κ = 0.5 κ = 1 κ = 2 τ = 50 τ = 450 τ = 1800 20 30 40 50 60 70 80 90 100

Dimensionless distance,

ρ

Dimensionless concentration,

G = C

/C

0 ρw= 10, ρ1 = 20

(5)

Fig. 3 shows the curves of dimensionless concentration

versus dimensionless distance forρw=10,ρ1=20, andκ=0.5,

1 and 2 whenτ=0.5, 450 and 1800, which is comparable with

the corresponding curves inFig. 2to represent the case of smaller

α2 under the same well radius condition. Additionally, the

dimensionless distance and time adopted inFig. 3are 10 and 100

times, respectively, those used inFig. 2for the same dimensional

values. InFig. 3, the differences in dimensionless concentrations

forκ=0.5, 1 and 2 are less significant than those inFig. 2. The

dimensionless concentration distributions forκ=0.5, 1 and 2

closely coincide when the time is large (τ=1800). In contrast to

Fig. 2, the highest dimensionless concentration does not always

occur in the case of greatestκ at a specific time.Fig. 3shows that

these dimensionless concentration curves cross near ρ=14

whenτ=50, near ρ=30 when τ=450, and near ρ=50 when

τ=1800.Fig. 3also indicates that the plumes migrate shorter

distances at corresponding times than those inFig. 2.

3.2. Effect of skin thickness

Fig. 4illustrates the dimensionless spatial concentration

distributions forρw= 1,ρ1= 4, andκ=0.5, 1 and 2 when

τ=0.5, 4.5 and 18 for the case of thicker skin. This figure indicates that the differences in dimensionless concentration

for κ=0.5, 1 and 2 are much more significant at larger

distances from the well than those inFig. 2. For the case of

κ=0.5, the dimensionless concentration does not change

significantly at the early time (τ=0.5) as the skin thickness

increases. In the case ofκ=0.5, the dimensionless

concen-tration forρ1= 4 (Fig. 4) is, however, much higher at short

distances and lower at long distances forτ=4.5 and 18 than

0 0.2 0.4 0.6 0.8 1 10 20 τ = 50 τ = 450 τ = 1800 κ = 0.5 κ = 1 κ = 2 τ = 50 τ = 450 τ = 1800 30 40 50 60 70 80 90 100

Dimensionless distance,

ρ

Dimensionless concentration,

G = C

/C

0 ρw= 10, ρ1 = 40

κ=0.5, 1 and 2 when τ=50, 450 and 1800.

1 2 3 4 5 6 7 8 9 10

Dimensionless distance,

ρ

1 2 3 4 5 6 7 8 9 10

Dimensionless distance,

ρ

1 2 3 4 5 6 7 8 9 10

Dimensionless distance,

ρ

-3 -2 -1 0 1 2

Dimensionless normalized sensitivity

,

u = U

/C

0 -3 -2 -1 0 1 2

Dimensionless normalized sensitivity

,

u = U

/C

0 -3 -2 -1 0 1 2

Dimensionless normalized sensitivity

,

u = U

/C

0 0 1 2 3 4 5

Dimensionless concentration, G = C/C

0 0 1 2 3 4 5

Dimensionless concentration, G = C/C

0 0 1 2 3 4 5

Dimensionless concentration, G = C/C

0

(a)

ρ

w

= 1, ρ

1

= 4, κ = 0.5

u_α1 u_α2 u_r₁ τ = 0.5 τ = 4.5 τ = 18 u_α1 u_α2 u_r₁ τ = 0.5 τ = 4.5 τ = 18 u_α1 u_α2 u_r₁ τ = 0.5 τ = 4.5 τ = 18

(b)

ρ

_w

= 1, ρ

₁

= 4, κ = 1

(c)

ρ

w

= 1, ρ

1

= 4, κ = 2

Fig. 6. Dimensionless concentration distributions and dimensionless nor-malized sensitivities (nornor-malized sensitivities per C0) of dimensionless

concentration in response to the change in the parametersα1,α2and r1for

(a)κ=0.5, (b) κ=1, and (c) κ=2 as well as ρw= 1 andρ1= 4 whenτ=0.5,

(6)

the case ofρ1= 2 (Fig. 2). In contrast, the dimensionless concentration is lower at short distances but higher at long

distances in the case ofκ=2 (thicker-skin) than those in

Fig. 2. Moreover, the spatial dimensionless concentration

distributions for the case of_{κ=2 have distinct changes of}

slope at the interface between the skin and formation zones

(ρ=4) when τ=4.5 and 18.

Fig. 5 shows the dimensionless spatial concentration distributions for the case that the skin thickness is twice that inFig. 3. The parameters used inFig. 5areρw= 10 andρ1= 40

forκ=0.5, 1 and 2 when τ=50, 450 and 1800. The increase in

skin thickness did not distinctly affect the dimensionless

concentration distribution at the early time (τ=50). The

differences in dimensionless concentration forκ=0.5, 1 and

2 are much more significant atτ=450 and 1800 than those in

Fig. 3. In comparison with Fig. 3, the same observation, as

discussed inFig. 4, applies to the case ofκ=0.5, where the

dimensionless concentration is much higher at short distances

and lower at long distances forτ=450 and 1800. In addition, a

lower dimensionless concentration is predicted at short distances but a higher concentration is predicted at long

distances in the case ofκ=2 as the skin thickness increases.

Figs. 4 and 5show that the spatial concentration distributions

change more sharply in the case ofκ=0.5.

3.3. Sensitivity analysis

Fig. 6illustrates the spatial dimensionless concentration distributions and the dimensionless normalized sensitivities of dimensionless concentration with respect to each of the

parametersα1,α2and r1forκ=0.5, 1 and 2 when ρw= 1 and

ρ1= 4. Generally, the values of dimensionless normalized sensitivity curves deviate from zero when the dimensionless concentration is between zero and one. The dimensionless

concentration in response to the change in α2 is more

sensitive than those inα1and r1. As shown inFig. 6, positive

perturbations inα2produce a negative concentration change.

Fig. 6b indicates that the dimensionless concentration is

insensitive to a change in r1when the skin is absent (κ=1).

In other words, the change of skin thickness in such a case

(κ=1) does not affect the dimensionless concentration since

the dispersivities of the skin and formation zones are the same. In addition, the dimensionless concentration seems

insensitive to the change in r1 when the plume has not

reached the formation zone, i.e.,τ=0.5, for κ=0.5 and 2, as

displayed inFig. 6a and c. These two figures also show that

significant decreases or increases of the sensitivity curves occur at the interface of skin zone and formation zone, i.e.,

ρ1= 4. Fig. 6a shows that a positive perturbation in r1

produces a positive concentration change within the skin

zone and a negative change beyond the skin zone whenτ=

4.5 and 18 forκ=0.5. In contrast, for the case of κ=2 shown

inFig. 6c, a positive perturbation in r1produces a negative change in the concentration within the skin zone and a

positive change beyond the skin zone whenτ=4.5 and 18.

4. Conclusions

A mathematical model has been developed to describe the temporal and spatial dimensionless concentration distributions due to the continuous injection of a solute in a radial two-zone

confined aquifer system. The Laplace-domain solution of the model equations was obtained in terms of Airy functions. A

numerical Laplace inversion, called DINLAP ofIMSL (2003b), was

applied to evaluate the dimensionless concentration in the time domain. The present Laplace-domain solution for contaminant

transport in a radial two-zone system can be reduced toMoench

and Ogata (1981)if the wellbore skin is absent.

The present solution was applied to investigate the skin effect on the temporal and spatial concentration distributions. The differences in concentrations between the skin-affected cases and non-skin case become less significant at longer time but are much more notable when the dispersivity in the aquifer formation zone is large. In addition, there may be a sudden change in the slope of the concentration distribution at the interface between the skin zone and formation zone. An aquifer with a smaller dispersivity in the aquifer formation zone was shown to produce a shorter travel distance of the plume. Results from the sensitivity analysis indicated that the dimensionless concentration was more sensitive to the change of dispersivity in the formation zone than those of two other parameters, i.e., skin thickness and dispersivity in the skin zone.

Notation

The following symbols are used in this paper:

Ai(⋅) Airy function

b aquifer thickness

b1 coefficient in Eq.(A.7), which is defined in Eq.(A.9)

b2 coefficient in Eq.(A.7), which is defined in Eq.(A.10)

Bi(⋅) Airy function

C0 initial constant concentration injected into the well

C1 contaminant concentration in the skin zone

C2 contaminant concentration in the formation zone

d1 coefficient in Eq.(A.8), which is defined in Eq.(A.11)

d2 coefficient in Eq.(A.8), which is defined in Eq.(A.12)

G1 dimensionless concentration in the skin zone,

which is defined as G1= C1/C0

G2 dimensionless concentration in the formation zone,

which is defined as G2= C2/C0

G1 dimensionless Laplace-domain concentration in

the skin zone

G2 dimensionless Laplace-domain concentration in

the formation zone

I(⋅) modified Bessel function of the first kind

K(⋅) modified Bessel function of the second kind

n aquifer porosity

Q constant injection rate at the well

r1 outer radius of the skin zone

rw radius of injection well

s Laplace variable

uα1 normalized sensitivity coefficient of dimensionless

concentration toκ, which is defined in Eq.(24)

uα2 normalized sensitivity coefficient of dimensionless

concentration to 1/ρw, which is defined in Eq.(25)

ur1 normalized sensitivity coefficient of dimensionless

concentration to_ρ1, which is defined in Eq.(26)

Uα1 normalized sensitivity coefficient of concentration

toα1, which is defined in Eq.(23)

Uα2 normalized sensitivity coefficient of concentration

(7)

Ur1 normalized sensitivity coefficient of concentration

to r1, which is defined in Eq.(23)

Ū1 function defined in Eq.(A.3)

Ū2 function defined in Eq.(A.4)

x argument

Yρ argument defined in Eq.(21)

Yw argument defined in Eq.(22)

Z1,ρ abbreviation of the function Z1(ρ, s) defined in

Eq.(A.5)

Z1,ρ1 abbreviation of the function Z1(ρ1, s), which can be

calculated from Eq.(A.5)

Z1,ρ_w abbreviation of the function Z1(ρw, s), which can be

Z2,ρ abbreviation from the function Z2(ρ, s) defined in

Eq.(A.6).

Z2,ρ1 abbreviation of the function Z2(ρ1, s), which can be

α1 longitudinal dispersivity in the skin zone

α2 longitudinal dispersivity in the formation zone

∇ equation defined in Eq.(A.14)

∇1 equation defined in Eq.(A.13)

κ ratio of longitudinal dispersivities in the skin and

formation zones, which is defined asκ=α1/α2

ρ dimensionless radial distance from the well, which

is defined asρ=r/α2

ρw dimensionless well radius, which is defined as

ρw= rw/α2

ρ1 dimensionless outer radius of the skin zone, which

is defined asρ1= r1/α2

τ dimensionless time since contaminant injection,

which is defined asτ=Qt/(2πbnα22₎

ξ argument defined asξ=2/3⋅x3/2

Acknowledgments

This study was partly supported by the Taiwan National Science Council under the grants NSC 99-2221-E-009-062-MY3, NSC 100-2221-E-009-106, and NSC 101-3113-E-007-008. The authors would like to thank the editor, Dr. Greg B. Davis, for editing the manuscript and two anonymous reviewers for their valuable and constructive comments that help improve the clarity of our presentation.

Appendix A. Development of Eqs.(14) and (15)

Eqs.(8) and (9)can be rearranged as the Airy differential equations: d2U1 dZ12 −Z1U1¼ 0; ρw< ρ ≤ ρ1 ðA:1Þ and d2U2 dZ22 −Z2U2¼ 0; ρ1< ρ < ∞ ðA:2Þ

where the new functionsŪ1andŪ2are respectively related to

G1and G2by U1ðρ; sÞ ¼ exp −₂ρ_κ G1ðρ; sÞ ðA:3Þ and U2ðρ; sÞ ¼ exp −ρ 2 G2ðρ; sÞ: ðA:4Þ

In addition, Z1and Z2are defined as

Z1ðρ; sÞ ¼ s κ 1=3 ρ þ₄1_κs ðA:5Þ and Z2ðρ; sÞ ¼ 4ρs þ 1 4s2=3 : ðA:6Þ

The general solutions to Eqs. (A.1) and (A.2) can be

respectively expressed as

U₁_ð_{ρ; s}_{Þ ¼ b}₁Ai Zð Þ þ b1 2Bi Zð Þ1 ðA:7Þ

and

U2ðρ; sÞ ¼ d1Ai Zð Þ þ d2 2Bi Zð Þ2 ðA:8Þ

where the coefficients b1, b2, d1 and d2 are the constants

determined by the associated boundary conditions(10) to

(13). The coefficients are then obtained as

b1¼ 1 sexp − ρw 2κ ₁ Ai Z1;ρw − ∇1 ∇ Bi Z1;ρw Ai Z1;ρw 2 4 3 5 _ðA:9Þ b2¼ 1 sexp − ρw 2κ _∇ 1 ∇ ðA:10Þ d1¼ 1 sexp ρ1−ρw 2κ −ρ 1 2 ₁ Ai Z_1;ρ w Ai Z_2;ρ 1 Ai Z1;ρ1 þ Ai Z1;ρw Bi Z_1;ρ 1 −Ai Z1;ρ1 Bi Z_1;ρ w h _{i ∇}₁ ∇ ðA:11Þ and d2¼ 0 ðA:12Þ where ∇1¼ κ 2=3_Ai′ _Z 1;ρ1 Ai Z2;ρ1 −Ai Z1;ρ1 Ai′ Z2;ρ1 ðA:13Þ and ∇ ¼ κ2=3_{Ai Z} 2;ρ1 Ai′ Z1;ρ1 Bi Z1;ρw −Ai Z1;ρw Bi′ Z1;ρ1 h i −Ai′ _Z 2;ρ1 Ai Z1;ρ1 Bi Z1;ρw −Ai Z1;ρw Bi Z1;ρ1 h i : ðA:14Þ

The notations Z1,ρw, Z1,ρ1and Z2,ρ1are the abbreviations for

the functions Z1(ρw, s), Z1(ρ1, s) and Z2(ρ1, s), respectively,

which can be calculated from Eqs. (A.5) and (A.6). The

Laplace-domain solutions of dimensionless concentrations in

skin and formation zones, Eqs. (14) and (15), are then

obtained by substituting Eqs.(A.7) and (A.8)into Eqs.(A.3)

(8)

References

Abramovitz, M., Stegun, I.A., 1972. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover Publications, New York.

Chen, Y.J., Yeh, H.D., 2009. Parameter estimation/sensitivity analysis for an aquifer test with skin effect. Ground Water 47 (2), 287–299. Chen, J.S., Liu, C.W., Chen, C.S., Yeh, H.D., 1996. A Laplace transform solution

for tracer tests in a radially convergent flow field with upstream dispersion. Journal of Hydrology 183 (3–4), 263–275.

Dagan, G., 1971. Perturbation solutions of the dispersion equation in porous mediums. Water Resources Research 7 (1), 132–142.

Gelhar, L.W., Collins, M.A., 1971. General analysis of longitudinal dispersion in nonuniform flow. Water Resources Research 7 (6), 1511–1521. Hoopes, J.A., Harleman, D.R.F., 1967. Dispersion in radial flow from a

recharge well. Journal of Geophysical Research 72 (14), 3595–3607. Hsieh, P.A., 1986. A new formula for the analytical solution of the radial

dispersion problem. Water Resources Research 22 (11), 1597–1605. IMSL, 2003a. IMSL Fortran Library User's Guide Math/Library Special

Function, Version 5.0. Visual Numerics, Houston, Texas.

IMSL, 2003b. IMSL Fortran Library User's Guide Math/Library Volume 2 of 2, Version 5.0. Visual Numerics, Houston, Texas.

Kabala, Z.J., 2001. Sensitivity analysis of a pumping test on a well with wellbore storage and skin. Advances in Water Resources 24 (5), 483–504. Liou, T.S., Yeh, H.D., 1997. Conditional expectation for evaluation of risk groundwater flow and solute transport: one-dimensional analysis. Journal of Hydrology 199 (3–4), 378–402.

Moench, A.F., Ogata, A., 1981. A numerical inversion of the Laplace transform solution to radial dispersion in a porous medium. Water Resources Research 17 (1), 250–252.

Ogata, A., 1958. Dispersion in porous media, Ph.D. dissertation, Northwest-ern University, Evanston, Illinois.

Stehfest, H., 1970. Numerical inversion of Laplace transforms. Communica-tions of the ACM 13 (1), 47–49.

Tang, D.H., Babu, D.K., 1979. Analytical solution of a velocity dependent dispersion problem. Water Resources Research 15 (6), 1471–1478. Wen, Z., Zhan, H.B., Huang, G.H., Jin, M.G., 2011. Constant-head test in a leaky

aquifer with a finite-thickness skin. Journal of Hydrology 399 (3–4), 326–334.

Yang, S.Y., Yeh, H.D., 2005. Laplace-domain solutions for radial two-zone flow equations under the conditions of constant-head and partially penetrat-ing well. Journal of Hydraulic Engineerpenetrat-ing ASCE 131 (3), 209–216. Yang, S.Y., Yeh, H.D., 2006. A novel analytical solution for constant-head test

in a patchy aquifer. International Journal for Numerical and Analytical Methods in Geomechanics 30, 1213–1230.

Yeh, H.D., Chen, Y.J., 2007. Determination of skin and aquifer parameters for a slug test with wellbore-skin effect. Journal of Hydrology 342 (3–4), 283–294.

Zhan, H.B., Wen, Z., Gao, G., 2009a. An analytical solution of two-dimensional reactive solute transport in an aquifer–aquitard system. Water Resources Research 45, W10501,http://dx.doi.org/10.1029/2008WR007479. Zhan, H.B., Wen, Z., Huang, G., Sun, D., 2009b. Analytical solution of

two-dimensional solute transport in an aquifer–aquitard system. Journal of Contaminant Hydrology 107, 162–174.