Recent advances in modeling of well hydraulics

(1)

Recent advances in modeling of well hydraulics

Hund-Der Yeh

a,⇑

, Ya-Chi Chang

a,b

a

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

Taiwan Typhoon and Flood Research Institute, National Applied Research Laboratories, Taiwan

a r t i c l e

i n f o

Article history:

Available online 15 March 2012 Keywords: Groundwater Aquifer test Well hydraulic Analytical model Numerical model

a b s t r a c t

Well hydraulics is a discipline to understand the process of flow to the well in an aquifer which is regarded as a source of groundwater. A variety of analytical and numerical models have been developed over the last few decades to provide a framework for understanding and quantifying the flow behavior in aquifer systems. In this review, we first briefly introduce the background of the theory of well hydraulics and the concepts, methodologies, and applications of analytical, semi-analytical, numerical and approx-imate methods in solving the well-hydraulic problems. We then address the subjects of current interests such as the incorporation of effects of finite well radius, wellbore storage, well partial penetration, and the presence of skin into various practical problems of groundwater flow. Furthermore, we also summa-rize recent developments of flow modeling such as the flow in aquifers with horizontal wells or collector wells, the capture zone delineation, and the non-Darcian flow in porous media and fractured formations. Finally, we present a comprehensive review on the numerical calculations for five well functions fre-quently appearing in well-hydraulic literature and suggest some topics in groundwater flow for future research.

1. Introduction

1.1. Background

Groundwater is an important alternative resource to surface water for agriculture, industry and domestic use. Studies in groundwater hydrology are therefore devoted to evaluate the occurrence, availability and quality of groundwater. Two inherent characteristics of the aquifer, referred to as the specific storage and hydraulic conductivity, generally provide the fundamental bases for the quantitative studies on groundwater hydrology/ hydraulics. These two parameters are commonly determined from the field data of aquifer tests. Many analytical solutions have been developed rapidly and steadily since Theis published the solution of transient flow equation, a milestone in well hydraulics, in 1935[1]. Some analytical solutions were used to develop the type curves so that the hydraulic parameters can be determined by the graphical fitting of the observed data to the type curves[2]. The analytical solutions are however restricted to the ideal cases, which should be under the conditions of simple aquifer boundary and homogeneous formation properties. Yet, the real-world well-hydraulic problems are often with complicated boundary

and/or heterogeneous aquifer properties. Under these circum-stances, numerical approaches such as finite-difference methods, finite element methods, or boundary element methods are then employed to develop the numerical solutions. Numerical ap-proaches have been extensively used in studying groundwater problems since the mid-1960s[3]. Often, newly developed numer-ical models are verified with analytnumer-ical models to check their correctness or accuracy. Analytical models are often developed for better modeling the aquifer systems in different subjects by accounting for the effects of wellbore storage, skin zone, and well partial penetration on the groundwater flow. In last two decades, many analytical solutions arisen from various types of aquifer tests with considering those effects have been developed. Some ground-water issues, such as pollution and ground-water resources management, also attract considerable concern of scientists and engineers. Many analytical models for describing groundwater flow induced by hor-izontal and collecting wells and for delineating the capture zone in a contaminated site are then developed in accordance with these issues. Many of the analytical solutions are in terms of complicated forms which may contain integrals of some special functions, e.g., Bessel functions or trigonometric functions, and these integrals are difficult to compute accurately because of the oscillatory nature and slow convergence in computing the integrand. Therefore, var-ious techniques are proposed to numerically calculate the solu-tions with high accuracy.

In this review, we address the subjects of current interests in re-gard to modeling groundwater ﬂow behavior in well hydraulics.

http://dx.doi.org/10.1016/j.advwatres.2012.03.006

⇑ Corresponding author. Address: 300 Institute of Environmental Engineering, National Chiao Tung University, No. 1001, University Road, Hsinchu, Taiwan.

E-mail addresses: [email protected] (H.-D. Yeh), rachel-chang@ttfri. narl.org.tw(Y.-C. Chang).

Contents lists available atSciVerse ScienceDirect

Advances in Water Resources

(2)

For each of these subjects, we provide a brief and historical over-view, summarize the recent developments, and suggest some top-ics for future research as well. In this chapter, we first present the groundwater flow equation and its associated boundary conditions in well hydraulics, then introduce the definitions of different types of aquifers, and finally review various aquifer tests in some detail.

1.2. Basic principle, ﬂow equations and initial and boundary conditions

1.2.1. Darcian ﬂow

Darcian flow, which obeys Darcy’s law and describes laminar flow behavior in a porous medium, can be expressed as q = –Ki in which q is the discharge velocity or Darcy velocity, K is the lic conductivity, h is the hydraulic head, and i = dh/dl is the hydrau-lic gradient. Instead of describing the flow state within individual pores, Darcy’s law represents the statistical macroscopic equiva-lent of the Navier–Stokes equations for viscous flow in porous media.

1.2.2. Three dimensional groundwater ﬂow equation

Based on Darcy’s law, a three-dimensional (3-D) equation in Cartesian coordinates for groundwater ﬂow in a heterogeneous and anisotropic aquifer can be written as

@ @x Kxðx; y; zÞ @h @x þ@ @y Kyðx; y; zÞ @h @y þ@ @z Kzðx; y; zÞ @h @z ¼ Ss @h @tþ R ð1Þ

where Kx, Ky, and Kzare the main components of the hydraulic

con-ductivity tensor, Ssis the speciﬁc storage, R is a recharge term which

is positive and defined as the volume of inflow to the flow system per unit volume of aquifer per unit of time, and t is the time from the start of test.

1.2.3. Radial ﬂow equation

Radial ﬂow toward the well or from the well occurs when a ver-tical extraction or injection is performed at the well in the aquifer. The differential equation for ﬂow in aquifers in cylindrical coordi-nates can be written as

Kr @2h @r2þ Kr r @h @rþ Kz @2h @z2¼ Ss @h @t ð2Þ

where r is the radial distance from the center of the pumping well.

1.2.4. Initial and boundary conditions

The initial condition should be speciﬁed when solving transient ﬂow problems in which the hydraulic head changes with time. The mathematical expression for initial condition is denoted as

h ¼ f ðx; y; z; 0Þ ð3Þ

It is important to select appropriate boundary conditions in devel-oping a mathematical model for an aquifer flow system. Boundary conditions for groundwater flow problems are of three types, i.e., the first-type (Dirichlet), the second-type (Neumann), and the third-type (Robin) boundary conditions.

A constant hydraulic head is the ﬁrst-type boundary which can be mathematically expressed as

hðx; y; z; tÞj_C¼ f1ðx; y; z; tÞ ð4Þ

where f1is a known function. Physically, a surface water body, e.g.,

river, lake, or reservoir could be treated as a constant head bound-ary if they hydraulically connect with the aquifer.

The mathematical expression of the water ﬂux at the boundary

Cis

qnðx; y; z; tÞjC¼ f2ðx; y; z; tÞ ð5Þ

where qnis a ﬂux component normal to the boundary surface. An

example of the second-type boundary can be selected at the top of an aquifer where there is recharge or discharge. The no-ﬂow boundary (qn= 0) is a special case of the second-type boundary

condition.

The third-type boundary condition is a linear combination of Dirichlet and Neumann boundary conditions, which relates hydraulic head to the ﬂux and can be mathematically expressed as

ahðx; y; z; tÞjCþ bqnðx; y; z; tÞjC¼ gðx; y; z; tÞ; ðx; y; zÞ

C

ð6Þ

In well hydraulics, an example of this kind of boundary is the case when taking into account of the skin effect in the wellbore as men-tioned in Section3.2.

It is worthy to mention that the mixed boundary value problem arise when a boundary condition of different types speciﬁed on dif-ferent subsets of the boundaryC. IfCis divided into two subsets

C1andC2, this mixed kind of boundary condition may be written

as

hðx; y; z; tÞjC1¼ g1ðx; y; z; tÞ ð7Þ

qnðx; y; z; tÞjC2¼ g2ðx; y; z; tÞ ð8Þ

The mixed boundary value problem is commonly encountered in the physical problems of ﬂuid dynamics, electricity, heat ﬂow, and others. In well hydraulics, it may appear in the case of constant head pumping in a well under partial penetration condition. Some issues relative to the mixed boundary value problem can refer to Section3.3.

1.2.5. Aquifer of finite or infinite extent in horizontal direction Analytical solutions can be used to predict the temporal or spa-tial drawdown distribution or to determine aquifer parameters. For instance, the most well-known transient and steady-state draw-down solutions for pumping with a constant rate conducted in confined aquifers are Theis and Thiem solutions, respectively. As indicated in Wang and Yeh[4], it is inappropriate to apply Theis equation to the case where an aquifer has a horizontally finite boundary or the pumping time tends to infinity. Most of existing analytical solutions in the well-hydraulic literature are developed by assuming that the aquifer is of infinite extent in horizontal direction. Studies on the groundwater problems in an aquifer with a horizontally finite boundary are very limited[5].

It is very common to consider that the aquifer is of ﬁnite extent in vertical direction in modeling groundwater problems. Some studies however treated the aquifer thickness as inﬁnite because the aquifer thickness is relatively large when compared with the screen depth of the well[6].

The analytical solutions may be developed for flow in finite aquifers with a regular shape such as rectangular and wedge-shaped configuration or irregular shape such as step change configuration. Several studies provided analytical solutions for groundwater flow problems in aquifers whose plan views are rect-angular[7–9], wedge-shaped[10–15]and triangular[16], and in aquifers whose cross sectional views are step-like[17,18].

1.3. Types of aquifer

An aquifer is a geological formation that is sufficiently perme-able to store and transmit groundwater. Related terminologies in-clude aquitard (also called as semipervious layer), defined as a geologic unit that has relatively low permeability compared with aquifers, and aquiclude, which is a formation not capable of trans-mitting a significant amount of groundwater. A groundwater flow system can be composed of different types of aquifers.

(3)

A confined aquifer, also known as an artesian or a pressure aqui-fer, is bounded both at the top and bottom by impervious or semipervious strata and thus hydraulically isolated from other geological formations. A piezometer or observation well installed in the confined aquifer has a hydrostatic pressure level which forms an imaginary surface called a piezometric surface. Generally, water in the piezometer or well rises above the top of the aquifer. An unconfined aquifer is also called as a water table or phreatic aquifer. The upper boundary of the unconfined aquifer is a free sur-face called the water table.

An aquifer that loses or gains significant water through the overlain, underlain, or both semipermeable strata is called leaky aquifers. It is common that an aquifer hydraulically connects with other formations in nature. A system of multilayered aquifers (or multiple aquifers) is composed by a series of aquifers separated from each other by semipervious layers. The mechanics of fluid flow through a multiple aquifer system becomes complicated be-cause of the hydraulic connection between individual aquifers. For groundwater flow associated with leaky and multilayered aqui-fers, the reader is referred to Section3.8.

1.4. Aquifer tests

Two different types of pumping test are commonly used in the ﬁeld for estimating aquifer and well characteristics; they are aqui-fer test and well test[19]. The aquifer test is used to determine the aquifer hydraulic parameters. On the other hand, the well test is utilized to provide the information about the speciﬁc capacity or the discharge–drawdown ratio of the well[20]. The aquifer test generally needs at least one observation well or piezometer to re-cord the water level response to the pumping, while the well test may require only the analysis of pumped well data. Kruseman and de Ridder [20] presented comprehensive discussion on the general set-up regarding the selection of test site, the design and construction of discharging well, and the performance including the water level and discharge rate measurements of the aquifer test. In addition, Driscoll’s book [21] provides a good reference on some practical well-hydraulic issues as well as some aspects involving well design, drilling, construction, maintenance and rehabilitation. In the petroleum literature, Gringarten[22]and oth-ers[23–25]gave the review of the evolution of aquifer test analysis over the past half century. The aquifer test may be further divided into four groups: constant rate test, constant head test, slug test, and recovery test.

1.4.1. Constant rate test

The constant rate test (CRT) (or constant discharge test) is the most common form of aquifer test in which a test well is pumped at a constant rate over a certain period of time. At the same time, the drawdown is measured in one or more nearby observation wells[19]. The flow rate to the well is general assumed constant and equal to the average value of the pumping rates which indeed vary with time and are difficult to maintain constant in the field tests. In the early stage of development for flow induced by pump-ing in confined aquifers, there are two classic articles devoted to that in well hydraulic literature. Thiem[26]was the first to derive a steady-state expression for a fully penetrating well with a con-stant rate pumping in a confined aquifer. Theis[1] analyzed the groundwater flow in a homogeneous, isotropic and infinitely radial extent confined aquifer with a fully penetrating well pumped at a constant rate. The diameter of the pumping well is infinitesimal and the well storage is negligible. In recent years, Yeh et al.[27]

presented a mathematical model for an aquifer having a fully pene-trating well with a finite-thickness skin and constant pumping rate. Later, Yang et al.[28] developed an analytical solution for transient flow in a confined aquifer with a partially penetrating

well pumped at a constant rate. For unconﬁned aquifers, Jacob

[29]adopted the Dupuit assumption of no vertical ﬂow to derive the steady-state drawdown in CRT. Boulton[30] introduced the drawdown equation of the water table near a pumped well based on the delayed response concept of unconﬁned aquifers. Neuman

[31] obtained a transient drawdown solution by shifting the boundary condition from the free surface to the horizontal plane at a fixed position. Moench[32]developed a semianalytical solu-tion by combining the Boulton and Neuman models for flow to a partially penetrating well in unconfined aquifers. Tartakovsky and Neuman [33] presented an analytical expression for draw-down in an unconfined aquifer caused by the pumping at a con-stant rate from a partially penetrating well. Pasandi et al. [34]

provided an analytical model for CRT conducted at a partially pene-trating well with a finite thickness skin in an unconfined aquifer. For leaky and multilayered aquifers, for instance, Hantush and Ja-cob[35]developed a mathematical model for aquifer dynamics un-der a transient CRT by assuming that the confined aquifer is overlain everywhere by a semipervious layer which has constant vertical hydraulic conductivity and a constant thickness. Moench

[36]took into account the effect of large-diameter well and devel-oped mathematical models for flow in an aquifer system, where semipervious layers are located above and below the main pumped aquifer. In addition, there are some other studies associated with the CRT in confined aquifers [17,37–40], in unconfined aquifers

[41,42], and in leaky and multilayered aquifers[39,43–46].

1.4.2. Constant head test

For a constant head test (CHT), a constant water level in the test well is maintained, while the discharge from the wellbore is mon-itored throughout the time. The CHT is commonly adopted to per-form in low permeability aquifers[47]. Advantage of CHT over CRT is that the effect of wellbore storage at the test well is minimized

[47]. Note that the articles on CHT mentioned below are associated with the transient ﬂow because the drawdown solutions tend to be identical in both CHT and CRT at sufﬁciently large time[48].

The literature on the solutions of CHT performed in confined aquifers is given in the following. Jacob and Lohman[48] consid-ered a radial flow problem in a confined aquifer during the CHT and gave the solution of transient flux across the wellbore. Han-tush[2]presented a transient drawdown solution in response to a fully penetrating well extracting groundwater from confined aquifers. Cassiani et al.[6]proposed a semianalytical expression for a flowing partially penetrating well with infinitesimal skin sit-uated in an anisotropic confined aquifer. Other researches with re-gard to CHT conducted in confined aquifers can also be found in literatures[5,49,50]. For unconfined aquifers, for example, Chen and Chang[51]considered the skin effect and developed a well-hydraulic theory for CHT in unconfined aquifers. Chang et al.

[52,53]gave hydraulic head solutions for CHT performed at a par-tially penetrating well in unconfined aquifers. The review of the re-cent development in the literature shows that there are a limited number of researches on the CHT in leaky and multilayered aqui-fers. Hantush[54]gave the transient solutions of drawdown and flux from the wellbore for leaky aquifers. Wen et al.[55] consid-ered the effect of finite-thickness skin and presented a mathemat-ical model for radial groundwater flow to a pumping well in a leaky aquifer during CHT.

1.4.3. Slug test

A slug test is a particular type of aquifer test in which a small amount of water in the control well is instantaneously added/ removed to/from the well and the response in drawdown is mon-itored through time in the control well or the surrounding observa-tion wells. The slug test has become popular in groundwater investigations because of its logistical and economic advantages

(4)

over other aquifer tests. The advantages include that the cost of test is low, the procedure is relatively simple and rapid, little or no water needs to be handled during the test, and the test provides information on local hydraulic properties[56]. The disadvantages include that the aquifer properties obtained from this test repre-sent the formation only near the test wellbore and the result of the test may be inﬂuenced by gravel or sand pack material in the borehole adjacent to the well screen.

The slug test is usually performed in confined or unconfined aquifers for most engineering practices. A number of mathematical models for analyzing slug test data have been published. For con-fined aquifers, for example, Cooper et al.[57] assumed the well diameter is finite and provided an analytical solution for the anal-ysis of slug tests conducted in fully penetrating wells. Faust and Mercer[58]simulated four slug test runs with a well skin hydrau-lic conductivity of 0.01, 0.1, 1.0, and 10 times the hydrauhydrau-lic con-ductivity of the formation. According to their calculation, the skin with a much lower permeability than that of the surrounding formation can lead to an incorrect estimate of the true hydraulic conductivity. Yeh et al.[59]presented a semi-analytical solution for a slug test performed at a partially penetrating well in a con-fined aquifer, accounting for the effect of finite-thickness skin. Other researches on the slug test conducted in confined aquifers can also be found in literatures[60–62]. For unconfined aquifers, Bouwer and Rice [63] proposed a procedure for calculating the hydraulic conductivity of aquifers with partially or completely pe-netrating wells. Dagan[64]presented simple numerical schemes for determining the hydraulic conductivity of unconfined forma-tions from the analysis of recovery, packer, and slug test data. But-ler et al.[65]presented a procedure for analyzing data from slug tests performed at partially penetrating wells in highly permeable formations which could be confined or unconfined. For multilay-ered aquifers, Butler et al.[66] employed a numerical model to evaluate the capability of multilevel slug tests in providing infor-mation about vertical variation in hydraulic conductivity in multi-layered aquifers.

1.4.4. Recovery test

At the end of a pumping/injection test, the water levels in the test and observation wells begin to increase/decrease. The process involved in the change in water levels are referred to as recovery, and the observed drawdown remaining during the recovery peri-od is named as the residual drawdown. The equations describing the residual drawdown in a piezometer[67]and an observation well [67] during the recovery period after constantly pumping at a partially penetration well in a homogeneous conﬁned aquifer have also been developed. The analysis of residual drawdown data can estimate the transmissivity and provide an independent check on the hydrogeological parameters found from the previous drawdown data analysis. There are many studies associated with the data analyses of water level recovery after a pumping at a constant rate in conﬁned aquifers with neglecting well radius

[68,69] or considering well radius [70], leaky aquifers [71], unconﬁned aquifers [72], as well as for head recovery after a packer test in unsaturated fractured media [73]. Shapiro et al.

[74] presented a model based on the solutions of Papadopulos and Cooper[75] and Cooper et al. [57] to determine the early-time recovering water level following the shut down of a constant pumping in wells with turbulent head losses. Recently, Samani et al.[76]used derivative analysis of pumping and recovery test data to estimate the parameters of Shiraz aquifer in Fars province, Iran. Yeh and Wang[77]developed a mathematical model to de-scribe the residual drawdown with consideration of the wellbore storage effect and the drawdown distribution occurring at the end of a previous CHT.

2. Solution methodology

A partial differential equation for describing groundwater ﬂow can be solved either analytically or numerically. In this section, we summarize the methodologies and application of analytical, semi-analytical, numerical and approximate approaches which are often made in the groundwater area.

2.1. Analytical methods

As mentioned in Section 1, the phenomena of groundwater flow can be mathematically described by a partial differential equation (PDE) or a system of several PDEs, with associated boundary and initial conditions. Analytical methods or numerical techniques can be used to solve the flow equations for problems arisen in well hydraulics. In this section, we introduce some mathematical ap-proaches commonly used in solving the groundwater problems analytically, including Laplace transform, Hankel transform, Fou-rier transform, Mellin transform, Green’s function, dual integral/ series equation method, and Boltzmann transform. Basically, the first four methods can be considered as different types of integral transform, which offers an easy and effective way in solving a vari-ety of problems arising in engineering and physical science[78]. The concept of integral transform is somewhat analogous to that of logarithmic transform. Their main aim is to transform the given problem into one with reduced number of dimensions. By taking an integral transform, a PDE can be reduced to an ordinary differ-ential equation (ODE) in the transformed domain. The solution can then be obtained by solving the ODE along with the associated ini-tial or boundary conditions.

2.1.1. Laplace transform

Laplace transform F(s) of a function f(t) can reduce initial value problems with certain types of ODEs to the algebraic solutions. It can also be used to transform initial boundary value problems with certain classes of PDEs to ODEs[79]. Generally, Laplace transform is proved to be most useful in reducing the ‘‘time-like’’ variables. Note that a necessary condition for the existence of the Laplace transform is that the value est_{f vanishes when t ? 1, otherwise}

the integral defined as the Laplace transform does not converge. In well hydraulics, Laplace transform is a powerful method to deal with transient groundwater problems when a first-order time derivative of function, i.e., the drawdown or hydraulic head, is in-volved. After taking Laplace transform with respect to the time var-iable, the original governing equation describing the transient groundwater flow becomes a boundary value problem. Notice that the transform procedure needs the information of the function va-lue at the initial stage, i.e., the initial condition. Plenty of researches were conducted using Laplace transform for problems in the well hydraulic literature[55,60,80–83].

2.1.2. Fourier transform

The existence of the exponential Fourier transform F(w) of a function f(x) presumes that f(x) vanishes when |x| is infinity. Nor-mally, such condition at either end of an infinite interval occurs when x is a ‘‘space-like’’ variable[79]. Hence, the exponential Fou-rier transform turns out to be helpful in reducing the ‘‘space-like’’ variables and it can be used to solve the boundary value problems for PDEs which describe groundwater flow problems with infinite boundaries. On the other hand, the Fourier sine (or cosine) trans-forms are well suited for solving flow problem of a semi-infinite domain represented by a second-order PDE when the function va-lue f (or its derivative fx) is known at the origin and the function

and its derivatives are required to vanish as x ? 1[79]. Moreover, the ﬁnite Fourier sine and cosine transforms are particularly useful

(5)

in solving the problem of a finite domain. The Fourier transforms and finite Fourier transforms are appropriate in solving the problems involved the Dirichlet or Neumann boundary conditions. Generally, the governing equation describing the transient ground-water flow after taking the Fourier transforms reduces to an initial value problem. Applications of these transforms can be found, for example, in Tsou et al.[84]who used exponential Fourier, Fourier sine, and finite Fourier cosine transforms to solve a 3-D groundwa-ter flow problem for a confined aquifer with a slanted well near a stream, in Huang et al.[85]who applied exponential Fourier, Fou-rier sine, and Laplace transforms to obtain an analytical solution for describing the head distribution in an unconfined aquifer with a single pumping horizontal well parallel to a fully penetrating stream, and in other literatures[10,86]employing the finite Fou-rier sine or cosine or both finite FouFou-rier sine and cosine transforms to solve a variety of groundwater flow problems.

2.1.3. Hankel transform

The Hankel transform F(

q

) of a function f(r) has been found to be the most appropriate for problems formulated with equations expressed in polar and cylindrical coordinates. In applying the Hankel transform to such equations involving the radial variable ranging from zero to inﬁnity, it is required that the value of f and frare bounded at the origin and

ffiffiffi r p

f and pffiffiffirfr vanish at r ? 1

[87]. Yeh and Chang[10]presented analytical solutions, developed via finite Fourier sine transform and Hankel transform, for describ-ing hydraulic head due to pumpdescrib-ing in a wedge-shaped aquifer with an infinite radius. With the aid of Hankel transform and Laplace transform, Malama et al.[46]presented a semi-analytical solution for the transient streaming potential response of an unconfined aquifer to the continuous pumping at a constant rate.

2.1.4. Mellin transform

The Mellin transform is closely related to the Fourier transform. Morse and Feshbach[88]demonstrated an example on the relation between the Mellin transform and Fourier transform, and men-tioned that the theorems given for the Fourier transform theory can be rewritten and applied to Mellin transforms. Sneddon[89]

showed that the Mellin transform can be expressed in a natural manner in the solution of boundary value problems regarding an infinite wedge. Chan et al.[90]used Mellin transform to find the solution for steady-state drawdown due to a single well pumping at a constant rate from a wedge-shaped confined aquifer. Craster

[91] solved several free boundary conditions encountered in groundwater ﬂow problem by the use of Mellin transform.

2.1.5. Green’s function

Instead of infinite series or eigenfunction expansions, solutions of ODEs and PDEs can be expressed in terms of an integral with a fundamental function. The fundamental function, usually named as Green’s function, is a part of a closed-form solution for linear models and represents the solution as the result of the influences of a system to a unit point load or source. The total responses are the linear superposition throughout the system; hence, the total solution can be obtained by adding the responses of influences in all parts of the domain[92]. Indelman and Zlotnik[93]used the mean Green function to develop a mathematical model of average nonuniform flow in heterogeneous stratified media. Using pertur-bation methods coupled with Green’s function techniques, Axness and Carrera[94]presented a second-order analytical solution for the steady hydraulic head flow in a two-dimensional (2-D) con-fined heterogeneous aquifer. Asadi-Aghbolaghi and Seyyedian

[16]developed Green’s function solution for groundwater ﬂow to a vertical well in a triangle-shaped aquifer using image well theory and double Fourier series.

2.1.6. Dual integral/series equation

A boundary value problem under a mixed boundary condition can be converted to dual integral/series equations using integral or finite transforms. The resulting dual equations are in the form of integral or series depending on whether the problem is infinite or finite. Once the dual integral/series equations are solved, the solutions of the mixed boundary value problems are obtained. More detailed development of the theory of dual integral/series equations is provided in Sneddon [95]. Cassiani and Kabala[96]

used the method of dual integral equations to develop semi-ana-lytical solutions for describing the well response to the pumping test and slug test performed at a partially penetrating well in a confined aquifer of semi-infinite vertical extent. Their solution is applicable for aquifers of a semi-infinite vertical extent or the sit-uations where the bottom boundary of the aquifer is far enough from the tested area. Since the thickness of aquifer is generally fi-nite in real world, Chang and Yeh[97]used the method of dual ser-ies equations to obtain a drawdown solution to the CHT performed at a partially penetrating well in an aquifer with a finite thickness. Chang and Yeh[98]further extended the work of Chang and Yeh

[97]to more general problems that the screen of a partially pene-trating well is allowed to install at varied depth of the aquifer. The way they solved the problem was to employ the method of triple series equations, which bears the similar concept to the method of dual series equations.

2.1.7. Boltzmann transform

The Boltzmann transform is a special case of a general method named similarity method. One can transfer a nonlinear diffusion equation into an ordinary differential equation by introducing a similarity variable, i.e., x=pffiffit, which is a combination of indepen-dent variables. In the well hydraulics, Boltzmann transform is com-monly used in solving two kinds of problems, i.e., Boussinesq equation for unconfined aquifers and non-Darcian flow problem. In solving the Boussinesq equation of unconfined aquifer, Boltz-mann transform was employed, for example, in Wang and Zhan

[99]who presented a solution for transient confined–unconfined flow due to a well of infinitesimal radius at a constant rate pump-ing. Among the published works in using Boltzmann transform to solve non-Darcian flow problems we may mention as follows. Sen [100,101] solved non-Darcian radial flow with and without considering the effect of wellbore storage. Wen et al.[102]derived the solutions for non-Darcian flows in a single vertical fracture to a well on the basis of Izbash power-law and Wen et al.[103] inves-tigated non-Darcian flow toward a finite-diameter pumping well with and without considering the wellbore storage and derived the solutions of non-Darcian well functions at the wellbore.

2.2. Semi-analytical methods (Laplace inversion)

The transient groundwater flow equation, Eq.(1), is a diffusion type of PDEs with a first order differential in time. In most ground-water problems, Laplace transform is suitable to apply to Eq.(1)for removing the time variable t and transform PDEs into ODEs in La-place domain. The LaLa-place-domain solutions may be obtained after solving the resulting ODEs with associated Laplace-domain bound-ary conditions. To acquire the inverse Laplace transform, one may use tables[104,105]together with rules or methods such as the shift theorem, partial fractions, and convolution theorem. The time-domain solutions (i.e., exact solutions or analytical solutions) can also be obtained by complex variable theory referred to as complex inversion integral or Bromwich’s integral[106], occasion-ally along with the residue theorem and/or Jordan’s lemma. How-ever, the inversion of Laplace transform is generally rather complicated and may not be tractable in some groundwater flow problems. Moreover, the time-domain solutions are usually in

(6)

terms of improper integrals with limits of integration from zero to inﬁnity and theirs integrands comprise singularity at the origin

[27]. The integrands of those time-domain solutions are in terms of oscillatory functions comprised of terms of the product of the Bessel functions of the first and second kinds of zero and first or-ders. The numerical calculations for those solutions are therefore time-consuming and very difficult to perform accurately. There-fore, numerous methods have been devoted to the numerical inversion of the Laplace transform. For comprehensive bibliogra-phies, the reader may refer to[107,108]. In the following we pres-ent a brief introduction to three inversion methods, namely Stehfest, Crump, and modified Crump methods, which are the most widely used approaches for numerical Laplace inversion in well hydraulic problems. Detailed algorithms of each can be found in Cheng[109]and de Hoog et al.[110].

2.2.1. Stehfest method

Since the complex analysis in Laplace inversion is difﬁcult to implement, Stehfest[111]proposed a series to approximate the La-place transform of a real-valued function f(t) using the following formula f ðtÞ ln 2 t PN n¼1 cnF n ln 2 t ð9Þ where cn¼ ð1ÞnþN=2 P minðn;N=2Þ m¼ðnþ1Þ=2 mN=2_ð2mÞ! ðN=2 mÞ!m!ðm 1Þ!ðn mÞ!ð2m nÞ! ð10Þ

The numbers of terms N in the series must be even. According to Stehfest[111], the accuracy can be improved by increasing N terms. However, round off error limits the value of N. The value of N is sug-gested to be in the range 6 6 N 6 20 for most of engineering pur-pose[109].

2.2.2. Crump method

Based on the trapezoidal rule, Crump[112]approximated the Laplace inversion by the following equation

f ðtÞ e ut Tp _Fð

u

Þ 2 þ P1 k¼1 Re F

u

þk

p

i Tp cos k

p

t Tp Im F

u

þk

p

i Tp sin k

p

t Tp ð11Þ A parameter / is introduced as /¼

u

ln E 2Tp ð12Þ

where E is the error tolerance,

u

is the real part of the leading pole of the function F(p), and Tpis a periodic function to approximate f(t).

The default value of

u

equals zero for functions without poles. Crump[112]used the epsilon-algorithm of Wynn[113](also called the Shanks method) to accelerate the convergence of the sequence in Eq.(9). The Shanks method developed by Shanks[114]is a non-linear sequence-to-sequence transformation. This transform is effective in accelerating the convergence of slowly convergent se-quences and inducing convergence in divergent sese-quences.

2.2.3. Modiﬁed Crump method

A signiﬁcant improvement over the Crump method is developed by de Hoog et al.[110]. They applied the Pade approximations to improve the acceleration procedure for approximating the trans-form of the sequence in Eq.(9). The Pade approximation can be ex-pressed by the quotient of two polynomials as[115]

f ðtÞ ﬃ RNðtÞ ¼

a0þ a1t þ a2t2þ þ antn

1 þ b1t þ b2t2þ þ bmtm

; N ¼ n þ m ð13Þ

where the parameters a0, a1, . . . , anand b1, b2, . . . , bm, are available for

the approximation of f(t), by RN(t). This method is designed to

re-duce truncation error that always occurs due to the fact that the Fourier series is not an inﬁnite series. A useful Fortran routine INLAP of IMSL[116]developed according to the work of de Hoog et al.

[110]is available for Laplace inversion.

It has been realized by researchers that it is impossible to devise a universal algorithm that performs accurately for all types of func-tions. Neither Stehfest method nor Crump method can perfectly implement in problems of well hydraulics. Each of them has their merit or defect in some particular problems. The execution time is less in Stehfest method, while the result from Crump might be more accurate if more terms are given. It can be easily found from the literatures that used Stehfest method [11–13,34,36,53,117], Crump method[118], and modiﬁed Crump method[59,115,119]

to numerically inverse the Laplace-domain solution into time-domain.

2.3. Numerical methods

In many practical problems of groundwater flow, analytical solutions are not possible due to the facts that the shapes of the boundaries might be irregular, the coefficients appearing in the governing equations and in the boundaries conditions might be space-variant, the dependent variables in the initial conditions might be non-uniformly distributed, and the source/sink term might be in nonanalytic forms. To deal with realistic situations, numerical techniques provide convenient, flexible, and powerful tools for solving groundwater flow problems in complex field situ-ations as discussed above. Four widely used numerical methods in well hydraulics are the finite difference, finite element, boundary element and analytic element methods addressed below.

2.3.1. Finite difference method

The finite difference method is probably the first numerical ap-proach used to solve PDEs. A PDE describing the head or drawdown distribution can be approximated by a single difference equation over a time interval or a system of algebraic difference equations over one or two time intervals at predetermined, finite number of discrete grid nodes in the problem domain. Forward, backward, and central difference schemes are commonly used for the finite difference approximation. The monographs of applications of finite difference methods to the problems in well hydraulics have been published and the interested reader may be referred to the book by Wang and Anderson[120]or Remson et al.[121]for the techni-cal details.

2.3.2. Finite element method

Another powerful numerical technique, known as the finite ele-ment method, has been widely applied to numerous groundwater flow problems as well. While the finite difference method is usu-ally implemented with rectangular cells, the finite element method is regarded as a flexible approach which can handle almost any shape of flow boundary and any combination of boundary condi-tions, inhomogeneous and anisotropic media, moving boundaries, free surfaces and interfaces, and multiphase flows[122]. For a con-cise yet comprehensive discussion of the application of finite ele-ment method in groundwater hydrology, the reader may be referred to Huyakorn and Pinder[123]or Lee[124]. The book by Yeh [125] is an advanced treatise on computational subsurface hydrology.

(7)

2.3.3. Boundary element method

The boundary element method or the boundary integral equa-tion method consists of formulating boundary value problems in terms of an integral equation. The solution of the integral equation can be obtained by approximating the boundary as a series of straight lines or elementary curves with the simpliﬁed assump-tions to the behavior of the solution along the boundary segments

[126]. In other words, the solution exactly satisﬁes the governing differential equation, but approximately meets the boundary con-ditions. Application of the boundary integral equation method in groundwater ﬂow is discussed in detail by Liggett and Liu[127].

2.3.4. Analytic element method

The analytic element method is usually applied to the aquifers of an inﬁnite extent and uses superposition to generate the analyt-ical solution of a certain problem, expressed as a sum of basic solu-tions, each with a number of possibly unknown parameters. These parameters are then determined from the boundary conditions. One of the primary distinctions between analytic and boundary element methods is that the boundary elements are always given in the form of integrals, which are often calculated numerically in boundary element method and analytically in analytic element method. The book associated with the analytic element method by Strack [126] provided a basic theoretical framework for the understanding of the analytic element method and the detailed mathematical descriptions of the analytic elements and their numerical implementation.

2.4. Approximate methods

The analytical solutions obtained from groundwater flow prob-lems are often in the form of infinite series for aquifers of a finite domain[4]or integrals for aquifers of an infinite domain [128]. Those solutions in terms of infinite series may converge slowly and are not suitable for numerical computation for very small val-ues of time[129]. On the other hand, the solutions are in the form of an integral mostly with the integration limits from zero to infin-ity and with the integration variable in the denominator of the integrand, posing the problem of singularity at the origin of the integration[50]. Due to the presence of singular point, the numer-ical calculation of those solutions is generally very time-consuming and rather difficult to achieve accurate results. Therefore, for solv-ing the problems in well hydraulics there is a need to develop approximate solutions which have simpler forms than the analyt-ical solutions and are much easier in describing the transient behavior of groundwater flow with desired accuracy. Without the aid of the numerical techniques such as finite difference methods and finite element methods, the approaches in the development of approximate solutions may be divided to three different catego-ries. The first is to solve the flow equation slightly different from the original one for finding the approximate solution. The second is to derive the approximate solution based on the Laplace domain solution along with the relationship of large p (Laplace variable) versus small t (hereinafter referred to as LPST) or small p versus large t (hereinafter referred to as SPLT). The last is to apply various types of approximate techniques into the time domain solution.

2.4.1. Approximation in governing equation

In the first category, the perturbation method is commonly used to find an approximate solution to nonlinear equations with vari-able coefficients or an irregular domain which cannot be solved exactly. By adding a ‘‘small’’ parameter to the mathematical description of the problem, perturbation method decomposes a tough problem into an infinite number of relatively easy ones. Hence, the perturbation method is most useful when the few first-order perturbation terms reveal the important features of

the solution and the remaining ones give small corrections of the approximation. The researchers who used the perturbation meth-od to solve the problems in well hydraulics are, for example, Batu and van Genuchten[130]and Moutsopoulos and Tsihrintzis[131]. The former adopted a singular perturbation method to solve the Boussinesq equation for the problem of a constant injection into a radial aquifer while the latter derived approximate analytical solutions for transient state, nonlinear ﬂows through porous media based on the perturbation method. The other approach used to de-velop an approximate solution of the diffusion equation was pre-sented by Fang et al. [132] for a problem in electrochemical. They extended the solution of the diffusion equation of steady-state process to the non-steady-steady-state solution. Perrochet [133]

used a similar concept to develop an approximation solution, which is exactly the same as that given in Fang et al.[132], for tran-sient wellbore ﬂux to a well subject to a constant drawdown. The approximate solution obtained from this approach is generally applicable for all values of the time.

2.4.2. Approximation in Laplace domain solution

The Laplace domain solutions of groundwater flow equation for describing flow in unbounded aquifers may be in a form with a quotient of the modified Bessel functions of the second kind of or-der zero or one[28]. For example, the Laplace domain solution for a CHT in aquifers of an infinite extent can be expressed, in our nota-tion, as h(r, p) = hwK0(

g

r)/[pK0(

g

rw)][49]in which rwis the well

ra-dius, p is the Laplace variable, hw is the hydraulic head at the

wellbore,

g

¼pffiffiffiffiffiffiffiffiffip=D, D is the hydraulic diffusivity, and K0() is the

modiﬁed Bessel function of the second kind of order zero. Gener-ally speaking, the approximate time domain solutions can be ob-tained by ﬁrst using series expansion of the Bessel functions and then inverting the quotients associated with polynomials in p based on the methods of partial fractions and table of Laplace transforms given in, for example, Spiegel[106]. Carslaw and Jaeger

[129]gave a systematic discussion of the methods in ﬁnding solu-tions useful for small or large values of the time. They mentioned that the method of Heaviside’s series expansion can be used to ex-pand the Laplace domain solution denoted as xðpÞ in a series of power of (1/p), and with this the corresponding time domain solu-tion x(t) can then be obtained as a power series in t. Addisolu-tionally, they also pointed out that some Laplace domain solutions after taking the asymptotic expansions may result in a series of negative exponentials and these results generally lead to the solutions use-ful for relatively small values of the time. Carslaw and Jaeger[134]

presented the analytical solutions to describe the temperature dis-tributions for a wide variety of heat flow problems. Their solutions can be applied to groundwater flow problems based on physical analogy. For the CHT, they obtained a drawdown solution for small values of time and both small-time and large-time solutions of wellbore flux as well by using asymptotic expansions of the Bessel functions. Furthermore, they also gave a large-time solution for drawdown distribution due to CRT at a pumping well. van Everdin-gen and Hurst[135]introduced the concept of symbolic relation between the derivative operator of time, i.e., d/dt, and p and in-ferred that p must be small if t is large, or inversely, p should be large if t is small. Accordingly, one might obtain a small time or a large time solution based on the LPST or SPLT relationship, respec-tively, to the Laplace domain solution. Yeh and Wang[136] pre-sented a short review on the application of LPST[137], SPLT[38], and both [138,139]to groundwater flow problems for obtaining the approximate solutions. In addition to those articles, a few works have been carried out to develop the approximate solutions in the groundwater literature. For example, among the published works on applying the relationships of LPST and SPLT we may men-tion the following. Chen and Chang[51]developed the early- and late-time solutions for a CHT in a homogeneous and anisotropic

(8)

unconfined aquifer with the skin effect. Hunt and Scott[140] pre-sented two approximate solutions for flow to a well in a leaky con-fined aquifer overlain by two layers; one is an unconcon-fined aquifer on the top and the other is an aquitard in between. Pasandi et al.

[34]also provided early time and late time drawdown solutions in addition to the Laplace domain solution for CRTs at a partially penetrating well in a phreatic aquifer with considering the effects of wellbore storage and ﬁnite thickness skin. Tsai and Yeh[141]

developed a large time solution for the wellbore flow rate induced by a constant-head test in two-zone finite confined aquifers by employing the relationship of SPLT and L’Hospital’s rule. It is of interest to mention that some approximate solutions had been developed in other areas such as the solute transport problems based on the relationships of LPST[142]and SPLT[143]and the dual-porosity media problem using the relationships of both LPST and SPLT[144]. Note that care must be taken when applying the relationship of SPLT to derive the large time approximation to the Laplace domain solution. The large time solution of the well-bore flux for the constant-head test in a homogenous confined aquifer can be written as Q ðrw;pÞ 1=½pInðp=kÞ [49] where

k¼ 4T=ðe2c_r2

wSÞ and

c

= 0.57722. . . is Euler’s constant. This

approx-imate solution introduces a pole at Re p = k which is artiﬁcial and must be excluded[136,145].

2.4.3. Approximation in time domain solution

A large amount of methods had been proposed to develop the approximation of the time domain solutions for various problems in well hydraulics. Tseng and Lee[146]reviewed six approximate methods used to estimate Theis well function and gave a concise introduction to each method. Those six methods are series expan-sion, asymptotic expanexpan-sion, continued fraction, polynomial, ra-tional and Chebyshev polynomial approximations, recurrence relations, numerical quadrature and inverse Laplace transform techniques. The application of each of the first four methods to a time domain solution yields an explicit mathematical expression, while that of the last two methods will give a procedure for the numerical computation of the time domain solution. Among those four methods, the first two, series expansion and asymptotic expansion may be the most widely used techniques to find approx-imate solutions in groundwater problems. The Cooper–Jacob equa-tion [147] is a good example of the application of the series expansion. This equation results from truncating the high order terms in the infinite series obtained after applying the series expansion to the Theis solution[1]. Some studies can be found in the well hydraulic literature that adopts the method of series expansion [148] or asymptotic expansion [149] to develop the approximate solutions. Interestingly, Barry et al.[150]presented an analytical expression for the Theis well function based on the interpolation between the series expansion of small u and the asymptotic expansion of large u.

3. Subjects in modeling groundwater ﬂow

Due to the increasing needs of irrigation, industrial, urban and suburban expansion, various subjects concerning groundwater flow attract scientists and engineers’ attentions. In this section, we first address topics of current interests regarding the effects of finite well radius, wellbore storage, well partial penetration, and the presence of skin zone on the groundwater flow problems. We then introduce the recent developments concerning the flow in aquifers with horizontal wells or collector wells, the capture zone delineation, and the non-Darcian flow in porous media and frac-tured formations. Finally, we give a critical synthesis of the body of work on numerical computations for five well functions often presented in the well-hydraulic literature.

3.1. Effects of ﬁnite well radius and wellbore storage

For small-diameter wells with radius varies between 0.05 m and 0.25 m, the groundwater flow is often modeled by neglecting the effect of well radius. In reality, large-diameter wells with ra-dius ranges from 0.5 m to 2 m are commonly installed in many countries to meet a large demand for domestic and irrigation water uses. The behavior of drawdown solution in wells with a finite well radius is different from that in wells with an infinitesimal radius because of the effect of well radius and the contribution of well-bore storage. Papadopulos and Cooper [75] presented an exact solution for the drawdown due to pumping in a large-diameter well. Their solution took into account the effects of finite well ra-dius and the water stored inside the wellbore, which were ne-glected in the Theis equation. Theoretical models are developed for incorporating the effect of finite well radius[151,152]and both effects of finite well radius and wellbore storage into models

[34,153]. Other issues of large-diameter wells such as seepage face, an application of discrete kernel approach and well function approximations are widely discussed in the literatures. For exam-ple, Ojha and Gopal[154]proposed a model based on flow velocity to describe the seepage face variation for large-diameter wells in an unconfined aquifer. Sakthivadivel and Rushton[155]considered the dynamic seepage face and gave a methodology to determinate the aquifer parameters from pumping tests with large-diameter wells. Mishra and Chachadi[156]extended the discrete kernel the-ory for analyzing flow behavior toward a large-diameter well dur-ing the recovery phase. Chachadi and Mishra[157]used discrete kernel approach to analyze the unsteady flow toward a large-diameter well in a confined aquifer while determining well loss component. Since the well function for large-diameter wells pro-posed by Papadopulos and Cooper[75]is an integral function of transmissivity and storativity, it is difficult to explicitly determine those two aquifer parameters from pumping test data. The devel-opment of well function approximation for large-diameter wells has been the subject of many studies[158].

3.2. Skin effect

The drawdown induced by pumping conducted in aquifers may be influenced by a region near the well referred to as the skin zone, or simply the skin. This zone has a lower or higher permeability than the adjacent formation. The skin is introduced during the well installation process including the drilling, construction, and instal-lation of annular fill of the pumping well. Two types of the well skins are classified by their permeability. If the permeability of the skin zone is less than that of the formation zone, this kind of skin will be called a positive skin. On the contrary, if the permeabil-ity of the skin zone is larger than that of the formation zone, the skin will be a negative skin. The approaches of quantifying the skin effect in the well hydraulic modeling by researchers may be grouped into two ways; one assumes the skin to be infinitesimally thin, while the other considers it to be of finite thickness.

3.2.1. Inﬁnitesimal skin

Studies of infinitesimal skin have been widely investigated in the areas of petroleum industry and well hydraulics. The skin effect is regarded as an additional pressure drop near the wellbore and proportional to the wellbore flow rate. Under this consideration, both the skin thickness and storativity are neglected in the skin zone. As such, a skin factor is introduced to represent the lumped properties of the skin and used as an energy loss term at the well-bore for mathematical convenience. The published works assum-ing the skin to be infinitesimal are, for example, Cassiani et al.

(9)

3.2.2. Finite thickness of skin

Novakowski [161] mentioned that the thickness of the skin zone may range from a few millimeters to several meters. The ap-proaches of developing the solutions by taking into account the skin thickness can be further divided into two categories. The ﬁrst category assumes that the skin thickness is a parameter associated with the skin factor. The other treats the skin as another formation, and therefore one should handle a two-zone ﬂow system for such an aquifer.

In an aquifer system subject to a pumping at a constant rate, Hawkins[162]took into account the skin thickness and developed a formula for the additional pressure drop to quantify its effect. Based on the pressure drop equation proposed by van Everdingen

[163], Hawkins deﬁned the formula of skin factor as

sf ¼ ½ðK=KskinÞ 1 lnðrskin=rwÞ ð14Þ

where rskindenotes the skin radius, K and Kskinare the hydraulic

con-ductivities of the formation and skin, respectively. The productivity ratio of a well deﬁned by Hawkin is In(re/rw)/[In(re/rw) + sf] with

drainage radius re. Thus the well productivity ratio becomes very

small when Kskintends zero and it depends only on drainage,

well-bore and skin radius when Kskintends to be a large value. Hawkin’s

skin factor has been adopted in CHTs, CRTs and slug tests by, for example, Faust and Mercer[58]and Chen and Chang[164].

Moench and Hsieh[165]considered an equation for describing the ﬂow in skin zone with negligible storativity and provided a hydraulic head solution for the slug test. They deﬁned the skin fac-tor as

sf ¼ ½ðK=KskinÞ lnðrskin=rwÞ ð15Þ

For CRTs, Moench[160]considered the ﬂow toward a partially penetrating well in a slightly compressible water table aquifer and assumed that the storage capacity of the skin is negligible. He assumed that the ﬂow rate through the skin zone is equal to that along the wellbore and gave another skin factor as

sf ¼ ðKds=KskinrwÞ ð16Þ

where dsis the skin thickness.

A ﬁnite skin region, however, has its storage capacity in nature. The newly introduced mathematical model is thus required for the aquifer system representing by two regions of radial ﬂow and each of it with individual transmissivity, storativity, and thickness. A two-zone aquifer may be named as a patchy aquifer [166,167]

and such an aquifer system may be referred to as a composite sys-tem[161]or a two-layered system[27]. By also considering the skin storativity, Moench and Hsieh[168] further extended their solution in Moehch and Hsieh[165] to solve the ﬂow equations in both skin and formation zones and developed the head distribu-tions for slug tests. Many other recent studies have modeled the two-zone aquifer system for CHTs[5,50,167], for CRTs[34,37,38, 161,169,170], and for slug tests[167,171]. In addition, Perina and Lee[153]further derived a generalized solution which is applicable to CHT, CRT, and slug test in the two-zone aquifer system.

3.3. Effect of well partial penetration

If the length of a well screen is less than the entire thickness of an aquifer, the well is called a partially penetrating well[172]. The well of partial penetration is often adapted for hydrogeological site investigation to characterize groundwater pollution problems. For example, the contaminant of light nonaqueous phase liquids gen-erally forms a pool on the top of water table in the subsurface. Under this circumstance, the screen of the well is only open at the upper part of aquifer. On the other hand, the well is usually screened at the lower part of aquifer if the contaminant is dense

nonaqueous phase liquid. The hydraulic head distribution in re-sponse to an aquifer test at a partially penetrating well will be influenced by the screen length and its location in the aquifer. The equation describing the groundwater flow should include the vertical flow component induced by pumping at a partially pene-trating well. In this section, mathematical approaches for dealing with the effect of well partial penetration in both CHT and CRT are discussed.

3.3.1. Effect of partial penetration in constant rate test

The theory of a partially penetrating well pumped at a constant rate was ﬁrst published in 1957 by Hantush for leaky aquifers

[172]. In his solution, the discharge is assumed to be uniformly dis-tributed over the well screen and the flow in the zone below the screen is neglected. The solution for leaky flow can be further sim-plified to the case in a non-leaky confined flow by making the hydraulic conductivity of the confining layer which overlays the confined aquifer equals zero. For the case that the screen does not install from the top of the aquifer, the solution can be found in Hantush[2]. The assumption of uniform discharge can also be applied to the case of unconfined aquifers. Neuman[173] devel-oped a solution of drawdown distribution in an unconfined aquifer by considering the effect of well partial penetration on drawdown during CRT. Yang et al. [28] took into account the radius of the pumping well and derived an analytical solution of drawdown in a confined aquifer with a partially penetrating well for the CRT. More references on the issues of partially penetrating effect during CRT can be found in the well hydraulic literature[34,39,44].

3.3.2. Effect of partial penetration in constant head test

In the well hydraulic modeling, the well water lever is main-tained constant at a CHT and the constant head is appropriate to specify at the screen section. If the aquifer is fully penetrated, the mathematical model describing the groundwater ﬂow can gener-ally be solved by the conventional integral transform techniques

[2]. If the aquifer is partially penetrated and the effect of well skin is negligible, the constant head condition will be used for the screen part and the no flow (or zero flux) will be specified at the casing. The flow problem induced by the partially penetrating well therefore becomes a mixed boundary value problem[98]. Some analytical approaches had been used to deal with the mixed boundary value problems happen in CHTs in the well hydraulic lit-erature. Those approaches may be classified into following five categories:

The ﬁrst approach is to replace the constant head along the well screen by a uniform radial ﬂux [151] which will transform the mixed boundary into a homogeneous Newman boundary and re-sults in an approximate solution.

The second approach is the well screen discretization method which was adopted by Chang and Chen[174]and Perina and Lee

[153]. Along the screen, the flux is assumed constant in the first ap-proach and considered to be non-uniform in the second apap-proach. The third is the domain decomposition method which divides the flow domain into two or several sub-regions and each region has its own flow equation and associated initial and boundary con-ditions [5,52,53]. The flow equations are then solved simulta-neously via the continuity conditions for the hydraulic head and flow rate at the interfaces of the sub-regions.

The fourth is to use the method of dual integral equation and di-rectly solve the mixed boundary value problems (MBVPs) for aqui-fers of inﬁnite vertical extent and the screen installed at the top of the aquifer[6,96]. For aquifers of ﬁnite vertical extent, Chang and Yeh[97]solved the MBVPs by the method of dual series equations in conjunction with perturbation method.

For cases of screen is not placed at the upper part of the aquifer, the previous approach becomes inapplicable. The use of triple

(10)

inte-gral equation method, i.e., the last approach, can solve the prob-lems with a well screen arbitrarily located at any portion of the aquifer of ﬁnite vertical extent[98].

3.4. Unconﬁned ﬂow problem

The modeling unconfined flow may be classified into five differ-ent approaches. The first approach is to use the confined flow equa-tion to model the unconfined flow problems. The second approach neglects the vertical flow and uses the Boussinesq equation to rep-resent the unconfined flow. The third is based on the radial con-fined flow equation embedded with a delayed yield term. The fourth employs Eq.(1)and adopts a free surface equation to repre-sent the top boundary condition. The last is to solve the unconfined flow equation by accounting for unsaturated flow above the water table.

3.4.1. Conﬁned ﬂow equation

In some field or practical problems, the changes in water table due to pumping or recharge is very small and thus the confined flow model becomes applicable to simulate the water table or drawdown distribution[45].

3.4.2. Boussinesq equation

Based on the approach of mass balance and the Dupuit assump-tion of neglecting the vertical flow[122], the Boussinesq equation, which describes the 2-D unconfined flow, can be developed and written as @ @x Kxh @h @x þ @ @y Kyh @h @y ¼ Sy @h @t R ð17Þ

where Syis the speciﬁc yield. Eq.(17)is nonlinear because it

con-tains the products of h and oh/ox. Therefore, the analytical solution for the Boussinesq equation is rather difficult to obtain. In the past, this equation has been employed in modeling and analyzing head distributions in many aspects of unconfined aquifer flow. Bear

[122] mentioned two methods of linearization to facilitate the development for the solution of Eq.(17). The first method is to re-place the variable thickness, h, in Eq.(17)with an average aquifer thickness, b, if the change of h due to pumping or recharge is very small compared with the saturated thickness. Then Eq.(17) be-comes a linear equation with the same form as the confined flow equation when replacing Kxh and Kyh by Tx and Ty, respectively.

The second method is to rewrite the ﬁrst term on the right-hand side (RHS) of Eq.(17)as (Sy/b)o(h2/2)/ot. Eq.(17)then becomes a

lin-ear equation in h2_{. The published works associated with the ﬁrst}

method are, for example, Verhoest and Troch [175], Bansal and Das[176], Parlange et al.[177], Li et al.[178]and with the second method, for example, Hantush [179], Yeh and Chang[180], Ilias et al.[181]. Different from those two methods given in Bear[122], Pulido-Velazquez et al.[182]also provided a linearization method to solve the Boussinesq equation based on the technique of change of variables.

In the past, some studies were also devoted to developing the analytical solutions from non-linear steady-state or transient Boussinesq equation without using the linearization method men-tioned above[183–187]. The literature on steady-state solutions for the Boussinesq equation is, for example, given in the following. Basha and Maalouf[188]discussed two particular cases of the ana-lytical solutions for surface and groundwater ﬂows and compared these two nonlinear cases with the solutions derived from the lin-earized Boussinesq equation using the Green’s function solution. Based on the Boussinesq equation, Batu and van Genuchten[130]

utilized a singular perturbation method to solve the transient ﬂow induced by a constant injection into a radial aquifer. Lockington

et al.[189] applied similarity transform to solve the problem of unconfined flow with a boundary condition in terms of a power of time. Semi-analytical approaches based on the Boltzmann trans-form have also been developed to solve problems involved the Boussinesq equation. For example, Guo [190] solved the Bous-sinesq equation for transient groundwater flow between a reser-voir and an unconfined aquifer of semi-infinite extent using the Boltzmann transform and Newton–Raphson method. Wang and Zhan[99]also developed a solution based on the Boltzmann trans-form along with the Runge–Kutta method for transient confined– unconfined flow driven by a pumping well of infinitesimal radius with a constant rate of pumping. Their paper also presented a short literature review on the issue of confined–unconfined conversion due to pumping using analytical and numerical approaches.

The Boussinesq equation is often used to model the hydrologi-cal problems such as the flow due to surface recharge and/or evapotranspiration in horizontal aquifer [179,181]or in sloping aquifers [175,188], hillslope flow [191,192], coastal aquifer hydraulics[193,194], designing drainage systems[195], the flow caused by a flood event[196]; and the hydraulics of the stream-aquifer interaction [176,197]. In the review paper by Winter

[198], he addressed the studies of the interaction of groundwater and surface water focusing on different landscapes from alpine to coastal.

Since Eq.(17)is nonlinear, numerical approaches such as finite-difference methods, finite-element methods, or hybrid numerical method are inevitable for solving the equation. Anderson and Woessner[199]mentioned that a few investigators in the 1970s had solved Eq. (17) with numerical techniques specifically de-signed to handle the problem of nonlinearity. Kim et al. [200]

developed a transient, mixed analytical and finite difference mod-els to simulate hydrological behavior for investigating the patterns of infiltration, evaporation, recharge and lateral flow across hill-slope. Stagnitti et al.[201]developed an explicit finite difference scheme for the solution of the nonlinear Boussinesq equation to examine the profiles of the water table height in a shallow sloping aquifer.

Taigbenu[202]developed a simplified Galerkin’s finite element model to solve the nonlinear Boussinesq equation with the advan-tages of computing efficiency and requiring less computer storage. Upadhyaya and Chauhan[203]presented a finite element solution of nonlinear Boussinesq equation to compare his analytical solu-tion derived from the linearized equasolu-tion for the problems of water table variation in a sloping aquifer caused by the sudden rise or fall of the water level in a nearby stream. By coupling the Boussinesq equation with the Richards equation, Hilberts et al.[192] devel-oped a tetrahedral finite element model to investigate the role of unsaturated storage in the relationship between rainfall and re-charge. Tang and Alshawabkeh[204]proposed a semi-analytical time integration approach for the simulations of 2-D transient unconfined flow described by the nonlinear Boussinesq equation.

3.4.3. Confined flow equation embedded with a delayed yield term Boulton[30]extended the theory of transient confined flow to a pumped well and added the term, the rate of delayed yield, to ac-count for slow drainage from the unsaturated zone. This term was expressed as

l

S0Rt

0elðtt

0_Þ

@s=@t0_dt0

where

l

is an empirical con-stant; S0 _{is the delayed yield per unit area, per unit drawdown;}

and s is the aquifer drawdown. He solved the drawdown equation by using Laplace transform method and the Faltung theorem. Later, he presented two approximated solutions based on his previous delayed-yield solution for unconﬁned aquifer ﬂow [205]; one is for the case in which

g

= (S + S0_{)/S tends to inﬁnity and the other}

is for the case that the pumping time is small. In addition, he also provided delayed yield type curves in the paper for the analyses of pumping test data. Notice that the mathematical model of using