Water Levels Near a Well 

Discharging From an 

Unconfined Aquifer

GEOLOGICAL SURVEY WATER-SUPPLY PAPER 1536-B
Water Levels Near a Well 

Discharging From an 

Unconfined Aquifer

By IRWIN REMSON, S. S. McNEARY, and J. R. RANDOLPH

GROUND-WATER HYDRAULICS

GEOLOGICAL SURVEY WATER-SUPPLY PAPER 1S36-B

A numerical solution of the differential 

equation for unsteady flow in an 

unconfined aquifer

UNITED STATES GOVERNMENT PRINTING OFFICE, WASHINGTON : 1961
UNITED STATES DEPARTMENT OF THE INTERIOR 

STEWART L. UDALL, Secretary

GEOLOGICAL SURVEY 

Thomas B. Nolan, Director

For sale by the Superintendent of Documents, U.S. Government Printing Office 

Washington 25, D.G. - Price 150 (paper cover)
CONTENTS

Page'

Abstract________________________________________________________ 27

Introduction._____________________________________________________ 27

Statement of the problem_________________________________________ 28

Restatement of the differential equation______________________________ 29

Derivation of the iterative equations_____--------_--------_----______ 30

Conversion to an ordinary differential equation__________________ 30

Finite-difference form of the differential equation __________________ 32

Evaluation of I close to 6=6n _____________________________ 33

Application to a hypothetical problem____________________________ 36

References_ _____________________________________________________ 39

ILLUSTRATIONS

FIGURE 9. Relation between 0 and $ for different values of Dn ---------- 37

10. Enlarged portion of curves relating 6 and <j) for different values

of £> __________________________________ 38

TABLES

TABLE 1. First iteration___________________________________________ 36

2. Final iteration___________________________________________ 38

m

GROUND-WATER HYDRAULICS

WATER LEVELS NEAR A WELL DISCHARGING FROM 

AN UNCONFINED AQUIFER

By IKWIN REMSON, S. S. MCNEARY, and J. R. RANDOLPH

ABSTRACT

The differential equation describing unsteady radial flow to a well in an unconfined aquifer has no simple solution in terms of elementary functions because the 

transmissibility decreases as the aquifer is dewatered. A numerical solution and 

type curves are presented showing the head of water in a homogeneous aquifer 

as a function of radial distance from the well and time elapsed after the water 

level in the well is lowered to and maintained at a given level. This solution is 

similar to one presented by J. R. Philip in 1955 for linear diffusion in a flow 

system in which diffusivity is concentration dependent.

INTRODUCTION

The U.S. Geological Survey carried on studies of the occurrence 

and movement of subsurface water, particularly in the zone of aeration, at Seabrook, N.J., from 1950 to 1960. Because of similarities 

between unsaturated flow and unconfined saturated flow, these studies 

led to a numerical solution of the differential equation describing 

unsteady radial flow in an unconfined aquifer.

The theory and the solution are developed in this paper. The 

methods are demonstrated by applying them to a convenient hypothetical problem. Type curves are presented and are applied to 

practical problems in a report in preparation.

Acknowledgment is made of the generous collaboration of the 

Seabrook Farms Co. These studies were under direction of Henry C. 

Barksdale, area chief, Ground Water Branch, U.S. Geological Survey, 

Arlington, Va.

27
28 GKOUND-WATEK HYDRAULICS

STATEMENT OF THE PROBLEM

The differential equation governing unsteady radial flow to a well 

in an aquifer (Jacob, 1950, p. 366) may be written

M_An^~LXM m

bt~dr[_S()r]+rSdr' ( }

where

h=head of water at any point in the aquifer, in feet; 

£=time, in days, since pumping started; 

r=radial distance, in feet, from the discharging well; 

T= coefficient of transmissibility of the aquifer, or rate of flow 

of water at the prevailing water temperature, through a 

vertical strip of aquifer 1 foot wide extending the full 

saturated height under a hydraulic gradient of 100 percent. 

T is normally expressed in gallons per day per foot of 

aquifer width (Theis, 1935). However, it is more convenient in the present development to express it in cubic 

feet of water per day per foot of aquifer width; 

S= coefficient of storage of the aquifer, or volume of water 

released from or taken into storage by an aquifer per unit 

surface area of the aquifer per unit change in the component 

of head normal to that surface. It is a dimensionless 

quantity (Theis, 1935).

In the following application of equation 1, it is assumed that the 

aquifer is infinite in areal extent, homogeneous, and iso tropic; that the 

discharging well penetrates and receives water from the entire thickness of the aquifer; and that the well has a negligibly small diameter. 

In a confined aquifer which is not dewatered during the discharging 

period, the transmissibility may be assumed to be constant at all times 

and at all places. In such a confined aquifer, the water is released 

from or taken into storage by virtue of the elastic response of the 

aquifer and the water to pressure changes induced by the pumping. 

It is assumed that the water is removed from storage instantaneously as 

head declines. Solutions of equation 1 are available for such confined aquifers, where T and S may be taken as constant (Theis, 1935; 

Jacob, 1940).

An unconfined aquifer is dewatered as the head is lowered by a discharging well. Because the dewatered portion of the aquifer cannot 

transmit water, the transmissibility decreases with the decline in head. 

Furthermore, the water is released from or taken into storage in 

response to head changes, largely by gravity drainage or refilling of the 

zone through which the water table moves. Although the coefficient
WATEE LEVELS NEAE A WELL IN UNCONFINED AQTJIFEE 29

of storage is usually assumed to be constant and to be virtually equal 

to the specific yield of an unconfined aquifer, the release of water from 

some rocks is not instantaneous, and some variation in this coefficient 

does occur in both time and space during a period of discharge.

Even if a constant value of S is assumed, the ratio T/S is variable 

because T is a function of the head in an unconfined aquifer. An exact 

solution of equation 1 with variable T/S has not been achieved in 

terms of elementary functions. It is the purpose of this paper to 

present a numerical solution for equation 1 applied to an unconfined 

aquifer for an assumed relation between T and head.

RESTATEMENT OF THE DIFFERENTIAL EQUATION

A saturated soil will drain very rapidly to a critical range of moisture 

content within which values of capillary conductivity become small; 

thereafter, the soil drains very slowly. At Seabrook, N. J., for example, 

sediments above the capillary fringe and below the root zone rapidly 

drain to, and then remain approximately at, this critical moisture 

content. Subsequent drainage is very slow (Remson, Randolph, and 

Barksdale, 1960, p. 150-151). Therefore, as a reasonable approximation, S can be taken as constant in equation 1 if applied over a time 

period as long as a few days, as is commonly the case in analysis of 

pumping-test data.

For convenience, let

D-l (2)

Because S is taken as constant, the hydraulic diffusivity, D, varies 

with T , which in turn varies with hydraulic head.

Let hn equal the constant hydraulic head throughout the aquifer 

before the well discharge begins, and let h 0 equal the head of water in 

a discharging well that commences to discharge at 2=0. It is assumed 

that the head is instantaneously lowered to h0 at t0 and is maintained 

constant throughout the discharging period; in other words, the well 

is "pumped at constant drawdown." It is assumed also that the head 

of water immediately outside the pumped well is hQ. Furthermore, let

/,_ h "o /o\ 0=7 T* vv nn n0

Then 6 is proportional to the difference between the hydraulic head at 

any point in the aquifer and the head at the discharging well, and is 

defined over the interval Q<,8^1. Thus, 6 is a function of head, 

which in turn is a function of radial distance from the well and time. 

At the well, 6=0} and at a remote point (r >»)J 0=1.
30 GROUND-WATER HYDRAULICS

Substituting equation 2 and appropriate differentials of equation 3 

in equation 1

r dr 

The boundary conditions for the case of a discharging well are

0=0« at t=0 and r>0 

and

0=00 at r=0 and 2^0 (5)

and 0«>00. It is assumed also that D is a single-valued function of 0. 

In applying equation 4 to unconfined flow, we are essentially adopting the assumptions of Dupuit (1863) as described by Jacob (1950, 

p. 378).

In discussing the flow of ground water in a vertical plane above a horizontal sole or 

impermeable bottom, he [Dupuit] assumed the velocity along the water table to be 

proportional to the slope or tangent of its angle of inclination, instead of to the sine 

as it actually is. Furthermore, he assumed the flow to be horizontal at the water 

table and everywhere below and the velocity to be uniform from top to bottom.

As shown by R. W. Stallman of the U.S. Geological Survey (written 

communication, 1960), equation 4 can be derived from Jacob's general 

equation for unsteady unconfined flow, which also is based upon the 

Dupuit assumptions (Jacob, 1950, p. 384).

Equation 4 has the same form as the radial-diffusion equation where 

diffusivity is a function of concentration. In the diffusion case, D 

would be the diffusivity and 6 would be the concentration. Philip 

(1955) has presented an iterative procedure for the numerical solution 

of the diffusion equation for one-dimensional flow, and his procedure 

may be adapted for solving equation 4. For clarity and convenience, 

much of Philip's development is reproduced below, most of his notation 

being retained.

DERIVATION OF THE ITERATIVE EQUATIONS

CONVERSION TO AN ORDINARY DIFFERENTIAL EQUATION

As a first step hi the derivation of the iterative equations, the partial 

differential equation 4 is converted into an ordinary differential equation by means of the Boltzman transformation

<£=rr 1/2. (6)

Q and D, both single-valued functions of r and t, are consequently 

single-valued functions of $. Therefore, equation 4 can be rewritten
WATER LEVELS NEAR A WELL IN UNCONFINED AQUIFER 31

When equation 6 and its derivatives are substituted into equation 7, 

the results can be simplified to the ordinary differential equation

_0 (de\_dL r dn D /de\ , .

2\dj>) <^L ^_r<#»w/ w

This equation is subject to the boundary conditions

0=00=0, at 0=0, 

and

0^»0n=l, as 0->oo, (9)

where 8o^d^6n for the case of a discharging well.

The boundary conditions state that, at any time, 6 equals 00 at the 

pumped well, where r=0 and 0=0. At any time, 6 remains at its 

original value, dn , at an infinite distance from the well, where r=«> 

and 0=oo. The second condition states also that, when 2=0, 

0=°= and heads are at their initial values throughout all parts of the 

system except at the well, where 0 is 00.

For a homogeneous aquifer, it may be assumed that 6 is a monotonically increasing function of 0, and equation 8 can be multiplied by

 TT and rearranged to give

2Z>_

 2

Because

as

the second boundary condition can be expressed in the form

Equation 10 can then be written in integral form incorporating the 

second boundary condition:

586293 61
32 GROUND-WATER HYDRAULICS 

FINITE-DIFFERENCE FORM OF THE DIFFERENTIAL EQUATION

A finite-difference form of equation 12 can be developed by dividing 

the interval from 00 to 6n into n equal steps of size 50. Let

ep=e0+p8e (13)

where p is the number of equal steps of size 50 between 00 and 6P .

The <£=/(0) curve can be replaced by a histogram with<£P and 6P as 

coordinates of the midpoint of step p. Then, by use of equation 13, 

the following finite-difference approximations are shown to hold:

f'»+* ,50 ,, 4 , 

##=0, T (14) Jov *

and

C6* Jn , 8d /1KN </)dO=(l>p -75- (15) 

 Vi 2

For 0=0p+i/2 , the finite-difference form for the first term in equation 

12 may now be written as

C gn C ev+* , C6n , 50 ,^N 0d0= (^^0 0^0= «^0 0P g- (16)

Let Dp+i/2 be the average value of Z? in the interval 6P to 0P+i, such 

that

J 9

"1 Die

de

Then, by analogy with equation 16, the finite-difference form for the 

second term in equation 12 is

The finite-difference form for the right side of equation 12 is

By substitution of equations 16, 18, and 19 into 12 and rearrangement,

, , _
WATER LEVELS NEAR A WELL IN TJNCONFINED AQUIFER 33

Let _ 

50 Jer L 9 J

i i i 

Then

In similar fashion, it can be shown that

 ip-l/2

where

50 J0,, L ° 

From equations 21 and 24,

p+i/2 _ -ip-i/2 

9P

When Dp+i/2+Dp-i/2 is replaced by 2DP,

Ip+i/2 Ip-i/2= 0p ~r^' (26)

Equations 22 and 26 form the basis of the iterative procedure for 

solving equation 12, the transformed version of the basic flow equation 4.

EVALUATION OF / CLOSE TO 0=0n

The iterative solution involves the calculation of 6 and 0 values 

over the interval 80 ^8^dn. At 0=00 , 9=0, and this value can be 

used directly. However, at 0=0re , $ ,= <*>, and this value can only 

produce other values of infinity when used to determine In-i/2- 

The following provides a suitable method for determining In-i/2-

A range of values, 0=0re_ e to 8=8n, can be selected such that the 

variation of D is small over this range. In such an interval, it can be 

assumed that D is constant. This is valid for small values of e because changes of water level are small in the vicinity of 6n . If e is 

taken as 1, 0w_ e is sufficiently remote from the well for D to be virtually constant. Therefore, Dn~i/2 will be taken as the constant value 

of D throughout the interval 8n-i to 8n, and it will be used to derive 

In-i/2, which is equal to /p+i/2 for p=n 1.

As indicated by equation 21, if -^"U
34 GROUND-WATER HYDRAULICS

can be evaluated, In-\i2 can be determined. Therefore, values of <j> 

between 0w_j and 6n will be determined to evaluate this integral. 

For constant diffusivity, equation 8 may be written

_ __ ^j? d?e Dn-i,2 de_m

In this case, the equation will be evaluated between the intervals

8=8n -it 4>=<l>n -i and e^en ,<f>=<i>rt . (28) 

Let

-=| <29>

(Frank, 1950). 

WT hen (29) is substituted in equation 27, one obtains

dw 0 7 , fa ,onv = = fa (oOJ 

w 2Dn -i/2 v

Upon integration, equation 30 becomes

(31)

where C\ is a constant of integration. When terms are rearranged, 

equations 29 and 31 become

A,-,/2 . (32) 

a<p <p

Upon integration, equation 32 becomes

Cp 4-Dn-l/Z

e=Ci\ - dfi+Cz (33)

J 9

Let

u=jf- (34)

4Z>ra -i/2 

When (34) is substituted in equation 33, one obtains

(35)

u 

where the boundary conditions are

e=Bn.ltu=un.l=-^- and e=en,u=un=-^-' (36)

4/>B _l/2 4^-1/2

If the second boundary condition is incorporated in equation 35, one 

obtains
WATER LEVELS NEAR A WELL IN TJNCONFINED AQUIFER 35

(37)

Because 0-»0n=l at a remote distance from the well, where

it may be assumed that tin=-=J^ = <» . Equation 37 can then be ex4ZV1/2 

pressed in terms of the exponential integral function

0-dn=C3Ei(u). (38)

When the first boundary condition is substituted for the beginning 

of the "n"th interval, one obtains the equation

4^ (39> 

or

sa

(40)

where 50=0« 0«_i. Thus, after both boundary conditions are satisfied, equation 38 becomes

. . 50

and, after terms are rearranged,

(42)

With known values of 0w-t> 6B, 0»_i, and Dn-i/2, equation 42 permits 

the computation of values of 0 corresponding to any values of 0 in 

the range #»_, to 0».

The integral in equation 21 cannot be evaluated because <£-»«> as 

0-»0w. However, assume that #»_3/4, 4>n-i/2, and 0«_i/4 have been 

evaluated by means of equation 42. The integral in equation 21 can 

then be replaced by the following open-quadrature formula:

J, (43)

'*-! «

(Kunz, 1957, p. 139),

where
36 GROUND-WATER HYDRAULICS

APPLICATION TO A HYPOTHETICAL PROBLEM

Assume that the pumping level in a discharging well is at the 

bottom of the aquifer. Although h0=0 is a poor choice so far as the 

Dupuit assumptions are concerned, it is convenient for illustrative 

purposes. If one lets n lQ and remembers that 60 Q and 0«=1, 

values of 6P can be tabulated as shown in table 1.

Assume that the aquifer has a transmissibility of 300 cubic feet 

per day per foot and its coefficient of storage is 0.30. From equation 

2, Dn= 1,000 cubic feet per day per foot. If D is assumed to be 

linearly related W 0, values of DP and DV-\K can be listed as shown 

in table 1.

The computation was begun by assuming a trial value of $9 . In

this case, 27.11 -j= was chosen arbitrarily. By use of equation 42,

tables of the exponential integral (Federal Works Agency, 1940), 

and the trial value of <£9 , it was found that 09M ==33.05, <£9M =41.36, 

and <£9M =55.02. These values were substituted into equation 43 and 

the result of that substitution into equation 21, giving /9M =42.36.

TABLE 1. First iteration

ep

Bo =0.0_ _ _ _ _ _ _

0t = .!.___._____. ___....._.

02 = .2-- - ___ ___ _ - ___ __

63 = .3 __ .__ ___ ._.__.__.

04 = .4_______. .___.__.____.

08 = .6__ ___ __ ___ _ ___

07 = .7 _ ___ ____ ___ _

08 = .8 ______ ___ ________

Bg= .9____- ._._._____._____

010=1.0 ________ __ _ __

Dp

0

100

200

300

400

500

600

700

800

900

1,000

VP-X

50

150

250

350

450

550

650

750

850

950

Ip-X

1,703

1,596

1,402

1, 155

890.3

640.5

427. 4

260. 1

135. 9

42. 36

4>>

1.84

1.90

2.09

2.45

3.06

4.07

5.79

8.83

14.60

27. 11

Substituting the values of /9^ and <£9 into equation 26, one obtains 

J8H=135.9, and substituting the values of <£9 and !$& into equation 

22, $8=14.60. By alternate use of these equations, the values of 

Ip-i/2 and <f>P in table 1 were determined.

Although $o in table 1 was computed as 1.84, it must equal 0 to 

satisfy the boundary requirements. Therefore, the original choice 

of fa was too large, and the computation was repeated with a smaller 

value of $9. When this procedure was repeated five times, $0 was obtained very close to 0 as shown in table 2. It is noteworthy that !<$& 

approaches 0 as <£0 approaches 0. A mathematical discontinuity 

makes it impossible to derive a value of 00 exactly equal to 0.
WATER LEVELS NEAR A WELL IN TJNCONFINED AQUIFER 37

The values of 0P and the corresponding values of <f> p lead to the 

solution. If time is fixed in <f>P, the variation of ground-water head 

with distance can be determined at a given time. If distance is fixed 

in <j>p, the variation of ground-water head with time can be determined 

at a given distance from the discharging well.

As discussed above, a linear relation is assumed between D and 6. 

Therefore, the solution of equation 4 forA0=0 and for any value of Dn

ft3 

= 1,000 X -I 57 can be determined by multiplying the respective

values of 4>p in table 2 by X1 '2 . The resulting curves relating 6 and 

<£ for A,0=0 and for different values of Dn are shown in figure 9. 

Values of <t> are plotted on a logarithmic scale to improve the usefulness of the curves. However, this method of plotting makes apparent the error in the approximation of the boundary conditions.

Assume that one or more observation wells are available in the 

vicinity of a well discharging from an unconfined aquifer. If the 

pumping level in the discharging well is at the bottom of the aquifer, 

the field curves relating 6 and $ can be compared with figure 9 and the 

value of Dn determined for the aquifer. Figure 10 is included to 

facilitate interpolation between the curves.

Da , IN FEET* PER DAY

0.8

0.7

A A

77

ID** I0~5 10'" 10"' 

<(>. IN FEET PER DAY*

FIGURE 9. Belation between 0 and <p for different values of D..
38 GROUND-WATER HYDRAULICS 

TABLE 2. Final iteration

«p

00 -0.0 ___ ----- __ -----

0i = .1

02= -2 _____ -- _______ --_ _

A ___ O 

ffn m\

04= .4___---_-__-___-_-_-__

05 = .5- ___ ____ - ____ __

0« = .6 _ -- _ _ __ -. _.

08= .8-_.-_ __ -__- __ --_-__

0» = .9. ___ _-_ _ _ _ ___

0io~1.0 -_ -

DP

0

100

200

300

400

500

600

700

800

900

1,000

D,-H

50

150

250

350

450

550

650

750

850

950

Ir-H

2, 634, 000

2, 132, 000

1, 390, 000

722, 300

294, 100

91, 180

20, 560

3,110

265. 1

.07112

*p

0.0004

.0004

.0005

.0009

.0019

.0049

.0170

.0802

.5625

6.9750

Dn , IN FEET 2 PER DAY

0.3

0.2

0.

0.0

IO'5 ~ 4

I0 , 

0, IN FEET PER DAY f

FIGURE 10. Enlarged portion of curves relating 0 and 0 for different values of D».
WATER LEVELS NEAR A WELL IN UNCONFINED AQUIFER 39

The assumptions made in applying equation 4 to unsteady unconfined flow limit its range of applicability. It must be emphasized that 

figures 9 and 10 apply only when the pumping level in the discharging 

well is at the bottom of the aquifer, or hQ =Q. This is a condition obviously difficult to attain and one that violates the Dupuit assumptions. However, similar curves can be developed for any values of 

AO. In fact, a single set of curves can be developed for use with any

value of At>-

REFERENCES

Dupuit, Jules, 1863, Etudes theoriques et pratiques sur le mouvement des eaux:

2d ed., Paris, Dunod, 304 p. 

Federal Works Agency, 1940, Tables of sine, cosine, and exponential integrals:

Work Projects Administration for the City of New York, 444 p. 

Frank, F. C., 1950, Radially symmetric phase growth controlled by diffusion:

Royal Soc. Proc., v. 201-A, p. 586-599. 

Jacob, C. E., 1940, On the flow of water in an elastic artesian aquifer: Am. Geophys. Union Trans., v. 21, p. 574-586. 

Jacob, C. E., 1950, Flow of ground water, ch. 5 in Rouse, Hunter, Engineering

hydraulics: New York, John Wiley and Sons, p. 321-386. 

Kunz, K. S., 1957, Numerical analysis: New York, McGraw-Hill Book Co.,

381 p. 

Philip, J. R., 1955, Numerical solution of equations of the diffusion type with

diffusivity concentration-dependent: Faraday Soc. Trans., v. 51, pt. 7, p.

885-892. 

Remson, Irwin, Randolph, J. R., and Barksdale, H. C., 1960, The zone of aeration

and ground-water recharge in sandy sediments at Seabrook, N.J.: Soil Sei.,

v. 89, no. 3, p. 145-156. 

Theis, C. V., 1935, The relation between the lowering of the piezometric surface

and the rate and duration of discharge of a well using ground-water storage:

Am. Geophys. Union Trans., v. 16, pt. 2, p. 519-524.

O