A new analytical solution solved by triple series equations method for constant-head tests in confined aquifers

(1)

A new analytical solution solved by triple series equations method for constant-head

tests in con

ﬁned aquifers

Ya-Chi Chang, Hund-Der Yeh

⁎

Institute of Environmental Engineering, National Chiao-Tung University, Hsinchu, Taiwan

a b s t r a c t

a r t i c l e i n f o

Article history: Received 21 May 2009

Received in revised form 18 March 2010 Accepted 18 March 2010

Available online 25 March 2010

Keywords: Constant-head test Aquifer

Mixed boundary value problem Partially penetrating well

The constant-head pumping tests are usually employed to determine the aquifer parameters and they can be performed in fully or partially penetrating wells. Generally, the Dirichlet condition is prescribed along the well screen and the Neumann type no-flow condition is specified over the unscreened part of the test well. The mathematical model describing the aquifer response to a constant-head test performed in a fully penetrating well can be easily solved by the conventional integral transform technique under the uniform Dirichlet-type condition along the rim of wellbore. However, the boundary condition for a test well with partial penetration should be considered as a mixed-type condition. This mixed boundary value problem in a confined aquifer system of infinite radial extent and finite vertical extent is solved by the Laplace and finite Fourier transforms in conjunction with the triple series equations method. This approach provides analytical results for the drawdown in a partially penetrating well for arbitrary location of the well screen in afinite thickness aquifer. The semi-analytical solutions are particularly useful for the practical applications from the computational point of view.

1. Introduction

Hydraulic parameters such as hydraulic conductivity, specific storage and leakage factor are important for quantifying groundwater resources. To determine these parameters, the constant-head pumping test is generally employed if the aquifer has low permeability. During the test, the hydraulic head at the well is kept constant throughout the test period and the transientflow rate across the wellbore is measured at the same time. A pumping test performed in a fully or partially penetrating borehole is influenced by the well skin effect. The fully penetrating well can be simulated as a Dirichlet (also called thefirst type) boundary condition, and the relative models can be solved by the conventional integral transform techniques[11]. If the well skin effect is negligible in the model, the Dirichlet and Neumann (second type) boundary conditions are appropriate for describing the drawdown along the well screen and casing, respectively. Thus, the boundary condition along the well face in the partially penetrating well is a mixed-type condition. The term“mixed-type” boundary condition is used to distinguish this boundary condition from the“uniform” Dirichlet and Neumann boundary condition or a combination of Dirichlet and Neumann boundary conditions (which is usually defined as third type or Robin type boundary condition). If the well skin effect is considered in the model, a more appropriate description for such an aquifer system is to treat the skin zone as a different formation zone instead of using a skin

factor. Thus, the aquifer system naturally becomes a two-zone formation (see, e.g.,[26,27,34,35]).

Mixed boundary conditions are widely used to describe many boundary value problems of mathematical physics. Such problems arise in potential theory and its numerous applications to engineering, fracture mechanics, heat conduction, and many others. Only limited analytical solutions to mixed boundary problems (MBPs) in theﬁeld of well hydraulics have been found so far by special solution techniques including the dual integral/series equation [8,22], Weiner–Hopf technique[18], and Green's function[12]. Most of the solutions to MBPs have been obtained numerically [29], or by approximate methods such as asymptotic analysis[2], or perturbation techniques

[8,25].

For the mathematical model under the mixed boundary condition in a con_{ﬁned aquifer of semi-inﬁnite thickness, Wilkinson and} Hammond[25]used a perturbation method to give an approximate solution for drawdown changes at the well. Cassiani and Kabala[4]

used the dual integral equation method to develop a Laplace-domain solution that accounts for the effects of wellbore storage, infinitesimal skin, and aquifer anisotropy. Cassiani et al.[5]further used the same method to develop the solutions in Laplace domain suited for constant-head pumping tests and double packer tests treated as the MBPs. Selim and Kirkham[21]used the Gram–Schmidt orthonorma-lization method to develop a steady state solution in a confined aquifer of_{finite horizontal extent. Similar problems under the mixed} boundary conditions also arise in thefield of heat conduction. Among others, Huang [13]used the Weiner–Hopf technique to develop a solution in a semi-infinite slab and Huang and Chang[12]combined ⁎ Corresponding author. 300 Institute of Environmental Engineering, National Chiao

Tung University, 1001 University Road, Hsinchu, Taiwan. E-mail address:[email protected](H.-D. Yeh).

Contents lists available atScienceDirect

Advances in Water Resources

(2)

the Green's function with conformal mapping to develop the solution in an elliptic disk.

In reality, the thickness of aquifer is generally_{ﬁnite. Cassiani et al.}

[5] have developed the Laplace-domain solutions to MBPs for constant-head test based on the infinite aquifer thickness assumption. Their solutions are appropriate for the early time condition when the pressure change caused by the constant-head pumping has not reached the bottom of the aquifer or for the special condition, where the screen length is significantly shorter than the aquifer thickness. Chang and Chen[6] removed such constraints by assuming finite aquifer thickness and treated the well skin effect as a skin factor. They also treated the boundary along the well screen as a Robin (third type) boundary condition and replaced the mixed boundary by homoge-neous Neumann boundary. They considered the wellbore flux entering through the well screen as unknown and discretized the screen length into M segments[7]. To avoid discretizing the well screen, Chang and Yeh [8] developed an analytical solution for a constant-head test performed in a partially penetrating well in an aquifer. They used the dual series equations (DSE) method and perturbation method to solve the MBP along the well which has a well screen extended from the top of the aquifer to any location of the well. According to the image well theory, their solution is applicable to the situation where the middle of the screen of the partially penetrating well is located right at the center of the aquifer. However, their solution can not apply to the case of a partially penetrating well with arbitrary location of the well screen.

This study aims to develop a new model describing a constant-head test performed in a flowing partially penetrating well for arbitrary location of the well screen in an aquifer of afinite thickness in depth. The solution of the model is based on the following assumptions: (1) the aquifer is homogeneous and infinite extent in the radial direction; (2) the well has a_{finite radius; (3) the initial head} is constant and uniform throughout the whole aquifer; and (4) the well loss is not considered in the system. The mixed-type boundary condition at the well is handled via the triple series equations (TSE) method. This solution contains infinite series involving the summa-tions of multiplication of integrals, trigonometric funcsumma-tions, and the modified Bessel functions of second kind, where the single and double integrals are presented in terms of trigonometric functions multiplying the associated Legendre functions. The infinite-series solution is difficult to accurately compute due to the oscillatory nature and slow convergence of the multiplied functions. Therefore, Shanks' transform method[19,20]is used to accelerate the evaluation of the Laplace-domain solution and the numerical inversion scheme, Stehfest algorithm[24], is used to_{find the time domain solution. To} the best of our knowledge, this is thefirst paper using the TSE method to solve the mixed boundary value problems in the area of water resources.

2. Mathematical model 2.1. Mathematical statement

Fig. 1shows a schematic representation of a partially penetrating well in a conﬁned aquifer with a ﬁnite thickness of b. The drawdown at the distance r from the well and the distance z from the bottom of the aquifer at time t is denoted as s(r, z, t). The well screen which extends from arbitrary location d1 to d2 is of length l under a

prescribed constant drawdown hw. The hydraulic parameters of the

aquifer are horizontal hydraulic conductivity Kr [L/T], vertical

hydraulic conductivity Kz [L/T], and speciﬁc storage Ss [1/L]. The

governing equation for the drawdown can be written as

Kr ∂ 2 s ∂r2 + 1 r ∂s ∂r ! + Kz∂ 2 s ∂z2 = Ss∂s ∂t: ð1Þ

The prescribed Dirichlet boundary condition for a constant draw-down along the well screen is:

s rðw; z; tÞ = sw d1≤z≤d2: ð2aÞ

A Neumann boundary condition of zeroﬂux is speciﬁed as: ∂s

∂r

j

r = rw

= 0 0≤z≤d1 and d2≤z≤b: ð2bÞ

In addition, the initial condition and other boundary conditions are: s rð; z; 0Þ = 0 ð3Þ s_ð_{∞; z; t}_{Þ = 0} _ð4Þ and ∂s ∂z= 0; z = 0; z = b: ð5Þ

The dimensionless parameters used hereafter are deﬁned in

Table 1. Eqs. (1)–(5) in dimensionless form are, respectively, ∂2 s* ∂ρ2 + 1 ρ∂s*∂ρ +α 2∂2s* ∂ξ2 = ∂s* ∂τ ð6Þ s*ðρ; ξ; τ = 0Þ = 0 ð7Þ s*ðρ = ∞; ξ; τÞ = 0 ð8Þ s*ðρ = 1; ξ; τÞ = 1; ξ1≤ξ≤ξ2 ð9aÞ ∂s* ∂ρ

j

ρ = 1= 0; 0≤ξ≤ξ1 and ξ2≤ξ≤β ð9bÞ ∂s* ∂ξ

j

ξ = 0;ξ = β= 0: ð10Þ

Note that Eqs. (6)–(10) constitute a MBP.

Fig. 1. Schematic representation of a partially penetrating well in a confined aquifer of finite extent with a finite thickness of b.

(3)

2.2. Laplace-domain solution

The detailed derivation for the solution of Eq. (6) with Eqs. (7)– (10) using Laplace transform,ﬁnite Fourier cosine transform, and TSE method is given inAppendix A. The solution for drawdown in an aquifer involving a partially penetrating well is obtained as:

s*ðρ; ξ; pÞ =1₂B 0ð ; pÞK0 pffiffiffipρ K0 pffiffiffip + ∑∞ n = 1 B nð ; pÞK0ðλnρÞ K0ð Þλn cosðηnξÞ ð11Þ

where K0is the modiﬁed Bessel functions of the second kind with

order zero,ηn= (nπ)/β and the coefﬁcients B(0, p) and B(n, p) are

expressed as B 0_ð _{; p}_{Þ = C 0; p}_ð _{Þ + D 0; p}_ð _Þ = C0+ D0 = B0 ð12Þ and B nð ; pÞ = C n; pð Þ + D n; pð Þ = Cn+ Dn = Bn: ð13Þ

The coefﬁcients C0, Cn, D0, and Dnin Eqs. (12) and (13) are calculated

by the following equations C0= 1 +pffiffiffipH0Ω1ð Þμ1 ₋₁ 4 pπΩ3ð Þ+μ1 2 p 1− μ1 π +∑∞ k = 1 2 kðIkCk−λkHkDkÞΩ2ðμ1; kÞ−D0 ffiffiffi p p _H 0Ω1ð Þμ1 ð14Þ Cn= ∑ ∞ k = 1 1 kðIkCk−λkHkDkÞ Ω2ðμ1; kÞf2ðn; μ1Þ−∫ μ1 0 Ω2ðy; kÞ df2ðn; yÞ dy dy + 1 2 ffiffiffi p p H0ðC0+ D0Þ ∫ μ1 0 Ω1ð Þy df2ðn; yÞ dy dy−Ω1ð Þfμ1 2ðn; μ1Þ + 2 pπ Ω3ð Þfμ1 2ðn; μ1Þ−∫ μ1 0 Ω3ð Þy df2ðn; yÞ dy dy −2 sin nμð 1Þ pnπ ð15Þ D0= 1 +ð pffiffiffipH0Ω1ð Þμ2Þ−1 ∑ ∞ k = 1 2 kð−1Þ k Ω2ðμ2; kÞ IðkDk−λkHkCkÞ−C0pffiffiffipH0Ω1ð Þμ2 ð16Þ and Dn=ð−1Þ n ∑∞ k = 1ð−1Þ k1 kðIkDk−λkHkCkÞ × Ω2ðμ2; kÞ⋅f2ðn; μ2Þ−∫ μ2 0 Ω2ðy; kÞ df2ðn; yÞ dy dy + 1 2 ffiffiffi p p H0ðD0+ C0Þ ∫ μ2 0 Ω1ð Þy df2ðn; yÞ dy dy−Ω1ð Þfμ2 2ðn; μ2Þ ð17Þ with μ1=ξ1π = β ð18Þ μ2=ðπ−ξ2π = βÞ ð19Þ Ω1ð Þ = ∫x x 0f1ðu; xÞudu ð20Þ Ω2ðx; kÞ = ∫ x 0f1ðu; xÞ sin kuð Þdu ð21Þ Ω3ð Þ = ∫x x 0f1ðu; xÞf3ðu; xÞdu ð22Þ λn= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi nπα β 2 + p s ð23Þ H nð ; pÞ = Hn = K1ð Þ = Kλn 0ð Þλn ð24Þ f1ðx; aÞ = ffiffiffi 2 p sin xð = 2Þ πpffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffifficos xð Þ− cosð Þa ð25Þ f2ðn; aÞ = P½ nðcosaÞ + Pn−1ðcosaÞ ð26Þ f3ðx; aÞ = 1

4ðln 1ð − cos a + xð ÞÞ− ln 1− cos a−xð ð ÞÞÞ ð27Þ where Pn(cos(⋅)) is the associated Legendre function ([1], p. 335) and

K1is the modiﬁed Bessel functions of the second kind with order one.

The deﬁnitions of functions or variables used in equations above are also listed inTable 2.

The determination of the values of C0, Cn, D0, and Dnin Eqs. (14)–

(17) from Eqs. (A29) to (A48) is described in detail inAppendix A. The coefﬁcients B0and Bnin the drawdown solution (11) can therefore be

determined based on Eqs. (12) and (13).

Theﬂux entering the well screen and the total well discharge obtained using Eq. (11) are respectively given as:

q* 1ð ; ξ; pÞ = −∂s* ρ; ξ; pð Þ ∂ρ

j

ρ = 1= 1 2B0 ffiffiffi p p K1 pffiffiffip K0 pffiffiffip + ∑∞ n = 1 Bnλn K1ð Þλn K0ð Þλn cos nξ βπ ð28Þ Table 1 Dimensionless expressions. Symbol Illustration s⁎(ρ, ξ, τ) s(r, z, t)/sw, dimensionless drawdown

s̄⁎(ρ, ξ, τ) Dimensionless drawdown in Laplace domain

ˆs* ρ; ξ; τð Þ Dimensionless drawdown in Laplace and Fourier domain α ffiffiffiffiffiffiffiffiffiffiffiffiffiffiKz= Kr

p

, anisotropy ratio

β b/rw, dimensionless aquifer thickness

ηn nπ/β

ρ r/rw, dimensionless radial distance

τ Krt/Ssrw2, dimensionless time

λ l/rw, dimensionless screen length

λn ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

nπα=β ð Þ2_{+ p}

q

ξ z/rw, dimensionless vertical distance

ξ1 d1/rw

ξ2 d2/rw

ω l/b, partial penetration ratio

Table 2 Symbol definitions. Symbol Illustration f1(x, a) ffiffiffi 2 p sin xð = 2Þ πpffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffifficos xð Þ− cosð Þa

f2(n, a) [Pn(cos a) + Pn− 1(cos a)], n≥1

f3(x, a) 14ðln 1− cos a + xð ð ÞÞ− ln 1− cos a−xð ð ÞÞÞ

Hn K1(λn)/K0(λn) In n−λnHn μ1 ξ1π/β μ2 π−ξ2π/β Ω1(x) ∫₀ x f1(u, x)udu Ω2(x, k) _∫ 0 x f1 (u, x) sin(ku)du Ω3(x) ∫₀ x

(4)

and Q pð Þ = 1_λ∫ ξ2 ξ1 q* 1ð ; ξ; pÞdξ = 1₂B0 ffiffiffi p p K1 pffiffiffip K0 pffiffiffip − ∑∞ n = 1 ðληnÞ −1_B nλnHn sin nð μ1Þ + −1ð Þ n sin nð μ2Þ ð29Þ

whereλ=l/rwis the dimensionless length of the screen.

2.3. Simpliﬁed solutions

2.3.1. Partially penetrating well (well screen extends from the top of the aquifer)

When the well screen extends from d1to the top of the aquifer, the

coefﬁcients in Eqs. (14)–(17) can be found by setting ξ2=β. The

drawdown can then be determined from solving Eqs. (A17a) and (A17b) which should be identical to the results obtained using

inﬁnity-order perturbation approach in solving DSE in Chang and Yeh

[8].

2.3.2. Fully penetrating well

When the well fully penetrates the entire thickness of the formation, i.e., _ξ1 is zero and ξ2 equals β, the drawdown and the

well discharge can be obtained using Eqs. (11) and (29), respectively, as[11] s*_ð_{ρ; ξ; p}_{Þ =} 1 p K0 pffiffiffipρ K0 pffiffiffip ð30Þ and Q pð Þ = K1 pffiffiffip ffiffiffi p p _K 0pffiffiffip ð31Þ

Eqs. (30) and (31) are identical to the solutions of drawdown and ﬂow rate in Laplace domain given in Yang and Yeh[28].

Fig. 2. The drawdown distribution at dimensionless timeτ=1, 100, 104

andτ=106

(5)

3. Results and discussion

Numerical calculations for the aquifer drawdown and well_flux are performed in PC using the FORTRAN code developed based on the present solutions. Thefirst step in the development of solutions is to determine the coef_{ficients of Laplace-domain solution in} Eq. (11) from using Eq. (A50). The single and double integrals involved in the elements are then computed using the subroutines DQDAG and DTWODQ in IMSL [10,15], respectively. Once the coefficients are known, the second step is to find the infinite summation in Eq. (11) by Shank's transform method. Then thefinal step is to transform the Laplace-domain solution of Eq. (11) into time-domain using IMSL subroutine LINV for the Stehfest method

[24]with eight weighting factors. The inﬁnite summation in the solution can be found more efﬁciently using Shank's transform which consists of a family of nonlinear sequence-to-sequence transformations[20]. Shanks[20]concluded that these transforma-tions are effective when applied to accelerate the convergence of some slowly convergent sequences and may also converge to some divergent sequences.

The solutions can be veriﬁed by calculating the values at the boundary along the test well in Eq. (11). Fig. 2 shows the dimensionless drawdown forβ=100, ξ1= 30,ξ2= 80 and variousρ

atτ=1, 100, 104_{and 10}6_{. As indicated in the}_{ﬁgure, the dimensionless}

drawdown is constant along the well screen and decreases with the increasing dimensionless radial distance atτ=1. In addition, the dimensionless drawdown increases with dimensionless time along the unscreened part of the well. The dimensionless drawdown has larger value in the screen part and smaller value along the unscreened part. Fig. 3 shows the plots of theﬂux along the well screen for β = 100, ξ1= 30 and ξ2= 80 at τ = 1, 100, 104 and 106. The

dimensionless _{flux is non-uniformly distributed and large at the} screen edge due to the verticalflow induced by the presence of well partial penetration. Fig. 4 exhibits the behavior of dimensionless drawdown versus dimensionless timeτ and illustrates the effect of screen length on the drawdown response, where the dimensionless radial distanceρ is 10, the vertical distance ξ is 50 and α=1 for different length of the screen. This figure indicates that the dimensionless drawdown increases with the length of the screen. To test the influence of anisotropy of the aquifer,Fig. 5is plotted for ρ=5, ξ=50, ξ1= 40,ξ2= 60 and various anisotropy α. As can be

observed, the drawdown increases withα. The spatial dimensionless drawdown contours atτ=100, 103_{and 10}4_{are plotted in}_{Fig. 6}_{. The}

dimensionless drawdown increases with dimensionless time at a ﬁxed radial distance and ﬂow is horizontal when the dimensionless

radial distance is large than 80 and the dimensionless time is 104_.

Fig. 7(a) and (b) shows the spatial dimensionless drawdown contours for various_α2_with

ξ1= 200 andξ2= 250 atτ=105and demonstrates

the inﬂuence of anisotropy on the dimensionless drawdown. The ﬂow is almost horizontal at the bottom of the aquifer when the dimensionless radial distance is large than 400 for_α2_{= 1; however,}

theﬂow is vertical at the bottom of the aquifer for α2_=0.5._{Fig. 8}_(a)

and (b) plots the spatial dimensionless drawdown contours for the same length of 50 but different locations of well screen. InFig. 8(a), the screen is symmetric withξ1= 12.5 andξ2= 37.5 and inFig. 8(b)

the screen extends from the top of the aquifer withξ1= 25 and

ξ2= 50 atτ=105. Since the screen is symmetric about the middle line

of the aquifer, the drawdown contours are symmetric as shown in

Fig. 6(a). Fig. 9 illustrates the spatial dimensionless drawdown contours forβ=200, ξ1= 100 andξ2= 150 atτ=107. The direction

of_{ﬂow is upward when the radial distance is far from the pumping} well and it is downward when the radial distance is close to the well screen.

4. Concluding remarks

This paper developed a new semi-analytical solution for describing the drawdown response for a constant-head test performed in a partially penetrating well in an aquifer of infinite radial extent andfinite vertical extent, where the well screen is installed within any part of the well. The Laplace andfinite cosine Fourier transforms in conjunction with TSE method are used to solve the mixed-type boundary and initial values problem for a partially penetrating well in an aquifer of afinite thickness. The present solutions can be reduced to the solutions given in Yang and Yeh[28]

for a fully penetrating well in an aquifer of afinite thickness. In addition, they are also equal to the results obtained using in finity-order perturbation approach for a partially penetrating well of a well screen extending from the top of the aquifer presented in Chang and Yeh[8]. The flux estimated from the solution is non-uniformly distributed along the screen and with a local peak at the edge, due to the verticalflow induced by the effect of well partial penetration.

These solutions are particularly useful for practical applications since they can be used to evaluate the sensitivities of the input parameters in a mathematical model (e.g.,[14]and[9]), to identify the hydraulic parameters if coupling with the extended Kalmanﬁlter (e.g.,[16]and[33]) or an optimization approach such as the nonlinear least-squares (e.g.,[30]and[31]) or simulated annealing (e.g.,[17]

and[32]) in the analysis of aquifer data, and to validate a numerical solution[36].

Acknowledgements

This study was supported by the Taiwan National Science Council under the grant NSC 96-2221-E-009-087-MY3. The authors would like to thank three anonymous reviewers for their valuable and constructive comments that help improve the clarity of our presentation.

Appendix A

The Laplace-domain solution for dimensionless drawdown can be obtained by taking the Laplace transform with respect to time and the ﬁnite Fourier cosine transform with respect to ξ. The deﬁnition of Laplace transform is:

s*ðρ; ξ; pÞ = Lp½s*ðρ; ξ; τÞ; τ→p = ∫ ∞ 0

s*ðρ; ξ; τÞe−pτdτ ðA1Þ

where s̄⁎(ρ, ξ, p) is the dimensionless drawdown in Laplace domain. Fig. 3. The distribution ofﬂux along the well screen at different dimensionless time for

(6)

Fig. 4. Type curve for drawdown forρ=10, ξ=50, α=1, β=100 and various penetration lengths.

(7)

Taking the Laplace transform of Eq. (6) and Eqs. (8)–(10) with the initial condition in Eq. (7), the problem reads:

∂2 s* ∂ρ2 + 1 ρ∂s*∂ρ +α 2∂2s* ∂ξ2−ps* = 0 ðA2Þ s*ðρ = ∞; ξ; pÞ = 0 ðA3Þ s*ðρ = 1; ξ; pÞ =1_p; ξ1≤ξ≤ξ2 ðA4aÞ ∂s* ∂ρ

j

ρ = 1= 0; 0≤ξ ≤ ξ1 and ξ2≤ξ≤β ðA4bÞ

(8)

Fig. 7. The spatial drawdown contours at dimensionless timeτ=105_for_{β=250 and various α}2_.

Fig. 8. The spatial drawdown contours at dimensionless timeτ=106

(9)

∂s*

∂ξ

j

ξ = 0;ξ = β= 0 ðA5Þ

Theﬁnite cosine Fourier transform with respect to ξ is then deﬁned as follows[23]:

ˆs* ρ; n; pð Þ = Fc½s*ðρ; ξ; pÞ; ξ→n = ∫ β 0

s*ðρ; ξ; pÞ cos ηð nξÞdξ ðA6Þ

where ˆs*_ð_{ρ; n; p}_{Þ is the dimensionless drawdown after ﬁnite cosine} Fourier transform. Substituting Eq. (A6) into Eqs. (A2), (A3) and (A5) results in the Bessel differential equation as

∂2_ˆs* ∂ρ2 + 1 ρ∂ˆs*∂ρ−λ 2 nˆs* = 0 ðA7Þ

with the boundary condition

ˆs* ρ = ∞; n; pð Þ = 0: ðA8Þ

The general solution to Eq. (A7) with the boundary condition (A8) is[3]

ˆs* ρ; n; pð Þ = A n; pð ÞK0ðλnρÞ ðA9Þ

where A(n, p) can be found from using the mixed-type boundary condition (A4a) and (A4b). The inverse of theﬁnite cosine Fourier transform is

s*ðρ; ξ; pÞ =_β1ˆs* ρ; 0; pð Þ +2_β∑∞

n = 1 ˆs* ρ; n; pð Þ cos ηð nξÞ: ðA10Þ

Thus, the solution in_{ξ domain obtained by inserting Eq. (A9) into} Eq. (A10) is s*ðρ; ξ; pÞ =_β1A 0ð ; pÞK0ðpffiffiffipρÞ + 2 β∑ ∞ n = 1 A nð ; pÞK0ðλnρÞ cos ηð nξÞ ðA11Þ

with its derivative with respect toρ given by ∂s* ∂ρðρ; ξ; pÞ=− 1 βA 0ð ; pÞ ffiffiffi p p _K 1ðpffiffiffipρÞ− 2 β∑ ∞ n = 1 A n_ð _{; p}_ÞλnK1ðλnρÞcos ηð nξÞ: ðA12Þ Substituting Eq. (A11) into Eqs. (A4a) and (A12) into Eq. (A4b) results in a system of TSE as

1 βA 0ð ; pÞK0ðpffiffiffipÞ + 2 β∑ ∞ n = 1 A nð ; pÞK0ð Þ cos ηλn ð nξÞ = 1 p; ξ1≤ξ≤ξ2 ðA13aÞ 1 βA 0ð ; pÞ ffiffiffi p p K1 ffiffiffi p p ð Þ + 2_β∑∞ n = 1 A nð ; pÞλnK1ð Þ cos ηλn ð nξÞ = 0; 0≤ξ≤ξ1; ξ2≤ξ≤β: ðA13bÞ Introduce B nð ; pÞ = 2A n; pð ÞK0ð Þ = βλn ðA14Þ and x =ξπ = β: ðA15Þ

Therefore,ηnξ=nx and the TSE of Eqs. (A13a) and (A13b) can be

rearranged as[22]: 1 2B 0ð ; pÞ ffiffiffi p p _{H 0}_ð _{; p}_{Þ + ∑}∞ n = 1 B nð ; pÞλnHncos nxð Þ = 0; 0≤x≤ξ_β1π ðA16aÞ 1 2B 0ð ; pÞ + ∑ ∞ n = 1 B nð ; pÞ cos nxð Þ = 1_p; ξ1 βπbx≤ξβ2π ðA16bÞ 1 2B 0ð ; pÞ ffiffiffi p p _{H 0} ; p ð Þ + ∑∞ n = 1 B n_ð _{; p}_ÞλnHncos nxð Þ = 0; ξ2 βπ≤x≤π: ðA16cÞ The symbol Hnis defined in Eq. (24) and H0is from Hnwhen n = 0.

Our goal now is to determine the coefﬁcients B(0, p) and B(n, p) in Eqs. (A16a)–(A16c). For convenience, the coefﬁcients B(0, p) and B(n, p) are expressed as B0and Bn, respectively, as in Eqs. (12) and (13). To solve

Fig. 9. The spatial drawdown contours as at dimensionless timeτ=107

(10)

the TSE in Eq. (A16a)–(A16c), we further split it into the following two DSE ([22], p. 192) 1 2ðC0+D0ÞpffiffiffipH 0ð ; pÞ+ ∑ ∞ n = 1 Cn+Dn ð ÞλnHncos nxð Þ = 0; 0≤ x ≤ μ1 ðA17aÞ 1 2C0+∑ ∞ n = 1 Cncos nxð Þ = 1 p; μ1bx≤π ðA17bÞ 1 2D0+∑ ∞ n = 1 Dncos nxð Þ = 0; 0b x ≤ π − μ2 ðA18aÞ 1 2ðC0+D0Þ ffiffiffi p p H 0ð ; pÞ+ ∑∞ n = 1 Cn+Dn ð ÞλnHncos nxð Þ = 0; π−μ2≤x≤π ðA18bÞ whereμ1andμ2are defined by Eqs. (18) and (19), respectively. With

Eqs. (12) and (13), Eqs. (A17a) and (A18b) are equal to Eqs. (A16a) and (A16c), respectively, and the sum of Eqs. (A17b) and (A18a) in the range ofμ1bx≤π−μ2 is equal to Eq. (A16b). Eqs. (A17a)–

(A17b) and (A18a)–(A18b) are regarded as dual series relations and by means of them, the coefﬁcients C0, D0, Cn and Dn can be

determined.

Eqs. (A17a) and (A17b) can be solved by the following procedure given in Sneddon ([22], p. 161). Assume that for 0≤x≤μ1

1 2C0+ ∑ ∞ n = 1 Cncos nxð Þ = cos x 2 ∫μ1 x h1ð Þdyy ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cosx− cos p y ðA19Þ

where h1(y) is an unknown function to be determined. Using

Eqs. (A17b) and (A19), for the full range 0≤x≤π, the coefﬁcients C0

and Cncan then be expressed as ([22], (5.4.56), (5.4.57))

C0= 2 π p ∫πffiffiffi2 μ1 0 h1ð Þdy + ∫y π μ1 1 pdy ðA20Þ and Cn= 2 π 2p ∫πffiffiffi2 μ1

0h1ð Þ Py½ nðcosyÞ + Pn−1ðcosyÞdy + ∫πμ1

1

pcos nyð Þdy

: ðA21Þ

The function h1(y) can be determined using Eq. (A17a) for

0≤x≤μ1. Integrating Eq. (A17a), one can obtain

1 2C0pffiffiffipH 0ð ; pÞx + ∑ ∞ n = 1 Cnsin nxð Þ =∫x 0 ∑ ∞ n = 1 CnIncos nuð Þ− 1 2D0 ffiffiffi p p H 0ð ; pÞ− ∑∞ n = 1 DnλnHncos nuð Þ du =∫x 0F uð Þdu ðA22Þ

where F(u) is defined as the integrand on the RHS of Eq. (A22). Substituting Eq. (A21) into Eq. (A22), one canfind that h1(y) satisfies

the following equation: ([22], p. 161, Eq. (5.4.58))

∫μ1 0 h1ð Þy p ∑1ffiffiffi₂ ∞ n = 1 PnðcosyÞ + Pn−1ðcosyÞ ½ sin nxdy =∫x 0F uð Þdu− 1 2 ffiffiffi p p _{H 0}_ð _{; p}_ÞC 0x− ∑ ∞ n = 1 2 π∫ π μ1 1

pcos nuð Þdu sin nxð Þ: ðA23Þ

The summation term on the left-hand side of Eq. (A23) can be expressed as ([22], p. 59, Eq. (2.6.31)) 1ffiffiffi 2 p ∑∞ n = 1 PnðcosyÞ + Pn−1ðcosyÞ ½ sin nx = cos x 2 Heavðx−yÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cosy− cos p _x ðA24Þ

where Heav(X) is the Heaviside unit step function deﬁned as

Heavð Þ =X 0 X_b0 1= 2 X = 0 1 XN 0: 8 < : ðA25Þ

Substituting Eq. (A24) into Eq. (A23) yields ∫μ1 0 h1_{ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi}ð ÞH x−yy ð Þ cosy− cos p _xdy = secx 2 ∫ x 0F uð Þdu− 1 2 ffiffiffi p p H 0ð ; pÞC0x− ∑ ∞ n = 1 2 π∫ π μ1 1

pcos nuð Þdu sin nxð Þ

( )

: ðA26Þ With Eq. (A25), Eq. (A26) can be expressed alternatively as ∫ x 0 h1ð Þy ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cosy_{− cos} p _xdy = secx 2 ×∫x 0F uð Þdu− 1 2 ffiffiffi p p _{H 0}_ð _{; p}_ÞC 0x−∑ ∞ n = 1 2 π∫ π μ1 1

( )

0≤ x b μ1:

ðA27Þ Then, the function h1(y) found based on Sneddon ([22], p. 162, Eq.

(5.4.60)) is h1ð Þ =y 2 π d dy∫ y 0 sin xð = 2Þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cosx− cos p _y × ∫x 0F uð Þdu− 1 2 ffiffiffi p p _{H 0}_ð _{; p}_ÞC 0x− ∑ ∞ n = 1 2 π∫ π μ1 1

( )

dx: ðA28Þ By integrating Eq. (A28) and substituting it into Eqs. (A20) and (A21), the coefﬁcients C0and Cncan then be expressed as Eqs. (14)

and (15), respectively.

For computational convenience, Eqs. (14) and (15) can be written as a vector equation

C = I−Xð Þ1YD + I−Xð Þ1Z ðA29Þ

where I is an (n + 1) × (n + 1) identity matrix; X = [xi,j] and Y = [yi,j]

are (n + 1)× (n + 1) matrices; CT_{= [C}

0, C1,…, Cn], DT= [D0, D1,…, Dn],

and ZT_{= [z}

1, z2,…, zn + 1] are column vectors. The elements in the

matrices and vectors are deﬁned as

x1;1= 0 ðA30Þ x1;j= 2 j−1 ð ÞIj−1Ω2ðμ1; j−1Þ 1 +pffiffiffipH0Ω1ð Þμ1 ðA31Þ xij= 1 j−1 ð ÞIj−1 Ω2ðμ1; j−1Þf2ði−1; μ1Þ−∫ μ1 0 Ω2ðy; j−1Þ df2ði−1; yÞ dy dy ðA32Þ xi1= 1 2 ffiffiffi p p _H 0 ∫ μ2 0 Ω1ð Þy df2ði−1; yÞ dy dy−Ω1ð Þfμ1 2ði−1; μ1Þ ðA33Þ y11= − ffiffiffi p p _H 0Ω1ð Þμ1 1 +pffiffiffipH0Ω1ð Þμ1 ðA34Þ

(11)

y1j= −2 j−1 ð Þλj−1Hj−1Ω2ðμ1; j−1Þ 1 +pffiffiffipH0Ω1ð Þμ1 ðA35Þ yi1= 1 2 ffiffiffi p p _H 0 ∫ μ1 0 Ω1ð Þy df2ði−1; yÞ dy dy−Ω1ð Þfμ1 2ði−1; μ1Þ ðA36Þ yij= 1 j₋₁ ð Þλj−1Hj−1 ∫ μ1 0 Ω2ðy; j−1Þ df2ði−1; yÞ dy dy−Ω2ðμ1; j−1Þf2ði−1; μ1Þ ðA37Þ z1= 4 pπΩ3ð Þ +μ1 2p 1−μπ1 1 +pffiffiffipH0Ω1ð Þμ1 ðA38Þ zi= 2 pπ Ω3ð Þfμ1 2ði−1; μ1Þ−∫ μ1 0 Ω3ð Þy df2ði−1; yÞ dy dy −2 sin iðð−1Þμ1Þ i−1 ð Þπp ðA39Þ where i and j goes from 2 to n and the functions f1(⋅), f2(⋅) and f3(⋅) are

deﬁned in Eqs. (25)–(27), respectively.

Similarly, Eqs. (A18a)–(A18b) can be solved by setting x′=π−x and Dn′=(−1)nDn. Eqs. (A18a)–(A18b) is rewritten as

1 2D ′ 0+ ∑ ∞ n = 1 D′ncos nxð ′Þ = 0; μ2bx′≤π ðA40aÞ 1 2 D′0+ C0 ffiffiffip_p_H 0+ ∑ ∞ n = 1 D′ n+ð−1Þ n Cn λnHncos nxð ′Þ = 0; 0≤x′≤μ2 ðA40bÞ and the vector equation for solving coefficients D0and Dnis expressed

as

D = I−˜X −1˜YC ðA41Þ

where X̃= [x̃ij] andỸ=[ỹij] are (n + 1) × (n + 1) matrices with the

elements ˜x1;1= 0 ðA42Þ ˜x1;j= 2ð−1Þj−1 j−1 ð Þ Ij−1Ω2ðμ2; j−1Þ 1 +pffiffiffipH0Ω1ð Þμ2 ðA43Þ ˜xi1= 1 2 ffiffiffi p p _H 0 ∫ μ2 0 Ω1ð Þy df2ði−1; yÞ dy dy−Ω1ð Þfμ1 2ði−1; μ2Þ ðA44Þ ˜xij= ffiffiffi 2 p π −1 ð Þi−1_ð₋₁_Þj−1 j−1 ð Þ Ij−1 Ω2ðμ2; j−1Þf2ði−1; μ2Þ−∫μ02Ω2ðy; j−1Þ df2ði−1; yÞ dy dy ðA45Þ ˜y11= − ffiffiffi p p _H 0Ω1ð Þμ2 1 +pffiffiffipH0Ω1ð Þμ2 ðA46Þ ˜y1j= −2 −1ð Þj−1 j−1 ð Þ λj−1Hj−1Ω2ðμ2; j−1Þ 1 +pffiffiffipH0Ω1ð Þμ2 ðA47Þ ˜yi1=ð−1Þ i−1 2 ffiffiffi p p H0 ∫ μ2 0 Ω1ð Þy df2ði−1; yÞ dy dy−Ω1ð Þfμ2 2ði−1; μ2Þ ðA48Þ ˜yij= ð−1Þ j−1_ð₋₁_Þi−1 j−1 ð Þ λj−1Hj−1 × ∫μ2 0 Ω2ðy; j−1Þ df2ði−1; yÞ dy dy−Ω2ðμ2; j−1Þf2ði−1; μ2Þ : ðA49Þ

If n tends to inﬁnity, Eq. (11) would give the exact solution for the drawdown. However, it would give a relatively accurate result even a

ﬁnite number of n is considered. Let n vary from 1 to N, where N is an arbitraryﬁnite number. Substituting Eq. (A41) into Eq. (A29), the elements in C column vector can be expressed as

c_j−1= N + 1∑ i = 1 φi;j

zi; j = 1; 2; 3; …; N + 1 ðA50Þ

withφi,jrepresents (i, j)th element in the matrix [I− (I − X)− 1Y(I−

X̃)− 1Ỹ]− 1_(I_{− X)}− 1_.

Once the coefﬁcients C0and Cnare known, the coefﬁcients D0and

Dncan then be obtained from Eq. (A41).

References

[1] Abramowitz M, Stegun IA. Handbook of Mathematical Functions. New York: Dover Publications; 1970.

[2] Bassani JL, Nansteel MW, November M. Adiabatic-isothermal mixed boundary conditions in heat transfer. J Heat Mass Transfer 1987;30:903–9.

[3] Carslaw HS, Jaeger JC. Conduction of heat in solids. 2nd Ed. Oxford: Clarendon; 1959.

[4] Cassiani G, Kabala ZJ. Hydraulics of a partially penetrating well: solution to a mixed-type boundary value problem via dual integral equations. J Hydrol 1998;211:100–11. [5] Cassiani G, Kabala ZJ, Medina Jr MA. Flowing partially penetrating well: solution to

a mixed-type boundary value problem. Adv Water Resour 1999;23:59–68. [6] Chang CC, Chen CS. An integral transform approach for a mixed boundary problem

involving aﬂowing partially penetrating well with inﬁnitesimal well skin. Water Resour Res 2002;38:1071.

[7] Chang CC, Chen CS. Aﬂowing partially penetrating well in a ﬁnite-thickness aquifer: a mixed-type initial boundary value problem. J Hydrol 2003;271:101–18. [8] Chang YC, Yeh HD. New solutions to the constant-head test performed at a

partially penetrating well. J Hydrol 2009,doi:10.1016/j.jhydrol.2009.02.016. [9] Chen YJ, Yeh HD. Parameter estimation/sensitivity analysis for an aquifer test with skin

effect. Ground Water2009;47:287–99,doi:10.1111/j.1745-6584.2008.00530.x. [10] GeraldCF,WheatleyPO.Appliednumericalanalysis.4thed.California:Addison-Wesley;

1989.

[11] Hantush MS. Hydraulics of wells. In: Chow VT, editor. Advances in hydroscience, Vol. 1. New York: Academic Press; 1964.

[12] Huang SC, Chang YP. Anisotropic heat conduction with mixed boundary conditions. J Heat Transfer 1984;106:646–8.

[13] Huang SC. Unsteady-state heat conduction in semi-inﬁnite regions with mixed-type boundary conditions. J Heat Transfer 1985;107:489–91.

[14] Huang YC, Yeh HD. The use of sensitivity analysis in on-line aquifer parameter estimation. J Hydrol 2007;335:406–18,doi:10.1016/j.jhydrol.2006.12.007. [15] IMSL. Math/library, special functions. Houston: Visual Numerics Inc; 1997. [16] Leng CH, Yeh HD. Aquifer parameter identiﬁcation using the extended Kalman

ﬁlter. Water Resour Res 2003;39:1062,doi:10.1029/2001WR000840.

[17] Lin YC, Yeh HD. Trihalomethane species forecast using optimization method: genetic algorithm and simulated annealing. Journal of Computing in Civil Engineering, ASCE 2005;19:248–57,doi:10.1061/(ASCE)0887-3801(2005)19:3(248).

[18] Noble B. Methods based on the Wiener–Hopf techniques. New York: Pergamon Press; 1958.

[19] Peng HY, Yeh HD, Yang SY. Improved numerical evaluation for the radial groundwaterﬂow equation. Adv Water Resour 2002;25:663–75.

[20] Shanks D. Non-linear transformations of divergent and slowly convergent sequence. J Math Phys 1955;34:1–42.

[21] Selim MS, Kirkham D. Screen theory for wells and soil drainpipes. Water Resour Res 1974;10:1019–30.

[22] Sneddon IN. Mixed boundary value problems in potential theory. Amsterdam: North-Holland; 1966.

[23] Sneddon IN. The use of integral transforms. New York: McGraw-Hill; 1972. 540 pp. [24] Stehfest, H., Algorithm 368, Numerical inversion of Laplace transforms, Comm.

ACM, 13 (1 and 10) (1970), p. 47–49 and p.624.

[25] Wilkinson D, Hammond PS. A perturbation method for mixed boundary-value problems in pressure transient testing. Trans Porous Media 1990;5:609–36. [26] Yang SY, Yeh HD. Solution forﬂow rates across the wellbore in a two-zone

conﬁned aquifer. J Hydraul Eng ASCE 2002;128:175–83.

[27] Yang SY, Yeh HD. Laplace-domain solutions for radial two-zoneﬂow equations under the conditions of constant-head and partially penetrating well. J Hydraul Eng ASCE 2005;131:209–16.

[28] YangSY,YehHD.Anovelanalyticalsolutionforconstant-headtestinapatchyaquifer.IntJ NumerAnalMethodsGeomech2006;30:1213–30,doi:10.1002/nag.523.

[29] Yedder RB, Bilgen E. On adiabatic-isothermal mixed boundary conditions in heat transfer. Wärme Stoffübertragung 1994;29:457–60.

[30] Yeh HD. Theis' solution by nonlinear least-squares andﬁnite-difference Newton's method. Ground Water 1987;25:710–5.

[31] Yeh HD, Han HY. Numerical identiﬁcation of parameters in leaky aquifers. Ground Water 1989;27:655–63.

[32] Yeh HD, Chang TH, Lin YC. Groundwater contaminant source identiﬁcation by a hybrid heuristic approach. Water Resour Res 2007;43:W09420, doi: 10.1029/2005WR004731.

(12)

[33] Yeh HD, Huang YC. Parameter estimation for leaky aquifers using the extended Kalman ﬁlter and considering model and data measurement uncertainties. J Hydrol 2005;302:28–45,doi:10.1016/j.jhydrol.2004.06.035.

[34] Yeh HD, Yang SY, Peng HY. A new closed-form solution for a radial two-layer drawdown equation for groundwater under constant-ﬂux pumping in a ﬁnite-radius well. Adv Water Res 2003;26:747–57.

[35] Yeh HD, Yang SY. A novel analytical solution for a slug test conducted in a well with aﬁnite-thickness skin. Adv Water Resour 2006;29:1479–89,doi:10.1016/j. advwatres.2005.11.002.

[36] Zheng C, Bennett GD. Applied contaminant transport modeling. Second Edition. New York: Wiley; 2002. p. 621.