Transient Flow into a Partially Penetrating Well during
the Constant-Head Test in Unconfined Aquifers
Ya-Chi Chang
1; Geng-Yuan Chen
2; and Hund-Der Yeh, Aff.M.ASCE
3Abstract: This study derives a semianalytical solution for drawdown distribution during a constant-head test at a partially penetrating well in an unconfined aquifer. The constant-head condition is used to describe the boundary along the screen. In addition, a free-surface condition is used to delineate the upper boundary of the unconfined aquifer. The Laplace-domain solution is then derived using separation of variables and Laplace transform. This solution can be used to identify the aquifer parameters from the data of the constant-head test when integrated with an optimization scheme or to investigate the effects of vertical flow caused by the partially penetrating well and free-surface boundary in an unconfined aquifer.DOI: 10.1061/(ASCE)HY.1943-7900.0000392. © 2011 American Society of Civil Engineers.
CE Database subject headings: Aquifers; Analytical techniques; Wells; Transient flow.
Author keywords: Unconfined aquifer; Semianalytical solution; Partially penetrating well; constant-head test.
Introduction
Aquifer tests, such as constant-head test (CHT) and constant-flux test (CFT), are usually performed to estimate aquifer parameters such as specific storage and hydraulic conductivity. For aquifers with low transmissivity, CHT is more suitable to apply than CFT. The wellbore storage at the pumping well has a large effect on the early drawdown behavior at pumping and observation wells in CFT (Renard 2005). If a CHT is established in a short period of time, the effect of wellbore storage is negligible if the aquifer has low transmissivity and the well radius is small (Chen and Chang 2003).
Many studies have been devoted to the solutions for CHT. Kirkham (1959) derived a steady-state solution for ground-water distribution in a bounded confined aquifer pumped by a partially penetrating well under CHT. They simplified the complexity of the geometry by dividing the model into two different regions. Javan-del and Zaghi (1975) considered the ground water in a confined aquifer pumped by a fully penetrating well that is radially extended at the bottom of the aquifer. The procedure used in their study is similar to that in Kirkham (1959), and the steady-state ground-water solution was obtained by separation of variables. Jones et al. (1992) and Jones (1993) discussed the practicality of CHTs on wells completed in low-conductivity glacial till deposits. Mishra and Guyonnet (1992) indicated the operational benefit of CHTs when the total available drawdown is limited by well construction and aquifer characteristics. They developed a method for analyzing observation well response under CHT. Issues involving CHT can be
found in the literature (e.g.,Uraiet and Raghaven 1980;Chen and Chang 2003; Yeh and Yang 2006; Singh 2007; Wang and Yeh 2008).
Considering a CHT performed in a partially penetrating well, Yang and Yeh (2005) developed a time-domain solution to describe the drawdown in a confined aquifer with a finite-thickness skin. The boundary conditions along the partially penetrating well are represented by a constant-head (first kind) boundary for the screen and a no-flow (second kind) boundary for the casing. They trans-formed the first-kind boundary along the screen into a second-kind boundary with an unknown flux that is time dependent; therefore, the boundary along the partially penetrating well became uniform. The solution was then solved by the Laplace and finite Fourier co-sine transforms. Chang and Yeh (2009) used the methods of dual-series equations and perturbation method to solve the mixed boun-dary problem for the CHT at a partially penetrating well. Chang and Yeh (2010) further developed an analytical solution for a partially penetrating well with arbitrary location of the well screen under constant-head tests in confined aquifers. However, the aforemen-tioned studies are only applicable for confined aquifers.
For unconfined aquifers, Chen and Chang (2003) developed a well hydraulic theory for CHT performed in a fully penetrating well and established a parameter estimation method. Chang et al. (2010) extended the work of Yang and Yeh (2005) to develop a mathemati-cal model for an unconfined aquifer system while treating the skin as a finite-thickness zone and derived the associated solution for CHT at a partially penetrating well. For other environmental appli-cations, light non-aqueous-phase liquids (LNAPLs) are usually re-covered by wells held at constant drawdown (Abdul 1992;
Murdoch and Franco 1994), and constant-head pumping is used to control off-site migration of contaminated ground water (Hiller and Levy 1994). At LNAPL contaminant sites, the pollutant forms a pool of LNAPL in the subsurface on top of the water table. In installing a well in unconfined aquifers, therefore, the screen goes from the top of the aquifer.
For CFT in unconfined aquifers, Neuman (1972) presented a new analytical solution for characterizing flow to a fully penetrat-ing well in an unconfined aquifer. He assumed that the drainage above the water table occurs instantaneously. Accounting for the effect of a finite-diameter pumping well, Moench (1997) developed a solution in Laplace domain for the flow to a partially penetrating 1Postdoctorate, Institute of Environmental Engineering, National Chiao
Tung Univ., Hsinchu, Taiwan. E-mail: rachel.ev91g@nctu.edu.tw
2M.S. Student, Institute of Environmental Engineering, National Chiao
Tung Univ., Hsinchu, Taiwan. E-mail: sunny0586.ev97g@nctu.edu.tw
3Professor, Institute of Environmental Engineering, National Chiao
Tung Univ., Hsinchu, Taiwan (corresponding author). E-mail: hdyeh@ mail.nctu.edu.tw
Note. This manuscript was submitted on March 9, 2010; approved on December 27, 2010; published online on December 29, 2010. Discussion period open until February 1, 2012; separate discussions must be submitted for individual papers. This paper is part of the Journal of Hydraulic En-gineering, Vol. 137, No. 9, September 1, 2011. ©ASCE, ISSN 0733-9429/ 2011/9-1054–1063/$25.00.
well in unconfined aquifers. Contrary to Neuman’s assumption, Moench used the free-surface boundary in Boulton (1955), assum-ing that the drainage of pores occurs as an exponential function of time in response to a step change in hydraulic head in the aquifer. Tartakovsky and Neuman (2007) presented an analytical solution for drawdown in an unconfined aquifer caused by pumping at a constant rate from a partially penetrating well. They generalized the solution of Neuman (1972,1974) by accounting for unsaturated flow above the water table and derived the solution from a linear-ized Richards’ equation in which unsaturated hydraulic conduc-tivity and water content are expressed as exponential functions of incremental capillary pressure head relative to its air entry value. Motivated by the aforementioned research, this paper aims to develop a mathematical model for constant-head pumping from a partially penetrating well in an unconfined aquifer. Without as-suming constant-head boundary along the screen as an unknown flux boundary, the system is separated into two different regions, and the solutions are directly obtained by separation of variables and Laplace transform. This new solution can be used to determine the aquifer parameters or to investigate the effects of vertical flow caused by the partially penetrating well and free-surface boundary on the drawdown distribution in unconfined aquifers.
Mathematical Model Laplace-Domain Solutions
The conceptual model for constant-head pumping in an unconfined aquifer system with a partially penetrating well is illustrated in Fig.1. The well screen starts from z¼ zlwith a finite well radius
rw, and the bottom of the screen is sealed. The domain is divided
into two different regions. Region 1 is defined by 0≤ r ≤ rw and
0≤ z ≤ zl whereas Region 2 is bounded within rw≤ r < ∞ and
0≤ z ≤ η, where η = saturated thickness. The aquifer is assumed to be homogeneous, with infinite extent in the radial direction, and the seepage face in Region 2 is neglected. Under this assumption, the governing equations in terms of drawdown in Regions 1 and 2 can, respectively, be written as
Kr ∂2s 1 ∂r2 þ 1 r ∂s1 ∂r þ Kz∂ 2s 1 ∂z2 ¼ Ss∂s 1 ∂t ; 0≤ r ≤ rw; 0≤ z ≤ zl ð1Þ and Kr ∂2s 2 ∂r2 þ 1 r ∂s2 ∂r þ Kz∂ 2s 2 ∂z2 ¼ Ss ∂s2 ∂t ; rw≤ r < ∞; 0≤ z ≤ η ð2Þ
The subscripts 1 and 2 denote Regions 1 and 2, respectively. The drawdown at distance r from the center of the well and distance z from the bottom of the aquifer at time t is denoted as sðr; z; tÞ, which is equal to h0 h, where h0 and h = initial and hydraulic
head, respectively. The aquifer has the horizontal hydraulic conduc-tivity Kr, vertical hydraulic conductivity Kz, specific storage Ss, and
specific yield Sy. Assuming that the drawdown is small compared
with the saturated aquifer thicknessη, the boundary at the free sur-face (z¼ η) can be approximated as z ¼ b, where b = initial satu-rated thickness. Therefore, the governing equation for Region 2 can be further expressed as Kr ∂2s 2 ∂r2 þ 1 r ∂s2 ∂r þ Kz ∂2s 2 ∂z2 ¼ Ss ∂s2 ∂t ; rw≤ r < ∞; 0≤ z ≤ b ð3Þ
The initial condition for saturated thicknessηðr; tÞ is equal to b; therefore, the drawdowns are assumed to be zero initially in Regions 1 and 2, that is,
s1ðr; z; 0Þ ¼ s2ðr; z; 0Þ ¼ 0 ð4Þ
The no-flow boundary condition at the bottom of the aquifer for both regions is ∂s1ðr; z; tÞ ∂z z¼0¼ ∂s2ðr; z; tÞ ∂z z¼0¼ 0 ð5Þ
The boundary at the top of the Region 1 can also be expressed as ∂s1ðr; z; tÞ ∂z z¼z l ¼ 0; 0< r < rw ð6Þ
The nonlinear boundary describing the free surface in Region 2 for the unconfined aquifer can be linearized to the form (Neuman 1972) Kz∂s2ðr; z; tÞ ∂z z¼b¼ Sy∂s2ðr; z; tÞ ∂t z¼b ð7Þ
In addition, the boundary at r¼ 0 attributable to symmetry along the center of the well is written as
∂s1ðr; z; tÞ
∂r
r¼0¼ 0; 0< z < z1 ð8Þ
When r approaches infinity, the boundary condition for Region 2 is
s2ð∞; z; tÞ ¼ 0 ð9Þ
The boundary condition specified along the well is
s2ðrw; z; tÞ ¼ sw; z1< z < b; t> 0 ð10Þ
where sw= constant drawdown in the well at any time.
At the interface between Regions 1 and 2, the continuities of the drawdown and flow rate must be satisfied:
s1ðrw; z; tÞ ¼ s2ðrw; z; tÞ; 0< z < zl; t> 0 ð11Þ
Fig. 1. Schematic of constant-head test in unconfined aquifer with par-tially penetrating well
and ∂s1ðr; z; tÞ ∂r r¼rw ¼ ∂s2ðr; z; tÞ ∂r r¼rw ; 0< z < zl; t> 0 ð12Þ To express the solutions in dimensionless form, the following dimensionless variables are defined: s1¼ s1=sw; s2¼ s2=sw;
σ ¼ Sy=Ssb; κ ¼ Kz=Kr; ρ ¼ r=rw; ρw¼ rw=b; αw¼ κρ2w;
α ¼ αwρ2; τ ¼ Krt=Ssr2w; ζ ¼ z=b; and ζl¼ zl=b, where s1 and
s2= dimensionless drawdowns for Regions 1 and 2, respectively; σ = ratio of specific yield Sy to the storativity Ssb;κ represents
the dimensionless conductivity ratio;ρ denotes the dimensionless radial distance;ρw = dimensionless radius of the pumping well;
α = dimensionless conductivity ratio times the square of the ratio of radial distance r from pumping well to aquifer thickness b; τ refers to the dimensionless time during the test; andζ and ζl=
di-mensionless vertical distance and the didi-mensionless distance from the bottom of aquifer to the bottom of the screen, respectively.
Taking the Laplace transform to the dimensionless governing equations of Eqs. (17) and (18) subject to the dimensionless boun-dary conditions of Eqs. (20)–(27), the Laplace-domain solutions for the dimensionless drawdowns in Regions 1 and 2 are, respectively, as follows:
Fig. 2. Dimensionless drawdown distributions at (a)τ ¼ 1, (b) τ ¼ 102, (c)τ ¼ 104, and (d)τ ¼ 106
~s1ðρ; ζ; pÞ ¼ X∞ m¼0 A01mðpÞ I0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω1m p ρÞ I0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω1m p Þ cosðΩ1mζÞ; t> 0 ð13Þ and ~s2ðρ; ζ; pÞ ¼ X∞ n¼0 A02nðpÞK0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω2n p ρÞ K0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω2n p Þ cosðΩ2nζÞ; t> 0 ð14Þ Applying the continuity conditions to Eqs. (13) and (14), the coefficients A01m and A02n are respectively obtained as
A01mðpÞ ¼ iom Ω1m cosðζlΩ1mÞ sinðζlΩ1mÞ þ ζlΩ1m ·X ∞ n¼0 A02nðpÞk0nΛmn ð15Þ and A02nðpÞ ¼ 2Ω2n sinð2Ω2nÞ þ 2Ω2n ·X ∞ m¼0 A01mðpÞΛmn þ4 p sinðΩ2nÞ sinðΩ2nζlÞ sinð2Ω2nÞ þ 2Ω2n ð16Þ where p = Laplace variable; and A01m, A02n,Ω1m,Ω2n,Λmn, I0ð·Þ,
I1ð·Þ, K0ð·Þ, K1ð·Þ, and k0are defined in the notation list. The
de-tailed derivations of Eqs. (13) and (14) are given in the appendix. Fully Penetrating Wells in Unconfined and Confined Aquifers
By settingζl¼ 0 in Eqs. (15) and (16), the drawdown solution of
Eq. (13) in Region 1 is equal to zero, and the Laplace-domain sol-ution in Eq. (14) for dimensionless drawdown in Region 2 with fully penetrating wells in unconfined aquifers is exactly the same as the solution given in Chen and Chang [2003, Eq. (7)] when the skin factor Sk equals zero after some algebraic manipulations.
Furthermore, by settingζl¼ 0 and σ ¼ 0, the Laplace-domain
sol-ution of Eq. (14) in Region 2 can be reduced to the solution in Hantush (1964) for drawdown with a fully penetrating well in con-fined aquifers.
Results and Discussion
The numerical inversion method given by Stehfest (1970) is adopted for calculating the dimensionless drawdown solutions in Eqs. (13) and (14) for Regions 1 and 2, respectively, in real-time domain. Because there may be nonconvergence issues when evalu-ating the infinite summations in Eqs. (13) and (14), the Shanks method is applied to accelerate convergence for these infinite sum-mations. This method has been successfully applied to compute the solutions arising in the ground-water area (e.g., Yang and Yeh 2002;Peng et al. 2002).
Fig. 2(a)demonstrates the dimensionless drawdown distribu-tions for the dimensionless distanceρ ¼ 1, 1.1, and 1.5 at the di-mensionless timeτ ¼ 1, Fig.2(b)forρ ¼ 1, 2, and 5 at τ ¼ 102, Fig.2(c)forρ ¼ 1, 2, and 7 at τ ¼ 104, and Fig.2(d)forρ ¼ 1, 5,
and 7 atτ ¼ 106. The aquifer parameters used in these figures are
as follows:κ ¼ 1, σ ¼ 103, andζl¼ 0:5. These figures show that
the dimensionless drawdown atρ ¼ 1 matches the boundary con-dition of the wellbore at different time periods. The dimensionless drawdown decreases with increasingρ at τ ¼ 1, 102, 104, and 106.
In addition, it is apparent that vertical flows occur at the water table because of the free-surface boundary, as shown in Figs.2(b)–2(d). Fig.3is plotted to examine the effectσ (i.e., Sy=Ssb) of Region 2
on the dimensionless drawdown during CHT. This figure shows the response of dimensionless drawdown in a 100-m-thick aquifer at ρ ¼ 50, κ ¼ 1, ζ ¼ 0:75, and ζl¼ 0:5 for σ ranging from 0 to
3 × 103. The dimensionless drawdown decreases with increasing
σ. The typical three-stage drawdown patterns can be observed. The water releases from the elastic behavior of the aquifer forma-tion at an early time, i.e., the first stage. During the second stage at moderate times, the gravity drainage almost stabilizes the water ta-ble. Finally, the effect of vertical flow vanishes at late times, and the flow behaves like the first stage again. Fig.3shows that the largerσ is, the longer the delayed yield stage will be, perhaps because a largerσ supply more water from the drainage. If σ ¼ 0, the top boundary represented by Eq. (7) becomes the no-flow condition, and the aquifer can therefore be considered as confined.
Fig. 4 illustrates that the distributions of the dimensionless drawdown at the well screen extends from ζ ¼ ζl to ζ ¼ 1 in
Region 2 whenτ ¼ 106. This figure shows that the dimensionless
drawdown increases with the length of well screen. In addition, large slopes of the drawdown distribution curves occur near the free-surface boundary and the edge of the screen. Therefore, ver-tical ground-water flows are obviously large at these two areas.
The response of dimensionless drawdown versus dimensionless time at different observed locations is plotted in Fig.5forρ ¼ 10 and 100 with κ ¼ 1 and ζl¼ 0:25. The σ is zero for a confined
aquifer, and there is no vertical flow [Fig.5(a)]. On the other hand, the vertical flow is apparent at moderate times for different radial distances [Fig.5(b)] whenσ ¼ 103. Figs.6(a)and6(b)illustrates
the spatial flow pattern forσ ¼ 0 and 103atτ ¼ 104with the same
parameter values as those in Fig.5. Apparently, the vertical flow occurs only near the bottom edge of the well screen when the aqui-fer is confined. However, for unconfined aquiaqui-fers, the flow at free surface is almost vertical, and obvious vertical flows occur near both the top and bottom edges of the well. It demonstrates that the vertical flow in the unconfined system is induced not only by the effect of partial penetration but also the effect of free-surface boundary.
Fig. 3. Effect ofσ on dimensionless drawdown during CHT
Fig. 7 demonstrates the effect of the conductivity ratio κð¼ Kz=KrÞ on the dimensionless drawdown during CHT. The
ver-tical axis represents the dimensionless drawdown, and the horizon-tal axis represents the dimensionless time. Theκ ranges from 102 to 1 with Kr¼ 104m= min, σ ¼ 103, andζl¼ 0:5. The
dimen-sionless drawdown decreases with increasingκ, indicating that the vertical flow from delayed gravity drainage becomes large for greaterκ. Fig. 8illustrates the effect of the well radius on draw-down distribution in a 10-m-thick aquifer. The consideredρwwell
radii are 0.1, 0.01, and 0.001 m withσ ¼ 103,ζ ¼ 0:75, ζ l¼ 0:5,
andκ ¼ 1. Drawdown is calculated at the dimensionless distances of 3.16, 10, or 31.6 from the pumping well forα ¼ 101, 1, and 101, respectively. The drawdown decreases with increasing dis-tance from the pumping well for different ρw, as demonstrated
in Fig. 2. The drawdown increases with ρw for different values
ofα, indicating that the well radius has significant effect on the drawdown distribution. The effect ofα on drawdown in the aquifer atζ ¼ 0 when the well is fully (ζl¼ 0) and partially penetrating
(ζl¼ 0:8) is plotted in Fig.9for σ ¼ 103 and κ ¼ 1. The
draw-down difference between the cases of full penetration and partial penetration decreases with increasing α. It is reasonable that α1=2is directly proportional to the radial distance from the pumping
well when the aquifer is isotropic, and the partial penetration effect vanishes when the radial distance becomes large. Because r¼ α1=2b=pffiffiffiκ, this proves that the radial distance influenced by the partial penetration in an unconfined aquifer under CHT is pro-portional to the aquifer thickness, as do the results from Hantush (1964) for a confined aquifer under CFT.
The error of estimated dimensionless drawdown along the screen calculated from Eqs. (13) and (14) for different number of terms of the infinite series and Shanks method is demonstrated in Fig.10. As illustrated in this figure, the error decreases with in-creasing number of terms used in calculating the infinite series in the solution, and the largest error occurs at the edges of the screen. The largest errors are 0:179 when n ¼ m ¼ 100; 9:72 × 102 when n¼ m ¼ 200; 4:22 × 102 when n¼ m ¼ 500; 2:02 × 102 when n¼ m ¼ 1;000; and 5:64 × 104 when applying Shanks method. On the other hand, the smallest errors are9:61 × 103 when n¼ m ¼ 100; 4:61 × 103 when n¼ m ¼ 200; 1:91 × 104 when n¼ m ¼ 500; 9:35 × 104 when n¼ m ¼
1;000; and 6:18 × 1010 when applying Shanks method. These results indicate that the Shanks method can be applied to effectively accelerate convergence for the infinite summations in dimension-less drawdown in Eqs. (13) and (14).
Fig. 5. Relationship for dimensionless drawdown versus dimensionless time with ζ ¼ 0:5, 0.75, and 1.0 at ρ ¼ 10 or 100 for (a) σ ¼ 0 and
(b)σ ¼ 103
Fig. 4. Dimensionless drawdown distributions at well screen extended
fromζ ¼ ζlto ζ ¼ 1 in Region 2
Concluding Remarks
A semianalytical solution of the drawdown distribution is devel-oped for CHT performed in an unconfined aquifer with a partially penetrating well. The Laplace transforms and the method of sep-aration of variables are used to derive the transient drawdown in the Laplace domain for CHT. The Stehfest method is used to invert the
solutions in time domain, and the Shanks method is applied to ac-celerate convergence in evaluating the infinite summations in the solution.
Large slopes of the drawdown distribution curves can be ob-served near the free-surface boundary and the edge of the screen, which indicates that the vertical ground-water flows occur at these two areas. The dimensionless drawdown decreases with increasing σ but increases with the length of well screen. For different ρw, the
drawdown decreases with the increase of radial distance from the pumping well, and it may produce a large error in drawdown if the radius of the pumping well is assumed infinitesimal.
Fig. 6. Spatial flow pattern in unconfined aquifer with partially penetrating well forκ ¼ 1, ζl¼ 0:25 at τ ¼ 104when (a)σ ¼ 0 and (b) σ ¼ 103
Fig. 7. Effect of conductivity ratio (κ) of Region 2 on dimensionless
drawdown during CHT
Fig. 8. Drawdown distribution for well with three different values of
dimensionless well radii (ρw¼ 0:1, 0.01, and 0.001) with σ ¼ 103,
ζ ¼ 0:75, ζl¼ 0:5, and κ ¼ 1
The present solution can serve as an invaluable tool to explore the effects of hydraulic parameters on flow behavior in unconfined aquifers. In addition, it can reduce the solution for a fully penetrat-ing well in confined or unconfined aquifers under CHT.
Appendix
The dimensionless governing equations of Eqs. (1) and (3) can be expressed as ∂2s 1 ∂ρ2 þ 1 ρ ∂s 1 ∂ρþ αw∂ 2s 1 ∂ζ2 ¼ ∂ s1 ∂τ ; 0≤ ρ ≤ 1; 0≤ ζ ≤ ζl ð17Þ and ∂2s 2 ∂ρ2 þ 1 ρ ∂s 2 ∂ρ þ αw ∂2s 2 ∂ζ2 ¼ ∂ s2 ∂τ ; 1≤ ρ < ∞; 0≤ ζ ≤ 1 ð18Þ
The dimensionless initial conditions for Regions 1 and 2 are s1ðρ; ζ; 0Þ ¼ s2ðρ; ζ; 0Þ ¼ 0 ð19Þ and the boundary conditions at the bottom and top of the aquifer for Regions 1 and 2 in terms of dimensionless form can be written as
∂s 1ðρ; ζ; τÞ ∂ζ ζ¼0¼ ∂s2ðρ; ζ; τÞ ∂ζ ζ¼0¼ 0 ð20Þ ∂s 1ðρ; ζ; τÞ ∂ζ ζ¼ζ l ¼ 0; 0< ρ < 1 ð21Þ and ∂s 2ðρ; ζ; τÞ ∂ζ ζ¼1 ¼ σ αw ∂s 2ðρ; ζ; τÞ ∂τ ζ¼1; 1≤ ρ < ∞ ð22Þ The dimensionless boundary conditions atρ ¼ 0 and infinity are respectively written as
∂s 1ð0; ζ; τÞ
∂ρ ¼ 0; 0< ζ < ζl ð23Þ
s2ð∞; ζ; τÞ ¼ 0 ð24Þ
The dimensionless boundary condition along the screen is expressed as
s2ð1; ζ; τÞ ¼ 1; ζl< ζ < 1; τ > 0 ð25Þ
In dimensionless form, continuity conditions become
s1ð1; ζ; τÞ ¼ s2ð1; ζ; τÞ; 0< ζ < ζl; τ > 0 ð26Þ and ∂s 1ðρ; ζ; τÞ ∂ρ ρ¼1¼ ∂s2ðρ; ζ; τÞ ∂ρ ρ¼1; 0< ζ < ζl; τ > 0 ð27Þ The solution for the dimensionless drawdown solutions can be obtained by taking Laplace transforms of governing equations Eqs. (17) and (18) using the initial condition (19), and the results are
Fig. 9. Effect ofα on drawdown in 100-m-thick aquifer when σ ¼ 103,κ ¼ 1 at ζ ¼ 0 for α ¼ 100, 101, and 102
Fig. 10. Error of estimated dimensionless drawdown along screen cal-culated from Eqs. (13) and (14) for different number of terms of infinite series and Shanks method
∂2~s 1 ∂ρ2 þ 1 ρ ∂~s 1 ∂ρ þ αw∂ 2~s 1 ∂ζ2 ¼ p~s 1; 0≤ ρ ≤ 1; 0≤ ζ ≤ ζl ð28Þ and ∂2~s 2 ∂ρ2 þ 1 ρ∂~s 2 ∂ρ þ αw∂ 2~s 2 ∂ζ2 ¼ p~s2; 1≤ ρ < ∞; 0≤ ζ ≤ 1 ð29Þ The transformed boundary conditions at the bottom and top of the aquifer for Regions 1 and 2 can be written as
∂~s 1ðρ; ζ; pÞ ∂ζ ζ¼0¼ ∂~s2ðρ; ζ; pÞ ∂ζ ζ¼0¼ 0 ð30Þ ∂~s 1ðρ; ζ; pÞ ∂ζ ζ¼ζl ¼ 0; 0< ρ < 1 ð31Þ and ∂~s 2ðρ; ζ; pÞ ∂ζ ζ¼1¼ σ αw · p ·~s2ðρ; ζ ¼ 1; pÞ ð32Þ
Likewise, the transformed boundary conditions atρ ¼ 0 and ∞ are ∂~s 1ð0; ζ; pÞ ∂ρ ¼ 0; 0< ζ < ζl ð33Þ and ~s2ð∞; ζ; pÞ ¼ 0 ð34Þ
After taking the Laplace transform, the boundary condition along the well screen is
~s2ð1; ζ; pÞ ¼
1
p; ζl< ζ < 1; τ > 0 ð35Þ and continuity conditions become
~s2ð1; ζ; pÞ ¼ ~s1ð1; ζ; pÞ; 0< ζ < ζl ð36Þ and ∂~s 1ðρ; ζ; pÞ ∂ρ ρ¼1 ¼ ∂~s2ðρ; ζ; pÞ ∂ρ ρ¼1; 0< ζ < ζl ð37Þ
Assume that~s1and~s2are the product of two distinct functions,
i.e., ~s1ðρ; ζ; pÞ ¼ F1ðρ; pÞG1ðζ; pÞ and ~s2ðρ; ζ; pÞ ¼ F2ðρ; pÞ
G2ðζ; pÞ, respectively. Eqs. (12) and (13) can be respectively
transformed as G1∂ 2F 1 ∂ρ2 þ G1 1 ρ∂F 1 ∂ρ þ αwF1∂ 2G 1 ∂ζ2 ¼ pF1G1 ð38Þ and G2∂ 2F 2 ∂ρ2 þ G2 1 ρ∂F∂ρ2þ αwF2∂ 2G 2 ∂ζ2 ¼ pF2G2 ð39Þ
Dividing Eqs. (38) and (39) by F1G1 and F2G2, respectively,
Eqs. (38) and (39) can then be separated into the following two systems of ordinary differential equations after some arrangements:
∂2G 1 ∂ζ2 þ ω 1m αw G1¼ 0 ð40Þ ∂2F 1 ∂ρ2 þ 1 ρ ∂F1 ∂ρ ½p þ ω1mF1¼ 0 ð41Þ and ∂2G 2 ∂ζ2 þ ω 2n αw G2¼ 0 ð42Þ ∂2F 2 ∂ρ2 þ 1 ρ∂F∂ρ2 ½p þ ω2nF2¼ 0 ð43Þ
where ω1m and ω1n = separation constants.
The solutions of Eqs. (40) and (42) subject to the boundary in Eq. (30) are, respectively, as follows:
G1ðζ; pÞ ¼ a1mðpÞ cosðΩ1mζÞ ð44Þ and G2ðζ; pÞ ¼ a2nðpÞ cosðΩ2nζÞ ð45Þ whereΩ1m¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ω1m=αw p ;Ω2n¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ω2n=αw p ; and a1mðpÞ and a2nðpÞ
= constants with respect toζ. In addition, substituting Eq. (44) into Eq. (31) yields the following equation:
sinðΩ1mζlÞ ¼ 0 ð46Þ
The eigenvaluesΩ1min Eq. (44) can then determined by solving
Eq. (46) and results in Ω1m¼
mπ ζl
; m¼ 0; 1; 2; … ð47Þ
Similarly, substituting Eq. (45) into Eq. (32) gives the following equation:
Ω2nsinðΩ2nÞ ¼ σα w
p cosðΩ2nÞ; n¼ 0; 1; 2; … ð48Þ
Eq. (48) can be solved to obtain the eigenvaluesΩ2nin Eq. (45).
The solutions of Eqs. (41) and (43) are, respectively, as follows: F1ðρ; pÞ ¼ c1mðpÞI0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pþ ω1m p ρÞ ð49Þ and F2ðρ; pÞ ¼ d2nðpÞK0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pþ ω2n p ρÞ ð50Þ
where c1mðpÞ and d2nðpÞ = constants. Note that d1mðpÞ and c2nðpÞ
equal zero when using the boundary conditions of Eqs. (33) and (34), respectively.
The product of Eqs. (44) and (49) gives the general solution of Eq. (36) as ~s1mðρ; ζ; pÞ ¼ A1mðpÞI0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pþ ω1m p ρÞ cosðΩ1mζÞ m¼ 0; 1; 2; … ð51Þ
where A1mðpÞ = product of a1mðpÞ and c1mðpÞ. On the other hand,
the product of Eqs. (35) and (50) forms the general solution of Eq. (39) as ~s2nðρ; ζ; pÞ ¼ A2nðpÞK0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pþ ω2n p ρÞ cosðΩ2nζÞ n¼ 0; 1; 2; … ð52Þ where A2nðpÞ = product of a2nðpÞ and d2nðpÞ. Accordingly, the
lin-ear combination of all solutions of m yields the complete solution for~s1ðρ; ζ; pÞ ~s1ðρ; ζ; pÞ ¼ X∞ m¼0 A1mðpÞI0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pþ ω1m p ρÞ cosðΩ1mζÞ ð53Þ
Similarly, the complete solution for~s2ðρ; ζ; pÞ can be obtained as ~s2ðρ; ζ; pÞ ¼ X∞ n¼0 A2nðpÞK0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pþ ω2n p ρÞ cosðΩ2nζÞ ð54Þ
The coefficients A1mðpÞ and A2nðpÞ are unknowns at this stage
and can be solved from the following equation obtained by substi-tuting Eqs. (53) and (54) into Eqs. (35) and (36), respectively, as
X∞ n¼0 A2nðpÞK0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pþ ω2n p Þ cosðΩ2nζÞ ¼ 1 p; ζl< ζ < 1 ð55Þ and X∞ m¼0 A1mðpÞI0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pþ ω1m p Þ cosðΩ1mζÞ ¼X∞ n¼0 A2nðpÞK0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pþ ω2n p Þ cosðΩ2nζÞ; 0< ζ < ζl ð56Þ
Eqs. (55) and (56) are organized and expressed as X∞ n¼0 A2nðpÞK0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pþ ω2n p Þ cosðΩ2nζÞ ¼X∞ m¼0 A1mðpÞI0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pþ ω1m p Þ cosðΩ1mζÞ; 0< ζ < ζl ¼1 p; ζl< ζ < 1 ð57Þ
To obtain concise solutions, A02nðpÞ ¼ A2nðpÞK0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω2n
p
Þ and A01mðpÞ ¼ A1mðpÞI0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω1m
p
Þ are further defined, and Eq. (57) can be rewritten as X∞ n¼0 A02nðpÞ cosðΩ2nζÞ ¼ f ðζÞ; 0< ζ < 1 ð58Þ where fðζÞ ¼X ∞ m¼0 A01mðpÞ cosðΩ1mζÞ; 0< ζ < ζl ¼1 p; ζ1< ζ < 1 ð59Þ
The term on the left-hand side (LHS) of Eq. (58) is a half-range Fourier cosine series of the function on the right-hand side (RHS) for the region 0< ζ < 1. The coefficient A02nðpÞ can then be ob-tained from the properties of the Fourier series as
A02nðpÞ ¼ R1 0RcosðΩ2nζÞf ðζÞdζ 1 0cos2ðΩ2nζÞdζ ð60Þ Carrying out the integration in Eq. (60) and simplifying the re-sult yields the coefficient A02nðpÞ as expressed in Eq. (16).
Similarly, substituting Eqs. (53) and (54) into Eq. (37), one can obtain X∞ m¼0 A01mðpÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω1m p I1ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω1m p Þ I0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω1m p ÞcosðΩ1mζÞ ¼ X∞ n¼0 A02nðpÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω2n p K1ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω2n p Þ K0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω2n p ÞcosðΩ2nζÞ; 0< ζ < ζl ð61Þ
From Eq. (61), the coefficient A01mðpÞ can be determined as Eq. (15).
Accordingly, based on the coefficients A01mðpÞ and A02nðpÞ, the complete solution for~s1 and~s2can be obtained as Eqs. (13) and
(14), respectively.
Notation
The following symbols are used in this paper: A01m = A1mðpÞI0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω1m p Þ; A02n = A2nðpÞK0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω2n p Þ; s1 = s1=sw; s2 = s2=sw;
I0ð·Þ = modified Bessel function of first kind of order 0;
I1ð·Þ = modified Bessel function of first kind of order 1;
i0m = I0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω1m p Þ=½ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω1m p · I1ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω1m p Þ
K0ð·Þ = modified Bessel function of second kind of order 0;
K1ð·Þ = modified Bessel function of second kind of order 1;
k0n =½ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω2n p · K1ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω2n p Þ=K0ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffipþ ω2n p Þ; β = b=rw; ζ = z=b; ζl= zl=b; κ = Kz=Kr; Λmn =fsin½ðΩ1mþ Ω2nÞζl=ðΩ1mþ Ω2nÞgþ fsin½ðΩ1m Ω2nÞζl=ðΩ1m Ω2nÞg; ρ = r=rw; τ = Krt=Ssr2w; Ω1m =β ffiffiffiffiffiffiffiffiffiffiffiffi ω1m=k p ; and Ω2n =β ffiffiffiffiffiffiffiffiffiffiffiffi ω2n=k p . Acknowledgments
Research leading to this work has been partially supported by the grants from Taiwan National Science Council under the contract numbers NSC 96-2221-E-009-087-MY3, NSC 98-3114-E-007-015, and NSC 99-NU-E-009-001.
References
Abdul, A. S. (1992).“A new pumping strategy for petroleum product re-covery from contaminated hydrogeologic systems: Laboratory and field evaluations.”Ground Water Monit. Rem., 12(1), 105–114.
Boulton, N. S. (1955).“Unsteady radial flow to a pumped well allowing for delayed yield from storage.” Gen. Assera. Rome, Tome II, Int. Assoc. Sci. Hydrol. Publ., 37, 472–477.
Chang, Y. C., and Yeh, H. D. (2009).“New solutions to the constant-head test performed at a partially penetrating well.”J. Hydrol., 369(1–2), 90–97.
Chang, Y. C., and Yeh, H. D. (2010).“A new analytical solution solved by triple series equations method for constant-head tests in confined aqui-fers.”Adv. Water Resour., 33(6), 640–651.
Chang, Y. C., Yeh, H. D., and Chen, G. Y. (2010).“Transient solution for radial two-zone flow in unconfined aquifers under constant-head test.”
Hydrol. Processes, 24(11), 1496–1503.
Chen, C. S., and Chang, C. C. (2003).“Well hydraulics theory and data analysis of the constant head test in an unconfined aquifer with the skin effect.”Water Resour. Res., 39(5), 1121–1135.
Hantush, M. S. (1964).“Hydraulics of wells.” Advances in hydroscience, Vol. 1, V. T. Chow, ed., Academic, San Diego, 281–432.
Hiller, C. K., and Levy, B. S. (1994).“Estimation of aquifer diffusivity from analysis of constant-head pumping test data.”Ground Water, 32(1), 47–52.
Javandel, I., and Zaghi, N. (1975).“Analysis of flow to an extended fully penetrating well.”Water Resour. Res., 11(1), 159–164.
Jones, L. (1993).“A comparison of pumping and slug tests for estimating the hydraulic conductivity of unweathered Wisconsin age till in Iowa.”
Ground Water, 31(6), 896–904.
Jones, L., Lemar, T., and Tsai, C. T. (1992).“Results of two pumping tests in Wisconsin age weathered till in Iowa.”Ground Water, 30(4), 529–538.
Kirkham, D. (1959). “Exact theory of flow into a partially penetrating well.”J. Geophys. Res., 64(9), 1317–1327.
Mishra, S., and Guyonnet, D. (1992). “Analysis of observation-well response during constant-head testing.” Ground Water, 30(4), 523–528.
Moench, A. F. (1997). “Flow to a well of finite diameter in a homo-geneous, anisotropic water table aquifer.”Water Resour. Res., 33(6), 1397–1407.
Murdoch, L. D., and Franco, J. (1994).“The analysis of constant drawdown wells using instantaneous source functions.”Water Resour. Res., 30(1), 117–127.
Neuman, S. P. (1972).“Theory of flow in unconfined aquifers considering delayed response of the water table.” Water Resour. Res., 8(4), 1031–1045.
Neuman, S. P. (1974).“Effects of partial penetration on flow in unconfined aquifers considering delayed aquifer response.”Water Resour. Res., 10(2), 303–312.
Peng, H. Y., Yeh, H. D., and Yang, S. Y. (2002).“Improved numerical evaluation of the radial groundwater flow equation.” Adv. Water Resour., 25(6), 663–675.
Renard, P. (2005). “Approximate discharge for constant head test with recharging boundary.”Ground Water, 43(3), 439–442.
Singh, S. K. (2007).“Simple approximation of well function for constant drawdown variable discharge artesian wells.” J. Irrig. Drain. Eng., 133(3), 282–285.
Stehfest, H. (1970). “Numerical inversion of Laplace transforms.”
Commun. ACM, 13(1), 47–49.
Tartakovsky, G. D., and Neuman, S. P. (2007).“Three-dimensional satu-rated-unsaturated flow with axial symmetry to a partially penetrating well in a compressible unconfined aquifer.”Water Resour. Res., 43, W01410, 10.1029/2006WR005153.
Uraiet, A. A., and Raghavan, R. (1980).“Unsteady flow to a well producing at a constant pressure.” J. Pet. Technol., 32(10), 1803–1812. Wang, C. T., and Yeh, H. D. (2008).“Obtaining the steady-state drawdown
solutions of constant-head and constant-flux tests.”Hydrol. Processes, 22(17), 3456–3461, 10.1002/hyp.6950.
Yang, S. Y., and Yeh, H. D. (2002).“Solution for flow rates across the wellbore in a two-zone confined aquifer.”J. Hydraul. Eng., 128(2), 175–183.
Yang, S. Y., and Yeh, H. D. (2005).“Laplace-domain solutions for radial two-zone flow equations under the conditions of constant-head and partially penetrating well.”J. Hydraul. Eng., 131(3), 209–216. Yeh, H. D., and Yang, S. Y. (2006). “A novel analytical solution for
constant-head test in a patchy aquifer.”Int. J. Numer. Anal. Methods Geomech., 30(12), 1213–1230.
