Non-Darcian flow to a well in an aquifer–aquitard system

Advances in Water Resources 31 (2008) 1754–1763

Contents lists available at ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

Non-Darcian flow to a well in an aquifer–aquitard system Zhang Wen a,b, Guanhua Huang a,b,*, Hongbin Zhan c a b c

Department of Irrigation and Drainage, College of Water Conservancy and Civil Engineering, China Agricultural University, Beijing 100083, PR China Chinese–Israeli International Center for Research and Training in Agriculture, China Agricultural University, Beijing 100083, PR China Department of Geology and Geophysics, Texas A&M University, College Station, TX 77843-3115, USA

a r t i c l e

i n f o

Article history: Received 12 April 2008 Received in revised form 31 August 2008 Accepted 18 September 2008 Available online 26 September 2008 Keywords: Non-Darcian flow Power law Finite difference method Laplace transform Linearization method

a b s t r a c t In this study, we use a linearization procedure and a finite difference method to solve non-Darcian flow to a well in an aquifer–aquitard system. The leakage effect is considered. Flow in the aquifer is assumed to be non-Darcian and horizontal, whereas flow in the aquitard is assumed to be Darcian and vertical. The Izbash equation [Izbash SV. O filtracii V Kropnozernstom Materiale. USSR: Leningrad; 1931 [in Russian]] is employed to describe the non-Darcian flow. The wellbore storage is also considered in this study. An approximate semi-analytical solution has been obtained by the linearization procedure, and a numerical solution has been obtained by using a finite difference method. The previous solutions for Darcian flow case and non-Darcian flow case without leakage can be described as special cases of the new solutions. The error caused by the linearization procedure has also been analyzed. The relative error caused by the linearization procedure is nearly 100% at early times, and decreases to zero at late times. We have also compared the results in this study with Wen et al. [Wen Z, Huang G, Zhan H. A numerical solution for non-Darcian flow to a well in a confined aquifer using the power law function. J Hydrol, 2008d [in revision]] in which the leakage effect is not considered, and Hantush and Jacob [Hantush MS, Jacob CE. Non-steady radial flow in an infinite leaky aquifer. Trans Am Geophys Union 1955;36(1):95–100] who investigated a similar problem in Darcian flow case. The comparison of this study and Wen et al. (2008d) indicates the dimensionless drawdown in the aquifer with leakage is less than that without leakage, and the leakage has little effect at early times. The comparison between the results of this study and that of Hantush and Jacob (1955) indicates that the dimensionless drawdown in the aquifer for non-Darcian flow is larger at early times and smaller at late times, than their counterparts for Darcian flow. A larger dimensionless non-Darcian conductivity kD results in a smaller dimensionless drawdown in the aquifer at late times, and leads to a larger dimensionless drawdown in the aquifer at early times. A smaller dimensionless leakage parameter BD results in a smaller drawdown at late times, and the leakage does not affect the early-time drawdown. The analysis of the dimensionless drawdown inside the well has also been included in this study when the wellbore storage is considered. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction Groundwater hydraulics in leaky aquifers has been investigated for decades (e.g. [13–15,24,29,41,57]). The first mathematical model for flow to a well in leaky aquifers was developed by Hantush and Jacob [13]. In that study, they presented an analytical solution for transient flow to a well in leaky aquifers based on a series of simplified assumptions, e.g. flow in aquitard is vertical, pumping rate of well is constant, storage of aquitard is negligible, etc. This solution has been extended by Hantush [14] to consider flow to a constant-head pumping well, and by Hantush [15] to

* Corresponding author. Address: Chinese–Israeli International Center for Research and Training in Agriculture, China Agricultural University, Beijing 100083, PR China. Tel.: +86 10 62737144; fax: +86 10 62737138. E-mail addresses: [email protected] (Z. Wen), [email protected] (G. Huang), [email protected] (H. Zhan). 0309-1708/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2008.09.002

consider the storage of aquitard. Recently, many researches about the hydraulics in leaky aquifers have been done (e.g. [1,16,22,23,30,54–56]). For instance, Li and Neuman [23] investigated flow to a well in a five-layer system, and validated their solution with a pumping test in the Oxnard Basin, California, USA. Hunt and Massmann [16] presented an analytical solution for vapor flow to a trench in a leaky aquifer, and the solution can be used to evaluate the permeability and the leakage parameter. Yeh et al. [54] used global optimization methods to determine the parameters for the leaky aquifers. Zhan and Bian [55] have provided some practically useful solutions for estimating the rate and volume of water leaked through a semi-confining layer within any given radial distance from the pumping well. All of these studies are based on the assumption that flow is Darcian. However, under some circumstances, flow can be non-Darcian because of the high velocities near a pumping well [2,37]. For instance, if the pumping rate of the well is sufficiently large, which will result in a relatively

Z. Wen et al. / Advances in Water Resources 31 (2008) 1754–1763

1755

Nomenclature r rw rc t q(r, t) s(r, t) sw(t) m S n k m1 k1 B

distance from the center of well (L) effective radius of well screen (L) radius of well casing (L) pumping time (T) specific discharge (L/T) drawdown (L) drawdown inside well (L) aquifer thickness (L) storage coefficient of aquifer power index, an empirical constant in the Izbash equation quasi hydraulic conductivity, an empirical constant in the Izbash equation ((L/T)n) thickness of aquitard (L) hydraulic conductivity of aquitard (L/T) leakage parameter, defined as m  m1/k1 (L/T)

large velocity near the pumping well, non-Darcian flow is likely to occur. Non-Darcian flow has caught the scientists’ attention for more than one hundred years. There are two general types of non-Darcian flow: e.g. pre- and post-linear flow. In hydraulic engineering, post-linear non-Darcian flow often occurs (e.g. [18,20,21,40]). There are two commonly used formulae to describe the relationship between the hydraulic gradient and the specific discharge for non-Darcian flow. The first one is the Forchheimer equation [12], which states that the hydraulic gradient is a second-order polynomial function of the specific discharge. The second one is the Izbash equation [17], which states that the hydraulic gradient is a power function of the specific discharge. Both these two equations describe the non-Darcian flow equivalently well [3,53]. In addition to these two equations, the Brinkman equation [4] can also be used to describe the non-Darcian flow. It stated that an extra viscous term should be added to the classical Darcy equation. The Brinkman equation is semi-empirical in nature, and it has been validated by a detailed numerical solution of the Stokes’ equations [25]. Non-Darcian flow has been incorporated into flow to a pumping well in confined aquifers since last century (e.g. [32–36,46,51,52]). Both analytical and numerical methods have been used to solve the non-Darcian flow problems. A similarity method named the Boltzmann transform has been used to solve such non-Darcian flow problems by Sen [32–36]. Wen et al. [46,49] have also used the Boltzmann transform to solve the non-Darcian flow problems in a fractured and a porous medium, respectively. However, the Boltzmann transform has been criticized by some scientists (e.g. [5,6,26,49]) because of its non-rigorous mathematical derivation. Ramey [31] was probably the first one to discuss such an issue in detail. He found that the Boltzmann approach only provided a late-time approximation drawdown in this configuration, and it could not be used for generation of the type curves. Recently, Wen et al. [47,48] developed a linearization procedure to solve non-Darcian flow to a well in a confined aquifer, and obtained some approximate analytical solutions. The linearization solution is shown to work very well at late times, whereas it will underestimate the drawdown at early times [26,50]. Up to now, some numerical solutions for non-Darcian flow have also been obtained based on the finite difference method or the finite element method (e.g. [7,10,11,19,26]). All these studies are based on a presumption that the non-Darcian flow can be

p Q

C(x) Km(x) rD rwD rcD tD qD sD swD kD BD

a

Laplace variable (1/T) pumping discharge (L3/T) Gamma function second kind m-order modified Bessel function dimensionless distance defined in Table 1 dimensionless radius of well screen defined in Table 1 dimensionless radius of well casing defined in Table 1 dimensionless time defined in Table 1 dimensionless specific discharge defined in Table 1 dimensionless drawdown defined in Table 1 dimensionless drawdown inside the well defined in Table 1 dimensionless non-Darcian flow hydraulic conductivity defined in Table 1 dimensionless leakage parameter defined in Table 1 a parameter used in Eq. (16)

described by the Forchheimer equation. For instance, Mathias et al. [26] developed a finite difference solution for non-Darcian flow to a well in a confined aquifer based on the Forchheimer equation. They found that both the Boltzmann solution [34] and the linearization solution [47] worked very well at late times, and some differences had been found at early and moderate times. Meanwhile, Wen et al. [50] have also obtained a numerical solution for non-Darcian flow in porous media on the basis of the Izbash equation. To our knowledge, there are few studies about non-Darcian flow to a well considering the leakage effect. Sen [37] investigated nonDarcian flow to a well in leaky aquifers with the Izbash equation by using a volumetric approach. This method was extended to solve a similar problem with the Forchheimer equation by Birpinar and Sen [2]. Both of these two studies are based on their previous solutions of Sen [35] and Sen [34] which were obtained by using the Boltzmann transform. As mentioned before, the Boltzmann transform cannot be used to solve the non-Darcian flow problems in a rigorous mathematical sense. Thus, the solutions of Sen [37] and Birpinar and Sen [2] for non-Darcian flow in leaky aquifers are questionable. In this paper, we will use the linearization procedure to solve the non-Darcian flow problem in leaky aquifers. The Izbash equation will be employed to describe the non-Darcian flow. A finite difference solution will also be obtained to verify the linearization solution.

2. Problem statement 2.1. Governing equation The aquifer–aquitard system investigated is schematically shown in Fig. 1. Similar to the assumptions used in Hantush and Jacob [13], the aquifer–aquitard system is homogenous and laterally infinite; the storage capacity of the aquitard is ignored; flow in the aquifer is horizontal, and assumed to be non-Darcian in the entire aquifer, whereas flow in the aquitard is vertical and assumed to be Darcian; the pumping well fully penetrates the aquifer and has a constant pumping rate; the system is hydrostatic before the pumping starts. Based on these assumptions, the problem discussed here can be described as follows when the wellbore storage is excluded for an infinitesimal-diameter well:

1756

Z. Wen et al. / Advances in Water Resources 31 (2008) 1754–1763 Table 1 Dimensionless variables used in this study

Land Surface

Land Surface

s

sw

Water-table Aquifer

rD ¼

Water-table Aquifer

r m

qD ¼ 

4pm2 qðr; tÞ Q

2rc

Aquitard

Aquitard

rcD ¼

rc m

sD ¼

rwD ¼

rw m

swD ¼

1=n

4 pk m sðr; tÞ Q

2rw

Confined Aquifer

r

1=n

1=n

BD ¼

Fig. 1. The schematic diagram of the system.

1 o½rqðr; tÞ sðr; tÞ S osðr; tÞ  ¼ ; r or B m ot sðr; 0Þ ¼ 0;

ð2Þ

sð1; tÞ ¼ 0;

ð3Þ

lim 2pmrqðr; tÞ ¼ Q;

ð4Þ

r!0

1=n

Bk m2 

ð1Þ kD ¼

tD ¼

1=n m2 k

1=n

where r is the distance from the pumping well; t is the pumping time; q(r, t) and s(r, t) are the specific discharge and drawdown at distance r and time t, respectively; S is the storage coefficient of aquifer; m is the thickness of aquifer; Q is the pumping rate, which is positive; B is the leaky parameter which is defined as m  m1/k1, where m1 and k1 are the thickness and the hydraulic conductivity of aquitard, respectively. Using Eq. (3) implies that the flow in the whole aquifer is assumed to be non-Darcian, and can be described by the Izbash equation. Actually, the flow velocity decreases as the distance increases. Non-Darcian flow might no longer hold when the distance is far enough. Further discussion about this issue can be found in the following discussion part. If the wellbore storage is included, the boundary Eq. (4) should be replaced by

2pr w mqðr; tÞjr!rw  p

r2c

dsw ðtÞ ¼ Q ; dt

ð5Þ

in which rw is the radius of well screen, and rc is the radius of well casing; sw(t) is the drawdown inside the well. In this study, we employ the Izbash equation to describe non-Darcian flow:

osðr; tÞ ½qðr; tÞ ¼ k ; or n

1 r2 dswD ðt D Þ ¼ 1: ðrD qD ÞjrD !rwD þ cD 2 dt D 4S

ðqD Þn ¼ kD

osD ; or D

o½r D qD  sD osD  ¼ ; r D orD BD ot D sD ðrD ; 0Þ ¼ 0; 

sD ð1; t D Þ ¼ 0; 1 lim r D qD ¼ 1: r D !0 2

ð7Þ ð8Þ ð9Þ

ð12Þ

where kD can be regarded as a dimensionless non-Darcian hydraulic conductivity. 3. Approximate analytical solution with linearization procedure 3.1. Solution without wellbore storage Substituting the Izbash equation Eq. (12) into the governing equation Eq. (7) results in 1=n

kD o2 sD kD osD þ n or 2D rD or D

! 

osD or D

1n n 

sD osD ¼ : BD otD

ð13Þ

Eq. (13) can be expressed in an alternatively way as

 n1  n1 o2 sD n osD sD n osD n n osD n osD þ   ¼  : 1=n or D orD ot D or 2D r D or D BD k1=n kD D The term

Defining the dimensionless variables in Table 1, then Eqs. (1)– (4) can be transformed as

ð11Þ

In terms of the dimensionless definition of the variables, the Izbash equation can be changed to

ð6Þ

2.2. Dimensionless transformation

1=n

1k . The boundary less leakage parameter defined as BD ¼ Bkm2 ¼ mmk 1 Eq. (5) with the wellbore storage can also been transformed as:

1=n

where n and k are empirical constants of the aquifer. While for postlinear flow, n is larger than 1 and less than 2. If n equals to 1, flow is Darcian. In this case, k becomes the hydraulic conductivity of the aquifer. Thus, k can be regarded as the quasi non-Darcian hydraulic conductivity of the aquifer, which reflects how easy the aquifer can transmit water.

k t Sm

1n

Q 4p

4pk m sw ðtÞ Q

ð14Þ

 n1 n osD  or in Eq. (14) is non-linear, which cannot be D

solved analytically. Fortunately, Wen et al. [47,48] have developed a linearization method to solve such a problem, which has been proven that it works very well at late times by Mathias et al. [26]. Similar to what have been done in Wen et al. [47,48], the  n1 n osD can be approximated as follows by comnon-linear term  or D bining Eq. (10) with Eq. (12)



osD qnD 2n ¼  n : or D kD r D kD

ð15Þ

With this approximation, Eq. (14) can be changed to

ð10Þ 2

Notice that a negative sign has used for defining qD in Table 1. The subscript ‘‘D” means the dimensionless variables. BD is a dimension-

o sD n osD sD 1n osD þ  ar  ar 1n ; D ot D or 2D r D or D BD D

ð16Þ

1757

Z. Wen et al. / Advances in Water Resources 31 (2008) 1754–1763

in which a ¼ 21nn k . With the Laplace transform, Eq. (16) becomes: D

  2 d sD n dsD sD a þ  n1 þ ap ¼ 0; 2 BD dr D rD dr D r D

ð17Þ

where sD is the drawdown in Laplace domain, p is the Laplace variable in respect to tD, and the over bar means the variables in the Laplace domain. It should be pointed out that the initial condition Eq. (8) has been used in the above transform. The boundary condition can be expressed as follows with the Laplace transform:

sD ð1; pÞ ¼ 0; lim rnD

rD !0

ð18Þ

osD 2n ¼ : or D kD p

ð19Þ

With considering the boundary conditions and some properties of the Bessel functions, the solution for the dimensionless drawdown without wellbore storage in the Laplace domain can be expressed as



2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi3n 1 a þ ap rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2nþ1 3n BD 1n 2 3n a qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r D2 K 1n sD ¼ þ ap ; r D2  2  3n a 3  n B D C 3n kD p BD þ ap

ð20Þ

where Kv(x) is the second kind modified Bessel function with order t, C(x) is the Gamma function. It seems that it is difficult, if not impossible, to invert Eq. (20) analytically. There are many developed methods to do this numerical Laplace inverse transform such as Stehfest [43,44], Crump [8], Talbot [45], and de Hoog et al. [9], among others. The de Hoog method has been commonly used because it is accurate for a wide range of functions and it also performs reasonably well in the neighborhood of discontinuity [27]. For the problem proposed here, both the Stehfest method we used before and the de Hoog method work very well. In this study, we employ the de Hoog method to do the numerical inversion. Notice that when BD goes to infinity, which means the leakage does not exist, the solution of Eq. (20) can be changed to

sD ¼

2  1 pffiffiffiffiffiffi3n   1n 2nþ1 3n ap pffiffiffiffiffiffi 2 3n  2  r D2 ap : 1n pffiffiffiffiffiffi r D2 K 3n 3n C 3n kD p ap

ð21Þ

Eq. (21) is the same as the solution of Wen et al. [50] in which the leakage is not considered. If the power index n is equal to unity, which means flow in the aquifer is Darcian, notice that kD and a are equal to unity as well, then Eq. (20) can be simplified to

kD

r n os r 2 dswD ðtD Þ D D jrD !rwD  cD  1: dt D 2 or D 4S

ð25Þ

With the Laplace transform, the solution with wellbore storage in the Laplace domain can be expressed as   1n 3n qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 rD2 K 1n 3n rD2 BaD þ ap 3n     : sD ¼ ffi 3n qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1n 3n qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2  n 1n qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 a þ ap þ rcD pr 2 K 2 a þ ap 2 2 rwD BaD þ apK 2 3n rwD p kD rwD wD 1n 3n r wD 2 BD 4S BD 3n

3n

ð26Þ Eq. (26) can also be evaluated by the numerical inversion of the Laplace transform. Similarly, we used the de Hoog method to do this calculation. 3.3. Approximate analytical solution at steady-state D ¼ 0. Then the governing equation Eq. At steady-state, one has os ot D (16) can be changed to:

o2 sD n osD sD 1n þ  ar ¼ 0; or 2D rD or D BD D

ð27Þ

The solution for Eq. (27) with considering the boundary conditions can be expressed as

2nþ1 sD ¼

C





1 3n

2 3n



2 qffiffiffiffi3n

kD

a

BD

qffiffiffiffi a



1n

rD2 K 1n

BD

3n

2 3n r2 3n D

rffiffiffiffiffiffi a : BD

ð28Þ

If BD goes to infinity, which means the leakage can be ig effect m nored, applying the approximation K m ðxÞ  Cð2mÞ 2x ; m > 0 when x goes to zero [42], Eq. (28) can be simplified as

sD 

r D1n 2n : kD ðn  1Þ

ð29Þ

Eq. (29) is the same as that has been obtained by Wen et al. [47] in confined aquifers. 4. Numerical solution

Eq. (23) is the same as the solution of Hantush and Jacob [13] for Darcian flow to a well in leaky aquifers.

In this study, we used the finite difference method to do the numerical simulation. Here we will consider the wellbore storage, because the case without wellbore storage can be easily obtained by giving a small value for the radius of the well. Because the drawdown changes greatly near the pumping well and changes progressively less when moving away from the well, it is better to have finer grid spaces near the well and to have progressively coarser grid spaces when moving away from the well. For the case of convenience, we discretize the dimensionless space rD logarithmically [26,51]. The dimensionless space domain [rwD, reD] was discretized into N nodes excluding the two boundary nodes rwD and reD, where reD is a relatively large dimensionless distance which is used to approximate the infinite boundary. For any node of ri, rwD < ri < reD, i = 1, 2, . . . , N, one has

3.2. Solution with wellbore storage

ri ¼ ðr i1=2 þ riþ1=2 Þ=2;

The boundary condition Eq. (11) with wellbore storage can be expressed alternatively as

where ri+1/2 is calculated as follows:

sD ¼

2 K0 p

sffiffiffiffiffiffiffiffiffiffiffiffiffiffi! 1 rD þp : BD

ð22Þ

Using the inverse Laplace transform, the solution of the dimensionless drawdown in real time domain can be expressed as

sD ðr D ; t D Þ ¼

Z

1

r2 D 4tD

  1 r2 exp y  D dy: y 4BD y

 1=n 1 osD r2 dswD ðt D Þ 1=n lim þ cD ¼ 1: r D kD  r D !r wD 2 dt D orD 4S

ð23Þ

ð24Þ

Eq. (24) is also a non-linear boundary condition. Similar to the linearization procedure used in the governing equation, Eq. (24) can be approximated as

i ¼ 1; 2; . . . ; N;

ð30Þ

 log10 ðr eD Þ  log10 ðr wD Þ log10 ðr iþ1=2 Þ ¼ log10 ðrwD Þ þ i ; N i ¼ 0; 1; . . . ; N:

ð31Þ

After the spatial discretization, the governing equation Eq. (7) can be reduced to the following differential equation in respect to the dimensionless time tD

1758

Z. Wen et al. / Advances in Water Resources 31 (2008) 1754–1763

r i1=2 qi1=2  r iþ1=2 qiþ1=2 si dsi   ; dt D ri ðriþ1=2  ri1=2 Þ BD

i ¼ 1; 2; . . . ; N;

2

10

ð32Þ

10

in which qi and si are the dimensionless specific discharge qD and dimensionless drawdown sD at node i, respectively. In terms of the Izbash equation Eq. (12), one can obtain:

si1  si r i  r i1

D

1=n ;

i ¼ 2; 3; . . . ; N;

−1

ð33Þ

10

r =10

s

qi1=2  kD

r =1

0

10

D

1=n



Hantush and Jacob (1955) Numerical solution Linearization solution

1

and 1=n



qiþ1=2  kD

si  siþ1 r iþ1  r i

1=n ;

i ¼ 1; 2; . . . ; N  1:

ð34Þ

At the boundaries, one has



n=1 k =1 D B =10

−4

ð35Þ

D

−5

10

−2

10

ð36Þ

in which swD is the dimensionless drawdown inside the well, which can be approximated as following with considering Eq. (11):

  dswD 4S 1  2 1  r wD q11=2 : 2 dtD r cD

−3

10 10

1=n

swD  s1 ; r1  rwD  1=n sN 1=n  kD ; reD  r N 1=n

q11=2  kD qNþ1=2

D

−2

10

ð37Þ

With these preparations, the problem discussed here can be solved numerically by using the stiff integrator ODE15s, which can be found in any standard version of MATLAB [26,38,39]. We have developed a MATLAB program named NFL (Non-Darcian Flow in Leaky aquifers) to do this calculation. It is found that sufficiently accurate solutions can be obtained by setting reD = 108 and N = 3000. For the case without wellbore storage, the dimensionless drawdown is found to be insensitive to the well radius when the dimensionless radius of well casing rcD and the dimensionless radius of well screen rwD are less than 104. 5. Results and discussion 5.1. Comparison with the solution of Hantush and Jacob (1955) for Darcian flow case We have compared our solutions obtained by the linearization procedure and numerical method with the solution of Hantush and Jacob [13] for Darcian flow case by setting the parameter n = 1, as shown in Fig. 2. It is evident to see that our linearization solution and numerical solution agree very well with the solution of Hantush and Jacob [13] for Darcian flow case. Actually, the linearization procedure is exactly right for the Darcian flow case. Therefore, it is no surprise that the linearization solution agrees well with the solution of Hantush and Jacob [13]. This also implies that the linearization method might work generally well when the non-Darcian effect such as flow turbulence is not very strong, meaning that n is very close to 1. The general agreement between the numerical solution and that of Hantush and Jacob [13] indicates that the numerical error is negligible. In the following, we will use the numerical solution to verify the linearization solution. 5.2. Dimensionless drawdown-time behavior without wellbore storage Fig. 3 is about the comparison of the numerical solution and the linearization solution with kD = 10, BD = 10 for rD = 1, and rD = 10, respectively: (a) n = 1.1; (b) n = 1.5; (c) n = 2.0 (fully turbulent flow). These values of dimensionless parameters are likely to be seen in real applications. It can be seen that for different n values the linearization solution agrees very well with the numerical solu-

0

10

2

10

4

10

6

10

tD Fig. 2. Comparison of the numerical solution and linearization solution with the result of Hantush and Jacob (1955) for Darcian flow case.

tion at late times, and the linearization method underestimates the dimensionless drawdown at early times. This feature is understandable, because the linearization procedure is less accurate at early times, more accurate at late times, and sufficiently correct at steady-state. To understand why the linearization method underestimates the drawdown at early-time (see Fig. 3), one has to look at the nature of the linearization method. To discuss drawdown at a radial distance r from the pumping well, let us consider a vertically cylindrical domain of aquifer surrounding the well with a radius of r and the well as the axis. At any given time, one can always conceive that the total pumping rate Q consists of two components. One is water released from the aquifer storage within that cylinder, denoted as Qs, and the other comes from influx outside that cylinder, denoted as Qf, i.e., Q = Qs + Qf. At early times, most of the water comes from release of the aquifer storage, i.e., Qs will be dominating. When water is released from aquifer storage, water pressure drops or drawdown increases in proportional to Qs. Now let us check with the linearization approach (Eq. (15)). This approach is essentially a quasi steady-state approximation and it does not fully consider the water released from the aquifer storage which occurs at early times of pumping. In another word, this approach will result in an approximated Qs term that is smaller than the actual Qs term at early times. Physically, it means that this approach underestimates the release of aquifer storage water at early times. Because of such an underestimation, a smaller drawdown will be resulted from the linearization approach. Such an approximation will become increasingly accurate when time increases since the aquifer storage release (Qs) will continuously decline with time; and when the steady-state is reached, Qs will drop to zero and the linearization approach will become sufficiently correct. Another feature that can be found in this figure is that the linearization solution approaches the numerical solution earlier when the dimensionless distance rD is smaller. This is because it will take shorter time to approach the steady-state as the distance is closer to the pumping well. Comparing Fig. 3a–c, it seems that the linearization procedure works better when n is smaller. This is because the error of the linearization procedure is smaller when n is smaller, as shown in Eq. (15). We have also analyzed the effect of the dimensionless leakage parameter BD on the dimensionless drawdown, as shown in Fig. 4. The parameters are given as: n = 1.5, kD = 10, rD = 1, BD = 1,

1759

Z. Wen et al. / Advances in Water Resources 31 (2008) 1754–1763

a

b

1

10

Numerical solution Linearization solution

1

10

Numerical solution Linearization solution 0

0

10

10

r =1

rD=1

D

−1

−1

10

10

D −2

−2

10

10

−3

−3

10

10

n=1.1 k =10 D B =10 D

−4

10

r =10 D

s

s

D

rD=10

n=1.5 k =10 D B =10 D

−4

10 −2

0

10

2

10

10

tD

c

−2

4

0

10

10

2

10

10

4

10

tD

1

10

Numerical solution Linearization solution 0

10

rD=1 −1

10

s

D

r =10 D

−2

10

−3

10

n=2.0 k =10 D B =10 D

−4

10

−2

0

10

2

10

10

4

10

tD Fig. 3. Comparison of the numerical solution with linearization solution for kD = 10, BD = 10, rD = 1 and rD = 10, respectively: (a) n = 1.1; (b) n = 1.5 and (c) n = 2.0 (fully turbulent flow).

1

10

BD=1 BD=10 B =100 D

0

Wen et al. (2008d)

s

D

10

−1

10

−2

10

n=1.5 k =10 D r =1 D

−3

10

−2

10

0

2

10

10

4

10

tD Fig. 4. Dimensionless drawdown versus dimensionless time without wellbore storage for n = 1.5, kD = 10, rD = 1 with BD = 1, 10, and 100, respectively.

10, and 100, respectively. It should be pointed out that we only analyze the numerical solution hereinafter for different parameters. Similar conclusions can be drawn from the linearization solutions. The corresponding solution without considering leakage effect obtained by Wen et al. [50] has also been depicted in this figure as a reference. The results indicate that all the curves for different BD values approach the same asymptotic value as that without leakage at early times. This means the leakage water has not arrived at the main aquifer yet at early times, and the pumping rate comes from the main aquifer only during this pumping period. While at late times, the dimensionless drawdown is smaller as BD is smaller. This is because a smaller BD implies a greater leakage effect, which ends in a smaller drawdown in the aquifer. All these features are similar to those reported by Hantush and Jacob [13] for Darcian flow case. The sensitivity analysis of the dimensionless non-Darcian conductivity kD is shown in Fig. 5 with n = 1.5, BD = 10, rD = 1, and kD = 10, 20, 50, and 100, respectively. The solutions with leakage effect for the Darcian flow case obtained by Hantush and Jacob [13] have also been included in this figure as a reference. The dimensionless hydraulic conductivity kD for the solution of Hantush and Jacob [13] is given as 1, because the power index n is equal to 1 for the Darcian flow. One can see that a larger kD results

1760

Z. Wen et al. / Advances in Water Resources 31 (2008) 1754–1763 1

1

10

10 k =10 D k =20 D k =50 D k =100

Steady−state solution tD=10 tD=1 t =0.1

0

D

D

10

Hantush and Jacob (1955)

n=1.5 kD=10 BD=10

0

s

D

sD

10

−1

10

−1

10

−2

10

n=1.5 BD=10 r =1 D

−2

−3

10

10

−2

10

0

10

−1

2

0

10

10

1

10

rD

tD

5.3. Dimensionless drawdown-distance behavior without wellbore storage

1

10

With wellbore storage Without wellbore storage r =0.1 D

0

10

rD=1

D

in a larger drawdown at early times, and leads to a smaller drawdown at late times. All the drawdowns for the non-Darcian flow case are larger at early times and smaller at late times than their counterparts for Darcian flow case. Actually, a larger kD means that the aquifer can transmit water easier, which will yield a smaller drawdown at late times. The feature of the early-time dimensionless drawdown can be explained as follows. At early times, the dimensionless discharge qD increases frome zero to its steady-state value rather than a constant, as the dimensionless time tD increases. For a larger kD value, qD increases with tD faster to its steady-state value. Thus, at a given time tD at early times, a larger kD will result in a greater qD. Then from Eq. (12), it is possible to see that a larger sD will be found for a larger kD. Similar results have also been found by Wen et al. [50] when the leakage is excluded.

−1

10

r =10 D

n=1.5 kD=10 BD=10

−2

10

S=0.001 r =0.1 wD r =0.1 cD

−3

10

−2

10

It is also important to analyze the dimensionless drawdowndistance behavior because sometimes it is possible to measure the drawdowns at different observation wells at the same time for a pumping test. The dimensionless drawdown-distance curves for different dimensionless time tD are given in Fig. 6. The steady-state solution, Eq. (28), is also depicted as a reference. The parameters are given as n = 1.5, kD = 10, BD = 10, tD = 0.1, 1, and 10, respectively. It can be seen that when tD is equal to 10, the curve is very close to the steady-state solution. This figure can be used to estimate the parameters of the aquifer when the drawdown-distance data are known. 5.4. Dimensionless drawdown-time behavior with wellbore storage In the following analysis, we use the same value for dimensionless radius of well casing rcD and dimensionless radius of well screen rwD although they are sometimes different for real cases. The results related to different rcD and rwD values can also been easily obtained. First, we compare the dimensionless drawdown with that without wellbore storage, as shown in Fig. 7. The parameters are given as n = 1.5, BD = 10, kD = 10, rwD = rcD = 0.1, S = 0.001 and rD = 0.1 (at the wellbore), 1 and 10, respectively. It is obvious to see when the pumping time is long enough, the dimensionless drawdowns with and without wellbore storage converge to a single curve. This

10

Fig. 6. Dimensionless drawdown versus dimensionless distance without wellbore storage for n = 1.5, BD = 10,kD = 10 with tD = 0.1, 1, and 10, respectively.

s

Fig. 5. Dimensionless drawdown versus dimensionless time without wellbore storage for n = 1.5, BD = 10, rD = 1 with kD =10, 20, 50, and 100, respectively.

2

10

0

2

10

10

4

10

t

D

Fig. 7. Comparison of the dimensionless drawdown with and without wellbore storage for n = 1.5, BD = 10,kD = 10, S = 0.001, rwD = 0.1, rcD = 0.1 with rD = 0.1, 1, and 10, respectively.

means the wellbore storage has no effect at large times. In addition, one can see as the dimensionless distance rD increases, the differences between the two dimensionless drawdowns become smaller and smaller. We also compare the linearization solution with the numerical solution with n = 1.5, BD = 10, kD = 10, rwD = rcD = 0.1, S = 0.001 for three typical distances rD = 0.1 (inside the well), rD = 1 and rD = 5, respectively, as shown in Fig. 8. It can also be found that the linearization procedure works very well at late times for both the dimensionless drawdowns inside the well and in the aquifer. Interestingly, the linearization procedure slightly overestimates the dimensionless drawdown at moderate times for rD = 0.1 (at the wellbore); slightly underestimates the dimensionless drawdown at early times and slightly overestimates the dimensionless drawdown at moderate times for rD = 1; and slightly underestimates the dimensionless drawdown at early and moderate times for rD = 5, when comparing to the numerical solutions in Fig. 8. The finite difference grid used in the numerical solution is sufficiently fine that the errors of the numerical solutions are probably negligible. The

1761

Z. Wen et al. / Advances in Water Resources 31 (2008) 1754–1763 1

60

10

Numerical solution Linearization solution

rD=0.1 r =0.5 D r =1

40

r =0.1 D

D

20

0

10

r =1

0

r =5

−20

D

−1

10

ε

sD

D

−40

−2

10

n=1.5 k =10 D B =10

−60

n=1.5 B =10 D k =10

D

−80

D

S=0.001 r =0.1 wD r =0.1

S=0.001 r =0.1 wD r =0.1

−100

cD

cD

−3

−120

10

−2

10

0

2

10

10

−2

4

0

10

10

10

D

Fig. 9. Relative error of the linearization solution for n = 1.5, BD = 10, kD = 10, S = 0.001, rwD = 0.1, rcD = 0.1 with rD = 0.1, 0.5, and 1, respectively.

The relative error,e, is defined as the ratio of the difference between the linearization solution and the numerical solution over the numerical solution. The parameters are given as n = 1.5, BD = 10, kD = 10, S = 0.001, rwD = 0.1, rcD = 0.1 with rD = 0.1, 0.5, and 1, respectively. It is obvious to see that the relative error is very small at late times for different rD values. However, at early times, the relative error appears to be large when the dimensionless distance rD is relatively large. At very early times, the relative error is almost equal to 100%, as can be seen for rD = 0.5 and rD = 1. As the pumping time increases, the relative error decreases gradually and approaches to zero when the pumping time is long enough. It is notable that at early times, the drawdown itself is very small, thus even a very small absolute difference might result in a large relative error. From a practical point of view, the early-time solution is probably less important than the late-time one, as the early-time drawdown is small and its measurements may be suffered from experimental errors. It is probably more reliable to estimate the parameters of the aquifer with the late-time solution. In this case, the linearization solution might be quite useful.

1

10

BD=1 B =10 D BD=100 Wen et al. (2008d) 0

10

wD

accuracy of the numerical solutions is also reflected in Fig. 2 for the case when n = 1. Therefore, the discrepancies between the numerical solutions and the linearization solutions observed in Fig. 8 are more likely to be caused by the linearization approximation. In the following, we explain why the drawdown-time curves in Fig. 8 should be that way. First, let us discuss the case of rD = 0.1, which is at the wellbore (notice that rwD is also 0.1). At early times, the pumping rate almost exclusively comes from the wellbore storage independent of the aquifer, so both the numerical solution and the linearization method yield nearly the same drawdown, shown as a straight line in Fig. 8. When time proceeds, the effect of the wellbore storage gradually declines and contribution from the aquifer continuously increases. As pointed out before, the nature of the linearization approach is essentially a quasi steady-state approximation. Physically, it treats the aquifer as a rather rigid medium with a much smaller storativity than its actual value. As a result, it leads to a quicker transaction from the early-time transient curve to the late-time or steady-state curve. Graphically, this is manifested by slightly greater drawdowns at moderate times and a quicker reach of the steady-state asymptote. Now let us see what happens for the case of rD = 5, which is rather far from the well. Because this point is far from the wellbore, the type curves are much less sensitive to the wellbore storage effect except that both curves (one is from the numerical solution and one is from the linearization method) are shifted rightward along the time axis when compared with their counterparts without the wellbore storage, reflecting the delayed aquifer response to pumping. The overall shapes of the type curves of this case are similar to what have been observed in their counterparts without the wellbore storage (see Fig. 3), i.e., the linearization method often yields a smaller drawdown than the numerical method, as explained in Section 5.2. For the case of rD = 1, which is somewhat between the two cases discussed above (rD = 0.1 and rD = 5), a somewhat mixing behavior of the type curves has been observed in Fig. 8. That is: at early times, the linearization method underestimates the drawdown while at moderate times, it overestimates the drawdown. This case indeed reflects the transaction from a near-well case of rD = 0.1 to a far-well case ofrD = 5. In order to quantify the error of the linearization procedure, we have plotted the relative error in semi-log scale, as shown in Fig. 9.

4

10

t

s

Fig. 8. Comparison of the dimensionless drawdown obtained by the linearization procedure with that obtained by the numerical method for n = 1.5, BD = 10, kD = 10, S = 0.001, rwD = 0.1, rcD = 0.1 at three different dimensionless distances rD = 0.1 (inside the well), rD = 1 and rD = 5, respectively.

2

10

tD

−1

n=1.5 kD=10

10

S=0.001 rwD=0.1 r =0.1 cD r =0.1 D

−2

10

−2

10

0

2

10

10

4

10

t

D

Fig. 10. Dimensionless drawdown versus dimensionless time with wellbore storage for n = 1.5, kD = 10, S = 0.001, rwD = 0.1, rcD = 0.1, rD = 1 with BD = 1, 10, and 100, respectively.

1762

Z. Wen et al. / Advances in Water Resources 31 (2008) 1754–1763 1

10

k =10 D kD=20 kD=50 kD=100

0

s

wD

10

−1

10

n=1.5 B =10 D

S=0.001 rwD=0.1 rcD=0.1 r =0.1 D

−2

10

−2

10

0

2

10

10

4

10

t

D

Fig. 11. Dimensionless drawdown versus dimensionless time with wellbore storage for n = 1.5, BD = 10, S = 0.001, rwD = 0.1, rcD = 0.1, rD = 1 with kD = 10, 20, 50, and 100, respectively.

When the wellbore storage is considered, it is important to analyze the drawdown inside the well. The effect of the leakage on the dimensionless drawdown inside the well is shown in Fig. 10. The parameters are given as n = 1.5, kD = 10, rD = rwD = rcD = 0.1, S = 0.001, BD = 1, 10, and 100, respectively. The corresponding solution without considering leakage obtained by Wen et al. [50] is also depicted in this figure as a reference. Several observations can be drawn from this figure. First, when BD is large enough, the dimensionless drawdown approaches the same value as that without leakage. The reason for this has been explained before. Second, all the curves for different BD values approach the same value at early times and are straight lines in log–log scale. This is the obvious effect of wellbore storage. Third, a larger BD value, which means a smaller leakage effect, results in a large drawdown inside the well at late times. Similarly, we have also analyzed the effect of the dimensionless non-Darcian conductivity kD on the dimensionless drawdown inside the well, as shown in Fig. 11. The parameters are given as n = 1.5, BD = 10, rD = rwD = rcD = 0.1, S = 0.001, kD = 10, 20, 50 and 100, respectively. At early times, all the curves converge to a single straight line in log–log scale. Once again, this is the wellbore storage effect. At late times, a larger kD results in a smaller drawdown. We have also analyzed the dimensionless drawdown in the aquifer when the wellbore storage is considered. It can be found that the features are similar to those without wellbore storage. Thus, the figures about the dimensionless drawdown in the aquifer will not be repeated here. 6. Discussion It is necessary to discuss the physical system of the mathematical model proposed here. We used the power law (Izbash equation) to describe the non-Darcian flow in the aquifer. In reality, the power index n in the Izbash equation might be changeable instead of a constant for the entire flow domain. This is because the flow velocity in the aquifer decreases as the distances increases. It is very difficult, if not impossible to capture such a continuously variable flow field. Thus, approximation might be inevitable for the analytical treatment. There are certain evidences to show that the Izbash equation might work very well in an approximate sense. As explained in our previous paper [48], a careful check of the Moody diagram of flow in rounded pipes with different surface

roughness indicates that the function of the so-called ‘‘friction factor” versus the Reynolds number can be approximated by a straight line with a slope other than unity when the relative roughness is small [28, p. 493]. Substituting this friction factor into the Darcy–Weisbach equation will lead to the Izbash type of non-Darcian flow equation [28, p. 491]. From the Moody diagram [28, p. 493], one can see the curve is nearly a straight line even the Reynolds number changes over several magnitudes. Therefore, it might be possible to use the Izbash equation to approximate the flow over a large spatial scale. The Izbash equation offers some advantages for obtaining an analytical solution, as shown in this study. For the case when the Izbash equation is incapable of approximating the non-Darcian flow within an acceptable accuracy, one should try a different non-Darcian flow equation such as the Forchheimer equation. If one can find the Forchheimer equation with two ‘‘constant” coefficients, one should use it instead of the Izbash equation. One of the advantages of the Forchheimer equation is that the second-order term in that equation is negligible when the flow velocity is very small, leading to the Darcian flow case with sufficiently small velocity. This can be done within the same framework of this study when a numerical approach is used. Mathias et al. [26] have done some work in this regard. Regarding the analytical solutions, even an approximate one with acceptable accuracy for the Forchheimer non-Darcian flow to a well in a leaky aquifer, we still do not find any from the literature.

7. Summary and conclusions In this study, we have used a linearization procedure and a finite difference method to solve non-Darcian flow to a well in an aquifer–aquitard system. The leakage effect is included in the mathematical model. We have also considered the wellbore storage effect. The Izbash equation has been used to describe non-Darcian flow in the aquifer, whereas flow in the aquitard is assumed to be Darcian. Semi-analytical solution and numerical solution for such a case have been obtained. The previous solution of Wen et al. [50] in which the leakage is not considered and the solution of Hantush and Jacob [13] who considered a similar problem based on Darcian flow in the aquifer, can also be described as a special case of the new solutions. We have compared our results with that of Wen et al. [50] who investigated a similar problem without considering leakage. The comparison of our results and the solution of Hantush and Jacob [13] is also included. Several findings can be drawn from this study: (1) The relative error of the linearization procedure is nearly 100% at early times, while it decreases to zero at late times. (2) A larger dimensionless non-Darcian conductivity kD results in a smaller dimensionless drawdown in the aquifer at late times, while it leads to a larger dimensionless drawdown in the aquifer at early times. (3) The dimensionless drawdown in the aquifer for non-Darcian flow with leakage in this study is larger at early times and smaller at late times than that for Darcian flow case [13]. (4) A smaller dimensionless leakage parameter BD results in a smaller drawdown at late times, and the leakage does not affect the early-time drawdown. Acknowledgements This research was partially supported by the National Natural Science Foundation of China (grant numbers: 50639040, 50428907, 50779067), the Elitist Program of the Ministry of Education of China (grant number: NCET-05-0125) and the Program for

Z. Wen et al. / Advances in Water Resources 31 (2008) 1754–1763

Changjiang Scholars and Innovative Research Team in University. We would like to thank Dr. Simon A. Mathias for his help on the MATLAB programs. The constructive comments from five anonymous reviewers and the editor are also gratefully acknowledged. References [1] Bakker M. Simulating groundwater flow to surface water features with leaky beds using analytic elements. Adv Water Resour 2007;30:399–407. [2] Birpinar ME, Sen Z. Forchheimer groundwater flow law type curves for leaky aquifers. J Hydrol Eng 2004;9(1):51–9. [3] Bordier C, Zimmer D. Drainage equations and non-Darcian modeling in coarse porous media or geosynthetic materials. J Hydrol 2000;228:174–87. [4] Brinkman HC. A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles. Appl Sci Res 1947;1:24–37. [5] Butler JJ, Liu W, McElwee CD. Type curves for two-regime well flow – discussion. ASCE J Hydraul Eng 1991;117(1):122–4. [6] Camacho RG, Vasquez M. Comment on ‘‘Analytical Solution Incorporating Nonlinear Radial Flow in Confined Aquifers” by Zekai Sen. Water Resour Res 1992;28(12):3337–8. [7] Choi ES, Cheema T, Islam MR. A new dual-porosity/dual-permeability model with non-Darcian flow through fractures. J Petrol Sci Eng 1997;17:331–44. [8] Crump KS. Numerical inversion of Laplace transforms using a Fourier-series approximation. J ACM 1976;23(1):89–96. [9] de Hoog FR, Knight JH, Stokes AN. An improved method for numerical inversion of Laplace transforms. SIAM J Sci Stat Comput 1982;3(3):357–66. [10] Ewing RE, Lazarov RD, Lyons SL, Papavassiliou DV, Pasciak J, Qin G. Numerical well model for non-Darcy flow through isotropic porous media. Comput Geosci 1999;3:185–204. [11] Ewing RE, Lin Y. A mathematical analysis for numerical well models for nonDarcy flows. Appl Numer Math 2001;39(1):17–30. [12] Forchheimer PH. Wasserbewegun durch Boden. Zeitschrift des Vereines Deutscher Ingenieure 1901;49:1736–49. and 50:1781–1788. [13] Hantush MS, Jacob CE. Non-steady radial flow in an infinite leaky aquifer. Trans Am Geophys Union 1955;36(1):95–100. [14] Hantush MS. Nonsteady flow to flowing wells in leaky aquifers. J Geophys Res 1959;64(8):1043–52. [15] Hantush MS. Modification of the theory of leaky aquifers. J Geophys Res 1960;65(11):3713–25. [16] Hunt B, Massmann J. Vapor flow to trench in leaky aquifer. J Environ Eng 2000;126(4):375–80. [17] Izbash SV. O filtracii V Kropnozernstom Materiale. USSR: Leningrad; 1931 [in Russian]. [18] Kohl T, Evans KF, Hopkirk RJ, Jung R, Rybach L. Observation and simulation of non-Darcian flow transients in fractured rock. Water Resour Res 1997;33(3):407–18. [19] Kolditz O. Non-linear flow in fractured rock. Int J Numer Methods Heat Fluid Flow 2001;11(6):547–75. [20] Lee JY, Lee KK. Analysis of the quality of parameter estimates from repeated pumping and slug tests in a fractured porous aquifer system in Wonju, Korea. Ground Water 1999;37(5):692–700. [21] Legrand J. Revisited analysis of pressure drop in flow through crushed rocks. J Hydraul Eng 1999;128(11):692–700. [22] Li J. Analysis of radial movement of an unconfined leaky aquifer due to well pumping and injection. Hydrogeol J 2007;15:1063–76. [23] Li Y, Neuman SP. Flow to a well in a five-layer system with application to the Oxnard Basin. Ground Water 2007;45(6):672–82. [24] Malama B, Kuhlman KL, Barrash W. Semi-analytical solution for flow in leaky unconfined aquifer–aquitard systems. J Hydrol 2007;346:59–68. [25] Martys NS. Improved approximation of the Brinkman equation using a lattice Boltzmann method. Phys Fluids 2001;13(6):1807–10. [26] Mathias SA, Butler AP, Zhan H. Approximate solutions for Forchheimer flow to a well. ASCE J Hydraul Eng 2008;134:1318–25. [27] Moench AF. Convergent radial dispersion: a note on evaluation of the Laplace transform solution. Water Resour Res 1991;27(12):3261–4.

1763

[28] Munson BR, Young DF, Okiishi TH. Fundamental of fluid mechanics. 3rd ed. New York: John Wiley & Sons; 1998. [29] Neuman SP, Witherspoon PA. Applicability of current theories of flow in leaky aquifers. Water Resour Res 1969;5(4):817–29. [30] Perina T, Lee TC. General well function for pumping from a confined, leaky, or unconfined aquifer. J Hydrol 2006;317:239–60. [31] Ramey HJ. Approximate solutions for unsteady liquid flow in composite reservoirs. J Can Petrol Technol 1970;9(1):32–7. [32] Sen Z. Non-Darcian flow in fractured rocks with a linear flow pattern. J Hydrol 1987;92:43–57. [33] Sen Z. Type curves for two-region well flow. J Hydraul Eng 1988;114(12):1461–84. [34] Sen Z. Analytical solution incorporating nonlinear radial flow in confined aquifers. Water Resour Res 1988;24(4):601–6. [35] Sen Z. Nonlinear flow toward wells. J Hydraul Eng 1989;115(2): 193–209. [36] Sen Z. Nonlinear radial flow in confined aquifers toward large-diameter wells. Water Resour Res 1990;26(5):1103–9. [37] Sen Z. Non-Darcian groundwater flow in leaky aquifers. Hydrol Sci J 2000;45(4):595–606. [38] Shampine LF, Reichelt MW. The MATLAB ODE Suite. SIAM J Sci Comput 1997;18:1–22. [39] Shampine LF, Reichelt MW, Kierzenka JA. Solving Index-1 DAEs in MATLAB and Simulink. SIAM J Sci Comput 1999;41:538–52. [40] Sidiropoulou MG, Moutsopoulos KN, Tsihrintzis VA. Determination of Forchheimer equation coefficients a and b. Hydrol Process 2007;21: 534–54. [41] Singh JP, Shakya SK. Drainage of ponded land through irrigation well in a leaky aquifer. J Hydrol 2006;321:377–89. [42] Spanier J, Oldham KB. An atlas of functions. New York, USA: Hemisphere Publishing; 1987. [43] Stehfest H. Algorithm 368 numerical inversion of Laplace transforms. Commun ACM 1970;13(1):47–9. [44] Stehfest H. Remark on algorithm 368 numerical inversion of Laplace transforms. Commun ACM 1970;13(10):624–5. [45] Talbot A. The accurate numerical inversion of the Laplace transforms. J Inst Math Appl 1979;23:97–120. [46] Wen Z, Huang G, Zhan H. Non-Darcian flow in a single confined vertical fracture toward a well. J Hydrol 2006;330:698–708. [47] Wen Z, Huang G, Zhan H. An analytical solution for non-Darcian flow in a confined aquifer using the power law function. Adv Water Resour 2008;31:44–55. [48] Wen Z, Huang G, Zhan H, Li J. Two-region non-Darcian flow toward a well in a confined aquifer. Adv Water Resour 2008;31:818–27. [49] Wen Z, Huang G, Zhan H. Non-Darcian flow toward a finite-diameter vertical well in a confined aquifer. Pedosphere 2008;18(3):288–303. [50] Wen Z, Huang G, Zhan, H. A numerical solution for non-Darcian flow to a well in a confined aquifer using the power law function. J Hydrol 2008, accepted for publication. [51] Wu YS. An approximate analytical solution for non-Darcy flow toward a well in fractured media. Water Resour Res 2002;38(3):1023. doi:10.1029/ 2001WR000713. [52] Wu YS. Numerical simulation of single-phase and multiphase non-Darcy flow in porous and fractured reservoirs. Transport Porous Med 2002;49: 209–40. [53] Yamada H, Nakamura F, Watanabe Y, Murakami M, Nogami T. Measuring hydraulic permeability in a streambed using the packer test. Hydrol Process 2005;19:2507–24. [54] Yeh HD, Lin YC, Huang YC. Parameter identification for leaky aquifers using global optimization methods. Hydrol Process 2007;21:862–72. [55] Zhan H, Bian A. A method of calculating pumping induced leakage. J Hydrol 2006;328:659–67. [56] Zlotnik VA. A concept of maximum stream depletion rate for leaky aquifers in alluvial valleys. Water Resour Res 2004;40:W06507. doi:10.1029/ 2003WR002932. [57] Zlotnik VA, Tartakovsky DM. Stream depletion by groundwater pumping in leaky aquifers. J Hydrol Eng 2008;13(2):43–50.