On the solutions of modeling a slug test performed in a two-zone confined aquifer

全文

(1)On the solutions of modeling a slug test performed in a two-zone confined aquifer Yang Shaw-Yang & Yeh Hund-Der. Abstract The slug test is a popular field method for estimating aquifer parameters. As a finite-thickness skin is present at a wellbore, the result of test-data analysis may represent average hydraulic properties of the skin and formation zones. Moench and Hsieh (Int. Assoc. of Hydrogeology, The IAH 17th International Congress on Hydrology of Rocks of Low Permeability, Int. Assoc. of Hydrogeology. Tucson, Arizona, 7–12 January) provided the Laplace-domain solution for a confined aquifer with a finite-thickness skin that satisfies the governing equations, boundary conditions, and continuity requirements at the interface of the skin and formation zones. Two new formulae presented here relate the modified zero-order Bessel functions of the first and second kinds to their first and second derivatives. Their solution is evaluated using the modified Crump approach for the Laplace inversion. The numerical values obtained by the modified Crump approach are compared with others based on finiteelement modeling. The results indicate that the values based on finite-element modeling exhibit poor accuracy, attributed to errors in the mass balance equation and interface condition in their mathematical model. Résumé Le slug-test représente une méthode populaire de terrain pour estimer les paramètres d’un aquifère. Du fait de la présence d’une couche à l’épaisseur limitée sur un site de forages, le résultat de l’analyse des données du test pourrait représenté le terrain mais aussi cette couche. Moench AF et Hsieh PA (Association Internationale d’Hydrogéologie, 7–12 Jan 1985) Analyses de données de slug test dans un puits présentant une couche d’épaisseur finie. 17ème congrès international de l’AIH Received: 3 July 2006 / Accepted: 8 August 2006 Published online: 18 October 2006 © Springer-Verlag 2006 Y. Shaw-Yang Department of Civil Engineering, Vanung University, Chungli, Taiwan Y. Hund-Der ()) Institute of Environmental Engineering, National Chiao-Tung University, Hsinchu, Taiwan e-mail: hdyeh@mail.nctu.edu.tw Fax: +886-3-5726050 Hydrogeology Journal (2007) 15: 297–305. sur l’Hydrologie des roches de faible perméabilité, Association Internationale d’Hydrogéologie. Tucson, Arizona, 7–12 Jan. apportent une solution du domaine de Laplace pour un aquifère confiné avec une épaisseur limitée qui satisfait l’équation en question, les conditions aux limites, et les besoins de continuité à l’interface de la couche et de l’encaissant. Deux nouvelles formules présentées dans cet article met en relation les fonctions de Bessel modifiées d’ordre zéro du premier et second type, avec leur dérivées première et deuxième. Leurs solutions sont évaluées en utilisant une approche de Crump modifiée pour l’inversion de Laplace. Les valeurs numériques obtenues avec l’approche de Crump sont comparées avec d’autres approches utilisant les modèles aux éléments finis. Les résultats indiquent que les valeurs basées sur les modèles aux éléments finis montrent une faible précision, du fait d’erreurs dans le bilan hydrologique et les conditions aux limites dans le modèle. Resumen Los ensayos slug o tipo pulso son métodos populares de campo para estimar parámetros de acuíferos. Debido a que en los sondeos existe una membrana de espesor finito, el resultado del análisis de los datos del ensayo pueden representar las propiedades hidráulicas promedio de las zonas de la formación y dicha membrana. Moench AF y Hsieh PA (Asociación Internacional de Hidrogeólogos, 7–12 Enero 1985) Análisis de datos de prueba slug en un pozo con membrana de espesor finito. 17 Congreso Internacional de la Asociación Internacional de Hidrogeólogos sobre Hidrología de Rocas de Baja Permeabilidad. Asociación Internacional de Hidrogeólogos. Tucson, Arizona. 7–12 Enero. han aportado la solución en el ámbito de Laplace para un acuífero confinado con una membrana de espesor finito que satisface los requerimientos de las ecuaciones dominantes, las condiciones restrictivas, y de continuidad en la interfase de las zonas de formación y membrana. Las dos nuevas fórmulas que se presentan aquí relacionan las funciones modificadas Bessel de orden cero de la primera y segunda clase con su primera y segunda derivadas. Se evalúa su solución usando el enfoque modificado de Crump para la inversión Laplace. Los valores numéricos obtenidos mediante el enfoque modificado de Crump se comparan con otros valores numéricos basados en modelos de elementos finitos. Los resultados indican que los valores basados en el modelo de elementos finitos DOI 10.1007/s10040-006-0100-x.

(2) 298. presentan poca precisión, que se atribuye a errores en la ecuación de balance de masa y a las condiciones de interfase en su modelo matemático. Keywords Slug test . Analytical solutions . Skin effect . Groundwater hydraulics . Shanks method. Introduction The slug test is a popular field method to evaluate in-situ aquifer parameters in low permeability aquifers. The test involves a sudden withdrawal from or injection into a test well, resulting in an instantaneous raising or lowering of the water level in the well. An instantaneous head change is thus imposed within a wellbore and the recovery (or falloff) of the water level in the well is measured, often by means of a pressure transducer. The aquifer parameters, e.g., transmissivity and storativity, can then be obtained from the analysis of test data. For engineering applications, the slug test is a simple, quick, and economical approach. An aquifer is considered a radial two-zone or composite system if the formation properties near the wellbore are apparently changed due to drilling or development of the well. Drilling may create a wellbore skin between the wellbore and formation. The invasion of drilling mud into the aquifer may produce a positive wellbore skin that has lower hydraulic conductivity than that of the undisturbed formation. On the other hand, the extensive well development and/or substantial fracturing of the borehole wall may increase the hydraulic conductivity as a negative wellbore skin. The objectives of this work are to give a numerical approach for evaluating the Moench and Hsieh (1985b) Laplace-domain solution and compare this with the numerical results obtained from the Yang and Gates (1997) finite-element model. The numerical approach includes the modified Crump approach (IMSL 1997) for performing the numerical inversion for the Moench and Hsieh (1985b) solution (hereinafter referred to as MHS), the series expansion for approximating the Bessel functions in MHS, and the Shanks method (Shanks 1955; Yeh et al. 2003; Yang and Yeh 2005) for accelerating the convergence when evaluating the series of the Bessel functions. It is also proven that MHS satisfies the mathematical model of the slug test in a well with a finite-thickness skin. In addition, two new formulae related the modified zero-order Bessel functions of the first and second kinds to their first and second derivatives, and these formulae may be of use in future studies.. does. Under this condition, the skin effect can lead to very low estimates of hydraulic conductivity if using the typecurve fitting method of Cooper et al. (1967). Moench and Hsieh (1985a) commented on the evaluation of slug tests in a finite-thickness skin by Faust and Mercer (1984). They stated that when the specific storage of the skin is negligibly small, the finite-thickness skin solution becomes equivalent to the infinitesimally thin skin solution published by Ramey et al. (1975). Under a finite-thickness skin condition, the skin properties in the test interval control the early time response, whereas the formation properties relate to the late time response. This causes the well response to depart significantly from the conventional type curves. Moench and Hsieh (1985b) provided a Laplace-domain solution for the mathematical model representing a drill-stem test in a well surrounded by a finite-thickness skin. They found that the standard methods of test-data analysis are adequate for open-well slug tests and may differ markedly in the pressure response for pressurized slug tests. Further, Sageev (1986) investigated the effects of well storage and wellbore skin in a confined aquifer system. He obtained a similar result to that of Moench and Hsieh (1985b). Karasaki et al. (1988) developed the various models of slug-tests and related solutions for linear flow, radial flow with boundaries, two zones, and concentric composite aquifer systems with different flow geometricies between the inner and outer regions; they provided type curves for each solution and stated that slug tests suffer the problems of non-uniqueness in matching the test data to type curves. Karasaki (1990) used the time convolution method of Duhamel’s theorem to derive the analytical solution of Moench and Hsieh (1985b). The systematized procedure and analysis method were proposed for a drill-stem test. Recently, Yang and Gates (1997) constructed a finiteelement model for a slug test in a confined aquifer system containing a finite-thickness skin. Their analysis shows that the effect of a low-permeability skin on the estimates of hydraulic conductivity could be minimized by the use of late-time test data. Butler and Healey (1998) investigated the estimate of hydraulic conductivities obtained through pumping or slug tests; they simply noted that the hydraulic conductivity estimate from pumping-test data is, on average, larger than that from a series of slug-test data in the same formation. Zlotnik and McGuire (1998a,b) and Hemker (1999) developed the analytical solutions of the transient well inflow in a multi-layered aquifer system. McElwee (2002) examined and used four techniques, including an extended capability nonlinear model, sensitivity analysis, correction for acceleration and velocity effects, and multiple slug tests to improve the analyses of slug tests.. Literature review Faust and Mercer (1984) provided an infinite-aquifer solution on the response of slug tests in wells containing a finite-thickness skin. They assumed that the skin has a much lower permeability than the adjacent formation Hydrogeology Journal (2007) 15: 297–305. Mathematical model The governing equations describing the hydraulic head distributions within the skin and formation DOI 10.1007/s10040-006-0100-x.

(3) 299. zones can respectively be written as (Moench and Hsieh 1985a,b). and 0. @ h1 1 @h1 S1 @h1 ¼ ; rw < r < rs þ r @r @r2 T1 @t 2. ð1Þ. and @ 2 h2 1 @h2 @h2 ¼ ; r > rs þ r @r @r2 @t. ð2Þ. where the subscripts 1 and 2 respectively denote the skin and formation zones, T is the transmissivity, and S is the storativity. The symbol h is the hydraulic head within the skin or formation zone, r is the radial distance from center of well, rs is the distance from center of well to outer boundary of skin zone, t is the test time, and rw is the well radius. The hydraulic heads are initially assumed to be zero within both the skin and formation zones, that is h1 ¼ h2 ¼ 0; t ¼ 0. ð3Þ. hw ¼ h0 ; t ¼ 0. ð4Þ. where h0 is the initial head. The conservation of mass at a well requires that @h1 2 dhw rc ¼ 2K1 brw ; r ¼ rw ð5Þ dt @r r¼rw and the outer boundary condition at infinitity is h2 ¼ 0; r ! 1. r ¼1. 0. 0. h2 ¼ 0; r ! 1 0. 0. 0. ð12Þ. 0. h1 ¼ h2 ; r ¼ rs. ð13Þ. and 0. dh1 dh 0 0 ¼ 20 ; r ¼ rs dr0 dr. ð14Þ. The Laplace-domain solutions for dimensionless hydraulic heads within both the skin and formation zones given by Moench and Hsieh (1985b) are respectively 0 0 Δ1 K0 qr Δ2 I0 qr 0 h1 ¼ ð15Þ c1 Δ1 c2 Δ2 and. ð6Þ. 0 0 0 Δ1 K0 qrs Δ2 I0 qrs K0 qr h2 ¼ ½c1 Δ1 c2 Δ2 K0 qrs0 0. where rc is the radius of the well casing in the region of changing water level, K is the hydraulic conductivity, and b is the aquifer thickness. The continuity conditions required at the interface between the skin and formation zones are h1 ¼ h2 ; r ¼ rs. ð10Þ. where h0 is the transform dimensionless hydraulic head, p is the Laplace transform variable of dimensionless time, q ¼ p1=2 , ¼ ðÞ1=2 , α=K 1 /K 2 , δ=S s 1 /S s 2 , 0 ¼ rc2 =ð2rw2 Ss2 bÞ, r'= r/rw, and rs ¼ rs =rw . The transform boundary and continuity conditions of Eqs. (5–8) are respectively, ! 0 dh1 0 ðphw 1Þ ¼ ð11Þ dr0 0. 0. The initial condition for a head in a well, hw, is. 0. d 2 h2 1 dh2 0 0 0 þ 0 ¼ q2 h2 ; r > rs r dr0 dr0 2. ð16Þ. ð7Þ. where the variables Δ1, Δ2, c1, and c2 are respectively defined as 0 0 0 0 Δ1 ¼ I0 qrs K1 qrs þ I1 qrs K0 qrs ð17Þ. ð8Þ. 0 0 0 0 Δ2 ¼ K0 qrs K1 qrs K1 qrs K0 qrs. ð18Þ. c1 ¼ pK0 ðqÞ þ qK1 ðq Þ. ð19Þ. and K1. @h1 @h2 ¼ K2 ; r ¼ rs @r @r. Applying the dimensionless transforms and the Laplace transforms to Eqs. (1) and (2) and using the initial condition, Eq. (3), gives the following subsidiary equations:. and 2. 0. 0. d h1 1 dh1 0 0 2 0 02 þ 0 0 ¼ ðqÞ h1 ; 1 < r < rs dr r dr Hydrogeology Journal (2007) 15: 297–305. ð9Þ. c2 ¼ pI0 ðqÞ qI1 ðqÞ. ð20Þ DOI 10.1007/s10040-006-0100-x.

(4) 300. When r′ = 1, the dimensionless head in a well obtained from Eq. (15) is h0w ¼. ½Δ1 K0 ðqÞ Δ2 I0 ðq Þ c1 Δ1 c2 Δ2. ð21Þ. As shown in this section, the MHS satisfies the governing equations, boundary conditions, and continuity requirements at the interface of the skin and formation zones. Taking the derivatives of Eq. (15) with respect to r and letting r′=1 yields ð22Þ. Substituting Eq. (21) into the left-hand side (LHS) of Eq. (11) and using simple algebraic manipulation, one can obtain Eq. (22). Thus, the boundary condition at a wellbore, Eq. (11), is satisfied. 0 With K0(∞)=0 as r′→∞, Eq. (16) yields h2 ¼ 0. Therefore, the solution satisfies the outer boundary condition, Eq. (12). 0 0 Letting r ¼ rs , Eqs. (15) and (16) become 0. h1 ¼. 0. 0. ½Δ1 K0 ðqrs ÞΔ2 I0 ðqrs Þ ¼ c1 Δ1 c2 Δ2. 0. h2. ð23Þ. This shows the continuity of dimensionless hydraulic head at the interface of the skin and formation zones. Taking the derivatives of Eqs. (15) and (16), respec0 0 tively, with respect to r′ and letting r ¼ rs gives. 0 0 2 K1 qrs dh1. ð24Þ. 0 0 ¼ dr0 r ¼rs rs0 ½c1 Δ1 c2 Δ2 and. 0 K1 qrs dh2. 0 0 ¼ dr0 r ¼rs rs0 ½c1 Δ1 c2 Δ2 . ð27Þ. and 0 0 0 2 ¼ I I I qr qr 1 0 2 qr qr0. Moench and Hsieh Solution (1985b). 0. dh1. q ½Δ1 K1 ðqÞ þ Δ2 I1 ðq Þ 0 ¼ c1 Δ1 c2 Δ2 dr0 r ¼1. 0 0 0 dI1 qr 1 ¼ I qr þ I 1 0 qr qr0 dr0. ð28Þ. Based on Eqs. (26–28), one can obtain d 2 h 0 i 1 d h 0 i þ 0 0 I0 qr I0 qr r dr dr0 2 2 2 0 ¼ q I0 qr. ð29Þ. Likewise, d 2 h 0 i 1 d h 0 i þ 0 0 K0 qr K0 qr r dr dr0 2 2 2 0 ¼ q K0 qr. ð30Þ. After taking the first and second derivatives with respect to r′ for Eq. (15) and using Eqs. (29) and (30), one can get 0. 0. d 2 h1 1 dh1 2 0 02 þ 0 0 ¼ ðq Þ h1 r dr dr Similarly, one can obtain 0 0 d 2 h2 1 dh2 2 0 þ 0 ¼ q h2 r dr0 dr0 2. ð31Þ. ð32Þ. According to Eqs. (31) and (32), the solutions of Eqs. (15) and (16) respectively satisfy Eqs. (9) and (10).. Numerical evaluation. 0. ð25Þ. Furthermore, Eq. (24) equals the product of Eq. (25) multiplied by α on both sides. Accordingly, this proves the continuity of flow rate between the skin and formation zones, i.e., Eq. (14). Following steps similar to those demonstrated above by taking the first and second derivatives with respect to r′ for Eqs. (15) and (16) will prove that MHS satisfies the governing equations, Eqs. (9) and (10). McLachlan (1955, pp 192–197) provided the formulae 0 0 dI0 qr ¼ ð q ÞI 1 qr dr0 Hydrogeology Journal (2007) 15: 297–305. ð26Þ. Bessel functions. Equations (15), (16), and (21) include the Bessel functions I0(u), I1(u), K0(u), and K1(u). These functions are approximated here by the formulae given in Watson (1958) and Abramowitz and Stegun (1964). The argument u in these formulae may be divided into two ranges; (0, 10) and (10, ∞) for I0(u) and I1(u), and (0, 2) and (2, ∞) for K0(u) and K1(u) in order to achieve better accuracy. In addition,, the formulae are essentially composed of infinite series and may converge slowly, especially when u is small. Therefore, the Shanks method (Shanks 1955; Yeh et al. 2003; Yang and Yeh 2005) is employed to accelerate the convergence when evaluating the series of the Bessel functions. Each function in Eqs. (15), (16), and (21) is calculated to ten decimal places, and thus they bear the same degree of accuracy as those listed in Abramowitz and Stegun (1964). DOI 10.1007/s10040-006-0100-x.

(5) 301. The Bessel functions of I0(u), I1(u), K0(u), and K1(u) may be evaluated by the formulae given in Abramowitz and Stegun (1964). For 0 ≤ u ≤10, the asymptotic expansions for I0(u) and I1(u) may be expressed respectively as (Abramowitz and Stegun 1964, p 375) 1 2 2 1 2 3 1 2 u u 4u I0 ðuÞ ¼ 1 þ þ 4 2 þ 4 2 þ ð33Þ 2 ð1!Þ ð2!Þ ð3!Þ and I1 ðuÞ ¼. " u 2. 1þ. 1 4. 2. 1. 2 2. 1. 2 3. u u u þ 4 þ 4 þ ð1!Þð2!Þ ð2!Þð3!Þ ð3!Þð4!Þ. where = ¼ 0:577215664901533 is the Euler’s constant. For 2 < u < ∞, K0(u) and K1(u) are respectively taken as (Abramowitz and Stegun 1964, p 378) eu K0 ðuÞ ¼ pffiffiffiffiffiffiffiffi 2u ( ) 12 12 32 12 32 52 1 þ þ ð1!Þð8uÞ ð2!Þð8uÞ2 ð3!Þð8uÞ3 ð39Þ. #. ð34Þ The Bessel functions of I0(u) and I1(u) for 10 < u < ∞ are respectively approximated as (Abramowitz and Stegun 1964, p 377). and eu K1 ðuÞ ¼ pffiffiffiffiffiffiffiffi 2u ( ) 4 12 ð4 12 Þ ð4 32 Þ ð4 12 Þ ð4 32 Þ ð4 52 Þ 1þ þ þ þ ð1!Þð8uÞ ð2!Þð8uÞ2 ð3!Þð8uÞ3. ð40Þ. Shanks method. I 0 ðu Þ. ) The Shanks method, which is also called the ε-algorithm (Kahaner 1972), is a non-linear sequence-to-sequence e 1 1 3 1 3 5 þ þ þ ¼ pffiffiffiffiffiffiffiffi 1 þ transformation. Shanks (1955) introduced a family of 2 3 ð1!Þð8uÞ ð2!Þð8uÞ 2u ð3!Þð8uÞ non-linear transforms designated as ek, em k , ek , and ed. ð35Þ Theorems are provided to show the effectiveness of this transform both in accelerating the convergence of slowly and convergent sequences and inducing convergence in 8 9 2 2 2 4 1 ð 4 1 Þ ð 4 3 Þ divergent sequences. The Shanks transform is an iterative > > > > þ >1 > 2 = algorithm based on the directly summing series. The u < ð1!Þð8uÞ ð2!Þð8uÞ e operator of the ith order transform of Sn is defined as I1 ðuÞ ¼ pffiffiffiffiffiffiffiffi > (Shanks 1955) ð4 12 Þ ð4 32 Þ ð4 52 Þ 2u > > > > > þ ; :. ð3!Þð8uÞ3. Snþ1 Snþm. Sn. ΔS. ð36Þ ΔS ΔS n nþ1 nþm. .. .. .. ... For 0 ≤ u ≤ 2, the asymptotic expansions for K0(u) and em ðSn Þ ¼. . . . . K1(u) may be written respectively as (Abramowitz and. Stegun 1964, p 375) ΔSnþm1 ΔSnþm ΔSnþ2m1. 1 h u i 1 2 u. 1 1 1. þ = I0 ðuÞ þ 4 2 K0 ðuÞ ¼ ln. ΔS 2 ð1!Þ ΔSnþ1 ΔSnþm. n. 2 .. .. .. .. . . . . . 1 14 u2. þ 1þ 2. 2 ð2!Þ ΔSnþm1 ΔSnþm ΔSnþ2m1. 3 ð41Þ 1 1 14 u2 þ 1þ þ þ ð37Þ where ΔSn ¼ Snþ1 Sn with n=1, 2,... ∞ as the number of 2 3 ð3!Þ2 series and m=0, 1, 2,... n-2k as the number of iterations. Wynn (1956) derived a mathematical expression from and the Shanks transform. This transform for the first and h u i second iterations can then be expressed for computation as 1 þ = I1 ðuÞ þ K1 ðuÞ ¼ ln follows (. u. 8 > > > <. u 4> > > :.

(6). 2. 2. 2. 2. 2. 2. 2. u 1 2 2 1 ð4u Þ. 9 12 3 > 1 1 ð4u Þ > > ð1!Þð2!Þ þ 1 þ 2 ð2!Þð3!Þ þ 1 þ 2 þ 3 ð3!Þð4!Þ þ =

(7). 1 u2 ð 1 u2 Þ2 > > þ 1 þ 1 þ 12 ð1!4Þð2!Þ þ 1 þ 12 þ 13 ð2!4 Þð3!Þ þ > ; 1 2 4u. . ð38Þ. Hydrogeology Journal (2007) 15: 297–305. e0 ðSn Þ ¼ Sn. ð42Þ. and e1 ðSn Þ ¼ ½ΔSn 1. ð43Þ. DOI 10.1007/s10040-006-0100-x.

(8) 302. where n=1, 2,.... The third and remaining iterations are expressed as eiþ1 ðSni Þ ¼ ei1 ðSniþ1 Þ þ ½ei ðΔSni Þ1. ð44Þ. where i=1, 2,... (n – 2). The Shanks transform requires use of the even-order terms when approximating the sum of Sn; the odd-order terms are purely intermediate quantities in the numerical computations. Applying the Shanks transform to evaluate a given series requires setting a convergence criterion, ERR, defined as. e2rþ2 ðSn1 Þ e2r ðSn Þ. ERR; for n ¼ 3; 4; . . . ð45Þ. e2rþ2 ðSn1 Þ Shanks (1955) indicated that the series converges if three successive values of n satisfy the criterion. However, such a suggestion may be too strict. For comparing the implementation of the Shanks transform with the other methods to the same problem, the convergence criterion is to terminate the summing process. When applying the Shanks transform to a given series, it is defined as (Huang et al. 2000). e2i ðSn2i Þ e2ði1Þ ðSn2iþ1 Þ ERR; ð46Þ for n ¼ 3; 4; . . . where r =1, 2,... (n–2)/2 when n is even and r= 1, 2,... (n–1)/2 when n is odd.. Numerical inversion. The Laplace transforms are commonly used to solve the differential and integral equations. In many engineering problems, the Laplace–domain solutions for mathematical models are tractable, yet the corresponding solutions in the time domain may be difficult or impossible to solve. Under such circumstances, the methods of numerical Laplace inversion such as the Stehfest method (Stehfest 1970), the Crump method (Crump 1976), or the Talbot method (Talbot 1979) may be used. de Hoog et al. (1982) mentioned that the different numerical methods of inversion give accurate answers for various cases. The Laplace transform of a real-valued function f(t), t ≥ 0, is defined as Z 1 F ð pÞ ¼ exp ðpt Þf ðt Þdt ð47Þ 0. and the inversion of F(p) is given by f ðt Þ ¼. 1 2i. Z. þi1. i1. exp ðpt ÞF ð pÞdp. where κ is to the right of any singularities of F(p). Hydrogeology Journal (2007) 15: 297–305. ð48Þ. Using the trapezoidal rule with step size Tp , Crump (1976) proposed the Laplace inversion formula using an infinite series as follows 1 f ðt Þ ¼ exp ðξt Þ Tp ( ð49Þ

(9) ) 1 F ðξ Þ X ik ikt þ Re F ξ þ exp 2 Tp Tp k¼1. in which F(·) is constructed to be an infinite set of even and odd period functions, each with a period of 2Tp, the parameter is an estimate of the maximum of the real parts of singularity functions, and the symbol Re represents the real part of the function. In fact, the Crump method approximates the Laplace inversion by expressing the transform function in a Fourier series. Crump (1976) used the ε-algorithm of Wynn (1956) (also called the Shanks method) to accelerate the convergence of the sum in Eq. (49) and gave the convergence of the sequence of partial sums as

(10) in int Sn ¼ Sn1 þ Re F ξ þ exp ; Tp Tp ð50Þ n ¼ 1; 2; . . . ; 2M The Crump’s accelerated estimate of f(t) is defined as 1 b exp ðξt Þ e2M ðS0 Þ ð51Þ f ðt Þ ¼ Tp where S0 ¼ F ðξ Þ=2 and the transform ej(Sn) is defined as Eq. (44). de Hoog et al. (1982) applied the Pade approximations to improve the acceleration procedure to the transform sequence in Eq. (51). The Pade approximations represent a function as the quotient of two polynomials a0 þ a1 t þ a2 t2 þ þ an tn ; 1 þ b1 t þ b2 t 2 þ þ bm t m N ¼nþm. f ðt Þ ffi RN ðt Þ ¼. ð52Þ. where the parameters, a0, a1,... and b1, b2,... bm, are available for the approximation of f(t), by RN(t). The Pade approximations are related to Maclaurin expansions in that the coefficients are determined in a similar fashion to make bf ðt Þ and RN(t) agree at t = 0 and the first N derivatives agree at t = 0 (Gerald and Wheatley 1989). de Hoog et al. (1982) also applied the Shanks method to accelerate the running sums of infinite series. The dimensionless hydraulic head and water level of a well are determined by inverting Eqs. (15), (16), and (21) with an appropriate numerical algorithm (e.g., Stehfest 1970; Crump 1976; or Talbot 1979). Herein, the modified Crump algorithm (Crump 1976; de Hoog et al. 1982) is adopted to invert the Laplace-domain solutions, Eqs. (15), (16), and (21), to an accuracy of five decimal places. This routine was developed based on an algorithm originally DOI 10.1007/s10040-006-0100-x.

(11) 303. proposed by Crump (Crump 1976) and modified by de Hoog et al. (de Hoog et al. 1982). A double-precision format is used for all evaluations and the convergence criterion for the Shanks method is set as 10−7. Then, the numerical inversion results are compared with the values obtained from those given by Yang and Gates (1997). The effect of dimensionless hydraulic head distribution is also discussed for a composite aquifer system.. Numerical results and discussion Based on a composite (two-zone) aquifer model, Yang and Gates (1997) used a finite-element method to simulate the dimensionless head distribution for a slug test. In their simulations, the hydraulic conductivity and storativity are respectively 1.8×10−6 cm/sec and 0.01 for a positive skin, 1.8×10−5 cm/sec and 0.01 for an undisturbed formation, and 1.8×10−4 cm/sec and 0.1 for a negative skin. Initially, MHS is evaluated using a numerical inversion routine INLAP (IMSL 1997). Table 1 gives the values of dimensionless heads in a tested well calculated using the parameters provided by Yang and Gates (1997) for rc=5.08 cm, rw=9.15 cm, and rs=30.5 cm when α (conductivity ratio)=0.1, 1, or 10. Note that the formation. has a negative skin if α=0.1 and a positive skin if α=10; when the formation has only a single zone, α=1. The differences of the tabular values for dimensionless heads between MHS and Yang and Gates simulation results are significant for all cases, especially for a large t and/or a small α. Figure 1 shows the curves of dimensionless head versus time obtained by MHS and those given by Yang and Gates (1997) when α=0.1, 1, or 10. Obviously, the curves of Yang and Gates (1997) are significantly different from those of MHS. In fact, their mass balance equation— Yang and Gates (1997); Eq. (4)—and interface condition— Yang and Gates (1997), Eq. (5)—are questionable. The LHS term of Eq. (5) is set r=rw in this paper to represent the mass balance equation at the wellbore. However, Yang and Gates (1997) neglected the well radius of the test well; in other words, they treated the well as a line source. Consequently, the estimated hydraulic heads near the well may have significant errors. In addition, the partial derivative of dimensionless hydraulic head in Eq. (8) was taken with respect to dimensionless time t, instead of r, as used in Yang and Gates (1997). Physically, Eq. (8) should be used to preserve the continuity of the flow rate at the interface of the skin and formation zones. The dimensionless well water levels estimated from MHS are compared with those given by Yang and Gates. Table 1 Dimensionless water level in a well, estimated by MHS and given in Yang and Gates (1997) for rc=5.08 cm, rw=9.15 cm, and rs=30.5 cm when α=0.1, 1, or 10 Time (sec) 2 4 6 8 10 15 20 25 30 35 40 45 50 60 70 80 90 100 140 180 220 260 300 400 500 600 700 800 900 1000. α=0.1 MH. YG. α=1 MH. YG. α=10 MH. YG. 0.87666 0.82822 0.79230 0.76289 0.73767 0.68617 0.64510 0.61070 0.58102 0.55490 0.53160 0.51057 0.49144 0.45780 0.42908 0.40428 0.38268 0.36377 0.30820 0.27395 0.25222 0.23795 0.22819 0.21380 0.20551 0.19946 0.19445 0.19003 0.18601 0.18230. 0.8179 0.6905 0.5979 0.5281 0.4738 0.3799 0.3195 0.2772 0.2457 0.2211 0.2014 0.1852 0.1716 0.1500 0.1336 0.1207 0.1103 0.1019 0.0802 0.0694 0.0636 0.0603 0.0583 0.0555 0.0539 0.0526 0.0515 0.0506 0.0497 0.0489. 0.98649 0.98046 0.97567 0.97153 0.96781 0.95966 0.95257 0.94616 0.94025 0.93472 0.92949 0.92452 0.91976 0.91076 0.90234 0.89438 0.88681 0.87956 0.85314 0.82971 0.80843 0.78881 0.77053 0.72940 0.69325 0.66088 0.63153 0.60468 0.57993 0.55701. 0.9801 0.9632 0.9484 0.9351 0.9230 0.8966 0.8739 0.8538 0.8356 0.8188 0.8032 0.7887 0.7750 0.7497 0.7268 0.7057 0.6862 0.6680 0.6056 0.5550 0.5126 0.4763 0.4447 0.3811 0.3328 0.2949 0.2644 0.2395 0.2187 0.2012. 0.99590 0.99415 0.99279 0.99164 0.99061 0.98839 0.98649 0.98480 0.98325 0.98181 0.98046 0.97918 0.97796 0.97567 0.97354 0.97153 0.96962 0.96781 0.96119 0.95531 0.94994 0.94495 0.94025 0.92949 0.91975 0.91075 0.90230 0.89429 0.88664 0.87928. 0.9978 0.9957 0.9937 0.9916 0.9896 0.9848 0.9802 0.9757 0.9714 0.9673 0.9633 0.9594 0.9556 0.9485 0.9417 0.9352 0.9291 0.9232 0.9017 0.8828 0.8657 0.8501 0.8356 0.8033 0.7750 0.7497 0.7268 0.7057 0.6861 0.6678. rc=radius of well casing; rw=well radius; rs is the distance from center of well to outer boundary of skin zone. α=K1/K2, where K1 is the hydraulic conductivity of the skin, and K2 is the hydraulic conductivity of the formation. MH and YG respectively denote the results obtained from MHS and given in Yang and Gates (1997). The hydraulic conductivity and storativity are, respectively, 1.8×10−6 cm/sec and 0.01 for a positive skin, 1.8×10−5 cm/sec and 0.01 for an undisturbed formation, and 1.8×10−4 cm/sec and 0.1 for a negative skin Hydrogeology Journal (2007) 15: 297–305. DOI 10.1007/s10040-006-0100-x.

(12) 304. Fig. 1 Curves of dimensionless head versus time obtained by MHS and those given by Yang and Gates (1997) when α=0.1, 1 or 10. (1997) based on a finite-element modeling. The result indicates that the Yang and Gates tabular values have very poor accuracy. This problem is mainly attributed to the problem that the interface condition for the flow rate in their mathematical model is incorrect. Figure 2 displays the curves of dimensionless hydraulic head versus dimensionless distance for rc=5.08 cm, rw=9.15 cm, rs=30.5 cm, and t=103 sec when α=0.1, 1 or 10. The storativity of the skin and formation zones is assumed as the same to simplify the problems and is set as. 0.01. For the case without a skin zone, the dimensionless hydraulic head gradually decreases with increasing radial distance as shown in Fig. 2. If a finite-thickness skin is present, the relation of dimensionless hydraulic head versus dimensionless distance exhibits two 0 curves with different slope joined at the interface rs ¼ 3:33 . A negative skin, which has a higher transmissivity than the formation, has a curve with relative mild slope in the skin zone and with steeper slope in the formation zone. In contrast, a positive skin has a very steep slope in the skin zone due to the lower transmissivity and a relative flat slope in the formation zone. In addition, the difference of dimensionless hydraulic head between the skin and formation zones for an aquifer with a positive skin is larger obviously than that with a negative skin. This indicates that the dimensionless hydraulic head distribution depends on the hydraulic properties of both the skin and formation zones. If ignoring the skin zone, the dimensionless hydraulic head in a composite aquifer system may be apparently overestimated or underestimated and the error of the test-data analysis will be inducted.. Conclusions A detailed derivation has been presented as proof that MHS for a slug test in a well surrounded by a finitethickness skin satisfies the governing equations, boundary conditions, and continuity equations for hydraulic head and flow. Two new formulae relating the modified zeroorder Bessel functions of the first and second kinds to their first and second derivatives are established. An efficient numerical approach is proposed to evaluate MHS. The numerical values of finite-element simulation given in Yang and Gates (1997) are apparently different from those evaluated based on the MHS. The problems in the mass balance equation and interface condition presented by Yang and Gates (1997) make their simulation results differ from the values of MHS. The plots for the spatial and temporal distributions of dimensionless hydraulic head demonstrate the dependence of the head on the hydraulic properties of both the skin and formation zones. In addition, the effect for an aquifer with a positive skin is significantly larger than that with a negative skin.. Notation b h h′. Fig. 2 Curves of dimensionless hydraulic head versus dimensionless distance obtained by MHS for rc=5.08 cm, rw=9.15 cm, rs=30.5 cm, and t= 103 sec when α=0.1, 1 or 10 Hydrogeology Journal (2007) 15: 297–305. Thickness of aquifer Hydraulic head in the skin or formation zone Dimensionless hydraulic head in the skin or formation zone hw Water level in a well h'w Dimensionless water level in a well I 0 (·) , Modified Bessel function of the first and K0(·) second kinds and order zero, respectively I 1 (·) , Modified Bessel function of the first and K1(·) second kinds and order one, respectively DOI 10.1007/s10040-006-0100-x.

(13) 305. K P q r rc rs rw r' r's Ss t t' α β γ δ. Hydraulic conductivity Laplace transform variable of dimensionless time =p1/2 Radial distance from center of well Standpipe radius Distance from center of well to outer boundary of skin Test well radius =r/rw Dimensionless radial distance from center of well =rs/rw Dimensionless distance from center of well to outer boundary of skin Specific storage Test time =K2t/(rw)2Ss2 =K2/K1 =(δα)1/2 =rc2/(2rw2Ss2b) =Ss1/Ss2. Acknowledgements This study was partly supported by the Taiwan National Science Council under the grant NSC 91-2211-E-009-020. The authors appreciate the comments and suggested revisions of anonymous reviewers that helped to improve the clarity of our presentation.. References Abramowitz M, Stegun IA (1964) Handbook of mathematical functions with formulas, graphs and mathematical tables, National Bureau of Standards. Dover, Washington, DC Butler Jr JJ, Healey JM (1998) Relationship between pumping-test and slug-test parameters: scale effect or artifact?. Ground Water 36(2):305–313 Cooper Jr HH, Bredehoeft JD, Papadopuls IS (1967) Response of a finite-diameter well to an instantaneous charge of water. Water Resour Res 3(1):263–269 Crump KS (1976) Numerical inversion of Laplace transforms using a Fourier series approximation. J Assoc Comput Mach 23(1): 89–96 de Hoog FR, Knight JH, Stokes AN (1982) An improved method for numerical inversion of Laplace transforms. Society for Industrial and Applied Mathematics. J Sci Stat Comput 3(3): 357–366 Faust CR, Mercer JW (1984) Evaluation of slug tests in wells containing a finite-thickness skin. Water Resour Res 20(4):504– 506 Gerald CF, Wheatley PO (1989) Applied Numerical Analysis, 5th edn. Addison-Wesley, Boston, MA. Hydrogeology Journal (2007) 15: 297–305. Hemker CJ (1999) Transient well flow in layered aquifer systems: the uniform well-face drawdown solution. J Hydrol 225: 19–44 Huang CP, Yeh HD, Yang SY (2000) Applications of accelerate methods on the evaluation of two types of drawdown solutions. Proceedings of the 11th hydraulic engineering conference, National Taiwan University, Taipei, Taiwan, pp D147–D152 IMSL (1997) Math/Library, vols. 1 and 2, Visual Numerics, Houston, TX Kahaner DK (1972) Numerical quadrature by the ɛ-algorithm. Math Comp 26:689–693 Karasaki K (1990) A systematized drillstem test. Water Resour Res 26(12):2913–2919 Karasaki K, Long JCS, Witherspoon PA (1988) Analytical models of slug tests. Water Resour Res 24(1):115–126 McElwee CD (2002) Improving the analysis of slug tests. J Hydro 269:122–133 McLachlan NW (1955) Bessel Functions for Engineers, 2nd edn. Oxford University Press, London Moench AF, Hsieh PA (1985a) Comment on “Evaluation of slug tests in wells containing a finite-thickness skin” by C.R. Faust and J.W. Mercer. Water Resour Res 21(9):1459–1461 Moench AF, Hsieh PA (1985b) Analysis of slug test data in a well with finite thickness skin. The IAH 17th International congress on Hydrology of rocks of low permeability, Int Assoc of Hydrogeol. Tucson, AZ, 7–12 Jan 1985 Ramey Jr HJ, Agarwal RG, Martin I (1975) Analysis of “Slug Test” or DST flow period data. J Can Pet Technol July/Sept:37–47 Sageev A (1986) Slug test analysis. Water Resour Res 22(8):1323– 1333 Shanks D (1955) Non-linear transformations of divergent and slowly convergent sequences. J Math Phys 34:1–42 Stehfest H (1970) Numerical inversion of Laplace transforms. Commun ACM 13(1):47–49 Talbot A (1979) The accurate numerical inversion of Laplace transforms. J Inst Math Appl 23:97–120 Watson GN (1958) A treatise on the theory of Bessel functions, 2nd edn. Cambridge University Press, Cambridge Wynn P (1956) On a device for computing the em(Sn) transformation. Math Tables Other Aids Comp 10:91–96 Yang YJ, Gates TM (1997) Wellbore skin effect in slug-test data analysis for low-permeability geologic materials. Ground Water 35(6):931–937 Yang SY, Yeh HD (2005) Laplace-domain solutions for radial twozone flow equations under the conditions of constant-head and partially penetrating well. J Hydrau Eng ASCE 131(3):209– 216 Yeh HD, Yang SY, Peng HY (2003) A new closed-form solution for radial two-layer drawdown equation under constant-flux pumping in a finite-radius well. Adv Water Resour 26(5):747– 757 Zlotnik VA, McGuire VL (1998a) Multi-level slug test in highly permeable formations: 1. Modification of the Springer-Gelhar (SG) model. J Hydrol 204:271–282 Zlotnik VA, McGuire VL (1998b) Multi-level slug test in highly permeable formations: 2. Hydraulic conductivity identification, method verification, and field applications. J Hydrol 204:283– 296. DOI 10.1007/s10040-006-0100-x.

(14)