Pumping from elastic artesian aquifers

Inr. J. Engng Sci., 1972, Vol. 10, pp. 409-423.

PUMPING

Pergamon Press.

FROM

Printed in Great Britain

ELASTIC HALLDOR

Department of Mathematics,

Institute of Hydrodynamics

ARTESIAN

AQUIFERS

ELIASSONS

University of Bonn, West Germany

JONAS ELIASSON and Hydraulic Engineering, Technical University of Denmark, 2800 Lyngby, Denmark

(Communicated

by P. I. POLOUBARINOVA-KOCHINA)

Ah&act-The pumping from a well penetrating an artesian aquifer of constant thickness is studied. For the case when the pumping rate is a known function of time, the solution to the governing differential equation is found to be expressible as an integral transform of the pumping function. We compute the exact expression of this integral transform (formulas 25, 26, 27) and obtain, as a good approximation, the formula (38): nz f(l-r)exp(-($+r))$ 0

h(x,t)=i/

for the deviation h of the pressure level from the undisturbed value at time I and distance x from the well, wherefdenotes the pumping rate function. This approximation is good ifx is large compared with the radius of the well. We also show that this result contains the result of Theis and Jacob as a special case. THE

FLOW

EQUATION

of unsteady flow of groundwater is in most textbooks on groundwater flow. It is derived from Euler’s equation combined with Darcy’s equation and the equation of continuity. We end up with an equation of the form

THE EQUATION

(1) where h is the deviation of the pressure level from the undisturbed value, A denotes the Laplacian operator, t is time and E, (Y,j3 coefficients. In the derivation of equation (1) all quadratic terms have been neglected, and the aquifer is assumed to be elastic. The different terms in (1) have a special physical meaning each, viz. as follows: (a 2h/a 2’) represents the acceleration. This term is very small except at the very start of the motion. It is usually thrown away in further treatment: (a h/at) is the transient term, or the ‘storage’ term. It originates either from the lowering of a free groundwater surface or elastic consolidation of the aquifer. h is the inflow term, representing inflow from other aquifers. In order to incorporate the inflow term in (l), one has to assume parallel flow, otherwise the inflow is a boundary condition. The inflow is assumed to be zero for t < 0, but proportional to h for t > 0. This is meant to model an aquifer, originally in a state of equilibrium, but it starts drawing water from surface reservoirs (rivers, lakes) as the pressure within it drops. The coefficients of (1) vary from case to case. tCurrent Address:Mathematics Institute, University of Warwick, Coventry, England. 409

H. ELIASSON

410

and J. ELIASSON

In the following we shall confine ourselves to the study of a single well, (Fig. 1). The time is made dimensionless by dividing by a time scale t’, and length coordinates are made dimensionless by dividing by a length scale 1’. The coefficients become Sl’Z E=-

ngatr2

To 112

P=r’,

2’

The following symbols are used: S: storage coefficient (m3/m”/m) it: porosity a: thickness of aquifer (m) T: transmissibility (m3/s*m) T,: inflow coefficient (inflow at boundary: qO= he TJa2) (m3/s/m) g: acceleration of gravity. The relative magnitude of the dimensionless constants can be found by choosing the length and time scale so that CY and 0 become equal to 1, which gives E=

1 2+w+V’

The coefficient W is lp=-

nga3S TOT .

Fig. 1.

1 C 4.

(2)

Pumping from elastic artesian

aquifers

In practice, q is very large, since T and, especially, usually have l=9-‘4 NATURE

411

T,, are small quantities,

so we

1.

OF THE EQUATION

By denoting the dimensionless distance from the well by x = r/l’ and assuming axial symmetry, the flow equation becomes: (3) We seek a solution h satisfying the following conditions: (4.0) h is continuous in its interior.

in the domain - ~4< c < 03,b c x < CQand twice differentiable

ah(w) ax

-

+ -ii(t)

wherefis a given piecewise, continuous bounded for f + - m.

for

x+

function in - CQ< t < m,f(t) hs0

3

b

aya"h(t,x)

XatY

ax”

(4.1) 2 0 and remains (4.2)

c M(T)

(4.3)

b < x < 03,- ~4< r c 7 and V, p equal 0, 1 and 2, where A4 is some increasing function in--co
s

- xhdx

b

representing the mean deviation of the pressure level at time 1. The function v is well defined and twice differentiable in - ~4< f < m as a consequence of (4.0) and (4.3). Now, all terms in equation (3) are absolutely integrable over b s x c CQ;so, by performing the integration we obtain the following differential equation for v:

Eg+$+v=

j;;(x$)dx=f(r).

412

H. ELiASSON

and J. ELIASSON

Moreover, u(r)

c

m?dx=M(T)

M(r) b

--m
using (4.3). The roots it 2 m of the characteristic equation eXz + A+ 1 = 0 are both negative as a consequence of 0 6 E < 4, so that a solution of the homogeneous equation is not bounded for t + - ~0. Since v is to be bounded for t -+ - 00 we obtain a unique solution for 0: If0 < E < $: t

t

e+@“(r) dr - emt

--oc

f

_ m e-“‘f(r) dr) ,

IfI = 0:

u(t) = e-l

It,

eTf(7)dT.

Consequently, if h is a solution of (3), (4) withf(t) = 0, then v(t) = 0, and since h 3 0 this implies h = 0. Thus we have proved the uniqueness of the solution of (3), (4). As an interesting a priori estimate on the solution h, we obtain u(t) s C, supf(r) --m
tax t z&Ix

=

const.

t O,b

s= t-fl(x-

b) Fig. 2.

Pumping from elastic artesian aquifers

413

equations any more than uniqueness, because of the boundary conditions. Nevertheless, we shall show that the solution of (3), (4), which was found by H. Eliasson, has this property. EIGENFUNCTIONS

We seek the eigenfunctions variables, so we put

of equation (3) and use the method of separating the

h = exp (iot) . R(x)

(3

which gives $$+i.E+AR=O

(6)

h2=~~2-l-iiW,i=~

(7)

where

(6) is Bessel’s equation, and A is a complex number. It is, of course, possible to regard w as complex, but here it is chosen to assume real values, - ~0< w < m. The complete solution of (6) is:

R(x)= A . Hf’(hx) + B . Z-Zff’(hx).

(8)

A and B are complex constants, Hhl) and Hhz’ the Hankel functions of zero order, first and second kind, respectively. They have the following asymptotic expansions for large arguments: H$‘(z_) = J&exp

Hf’(z)

(iz-i$(l+O($)

= zexp(-iz+i~)(l+O(~)).

With z=hx,A=A,+iA2,Az we obtain Hh2’ (Ax) --, m for x + m. Therefore, function is

> 0

(9)

B must be zero, and thus the eigen-

h = A exp (iwt)Hg’(Ax)

(10)

and the relationship between A and o is A2= ew* - 1 First, we investigate

iw,

A = A1+ iA,, A2 > 0.

the case o = 0. This makes h independent

(11) of time, and thus

H. ELIASSON

414

and J. ELIASSON

f(t) must be a constant. This is the case of a constant pumping rate in time. For (11) gives A = i, and we have the stationary solution: h =A

6.1 =

* K”(X)

0

(12)

K,(x) being the modified Bessel function of zero order and 2nd kind. With proper dimensions, Q the pumping rate and r,, the radius of the well, this becomes:

(13)

The constant A is determined by (4. I), here KA denotes the derivative of K,,. K,(x) the following asymptotic expansions for small and large arguments:

has

K,,(x) = (ln2-C,,)-lnx+O(x)

zexp (-x)(1+0(-j). K,(x) = d

(14)

Next, we consider the case of o 4 0. Then we obtain for large x: h=A

This is a damped wave travelling with the wave celerity:

c,+>o

(15)

1

Wavelength. k

Frequency* M ‘27T’ The amplitude decreases exponentially with exp (-h2x). This is a good approximation far away from the well, but closer to it the wave characteristics are a slowly varying function of A, * x. - lhl’

*K,(x) Fig. 3.

415

Pumping from elastic artesian aquifers

At and AZare calculated from: g-g

= (5”2--

1,

h,A,=-;,

A*

>o

This gives (A, with opposite sign of 0): (V(Ew2-1)*+0~+~&-1)

A, =?

(16) A2=

(X+X+l)~++-Ew~+1)

A, and AZare the two length scales of the problem. A, is the wave number, and A, the exponential damping factor. We shall consider the two limiting cases: o + CQand o + 0, the case o + -m having the same character as w 4 m. First limiting case: o + m A, --+ -

AGO

X2-+ -

1

2VS

As E is a small parameter, we have oscillations with large wave numbers, great damping, and high wave celerity. This is the case where acceleration forces are strong, and is without interest in ordinary groundwater theory. Second limiting case: o + 0. AI+-?

2

AZ-+ 1 C, -+ 2.

(17)

We have oscillations with small wave number, moderate damping, and low wave celerity. The motion is independent of E, so this is the case where the storage effect is the predominant characteristic of the aquifer. When the inflow also is small, this is the case of ordinary groundwater theory. An interesting result is that C, = 2; this case will be discussed in greater detail in a later chapter. GENERAL

SOLUTION

We have seen that (10) is a solution to the flow equation (3), when - tQ< o c m, A related to w according to (11) and (16). Now we try to find a function A(w) so that the general solution may be written: h(r,x)

=

Jrn A(o) -uJ

exp (~#~)~(Ax) do

(18)

416

H. ELIASSON

and J. ELIASSON

where H is written instead of Hb”. (4.1) gives

-f(t)

= b

I_“,A (w) exp

(~~~)~’ (Ab) A dw.

By using the Fourier transform we get

By inserting this result in (18) we get h(t,x) =

j-w f(s)G(r-s,x)ds --m

(19)

where the kernel G is given by: G(T,x)

=-&

H(W

hbH’(hb)

exp (iW) do.

(20)

The computations have so far been quite formai, in the sense that no attempt has been made to prove the validity of (19) and (20). It remains to investigate the kernel G in order to justify the above computations and thus show that (19) is the solution to the boundary value problem (3), (4). The Hankel function H(z) is analytic in the z-plane with, e.g. the lower part of the imaginary axis z = iy, y 6 0 cut away. If we now define --1 Ao= J 4E

1

(21)

which is positive (>O) and real by (2) and then cut from the z-plane the interval of the real axis between -A, and A0

we have a conformal mapping of the cut A-plane on to the cut o-plane, where w is given in terms of A according to

Figure 4 shows how the square root should be taken. It also shows the path A given by the relationship (11) and (16). Using this transformation we obtain G(T.x)

=-&

exp (&)F(x,

with F (x. A) = H (Ax) /H’ (Ab) and o given by (22).

dA A) ~ EW- $2

(23)

417

Pumping from elastic artesian aquifers

1

=o =rX

=o

w

me

0

Re

e

Fig. 4.

For large values of IA1x we have the asymptotic F(x.h)

=

expression

db.1 x

i exp (ih(x--b))(l+O(&)).

(24)

It follows that the absolute value of the integrand in (23) is bounded by a constant multiplied by ~

1

1+ iA11

exp (A,Y) 3 with

y=$-(x-b).

First,weexaminethecaseofy t-&(x-_)in(19) (see Fig. 1). The integrand in (23) is an analytic function of A in the upper half plane, so the integral over the part A, of the curve A with IA 1 I s p is by Cauchy’s rule of integration equal to the integral over the curve T,,+ (Fig. 4), which consists of the line segments A1 = p, AZ= p, and A1 = p, starting at the point of intersection of A with Al = p and ending at the point of intersection of A with AI = -p. But the integral over I’,+obviously tends to zero for p --, CCin the above estimate of the integral, since y < 0. Thus we have shown: G(7.x)

=O.

for-7 < &(x-b).

This important result shows that the drawdown h (t. x) in an observation well is affected only by the pumping done before the time f - de (x - b). Now, let us examine the case of y > 0. Then we can replace the integral over AP with the integral over the curve - Tp, where Tp consists of the line segment A1 = - p starting at the point of intersection with A, the line segment AZ= -p and the line segment AI = p, ending at the point of intersection with A, where the negative, imaginary axis and the real interval -A, < A s A,,are avoided by circumferring it as described in

418

H. ELIASSON

and J. ELIASSON

Fig. 2 (at a distance tending to zero with p + m). Our estimate of the integrand in (23) shows that the integral over the line segments in IQ tends to zero for p -+ “Q;so in the limit we are left with the integral over the part of r; circumferring the negative imaginary axis and the real interval -A,, 6 A s ho. Denoting the first integral G,, and the second Gz, we get, after a straightforward computation, G = G, + GZ with

. K(xY)~‘(~Y)

-K’(~Y)~(xY)

.

K'(by)2+7?z'(by)2

.J(xY) Y’(~Y) - Yby)J’(by) .I'(by)Z+

Y'(by)Z

dy VgTy

du

.

vxgq.

(25)

(26)

Here we use the notations J and Y for ordinary Bessel functions of the first and second kind I and K for modified Bessel functions of the first and second kind. Before we proceed with our discussion of the kernel G, we shall investigate what properties a kernel function G must have in order that

h(t,x) =

j:l-‘)

- s,x)

f(s)G(r

ds

(27)

is the solution to our problem. Lemma. Letf(t) be a differentiable function in --co < t < CQ,which is bounded as t --, - w. Let G (t,x) be a twice differentiable function in the closed domain x 2 6, t 2 ~E(X - 6) and satisfy there the differential equation (3) and the boundary conditions aG (t, b) =O.

ax

for

r>O

G(r, x) s const e+, for some c > 0. Then the function h(r. x) defined by (27) satisfies the flow equation boundary condition (4.1) if - and only if G(&(x-b).x) Proof. Condition (29) guarantees sincefis differentiable. We have

cw (29) (3) and the

= (ebx)-:exp that h is well defined and twice differentiable,

419

Pumping from elastic artesian aquifers

ah

--&f(t-&(x-b)) ax-

+ Moreover,

t-V&'-b) f(s)

*G(&(x-b).x)

W-A

ds

ax

computing the other partial derivatives

*

of h, we obtain

whereu(x) = G(&(x-_),x). It follows that aGO, ax

b)

= 0, u(b)

= G(0,

b) = z

1

and

are necessary and sufficient conditions, in order that h may satisfy (4.1) and (3). But (30) is the unique solution for u Q.E.D. We shall now continue with our discussion of the kernel G = G1 + Gz. From the representation of G by (25) and (26), and using the asymptotic expansions for K and I, it is easily seen that G is infinitely differentiable in 7 > &(x - b), x 3 b, and satisfies there equation (3). Moreover, conditions (28) and (29) follow immediately. If we take forf(t) some infinitely differentiable function with compact support, then we can easily justify the computations leading to (19) and (20) and thus prove that h, given by (27), is the solution to (3) and (4.1). Then, by the necessity in the lemma above, it follows that G satisfies condition (30). Now, h is positive, since G is easily seen to be positive, and h satisfies (4.3), since h tends exponentially to zero for x + w as will be shown by the estimate of h in the next chapter. Thus we have proved: Theorem. The function h(t, x) given by (27), with the kernel G given by (25), (26), is the unique solution of the boundary value problem (3), (4), if we assume additionally thatf(t) is a differentiable function.

PRACTICAL

CALCULATIONS

The formulas (25), (26) for the kernel G are rather complicated and not well suited for numerical evaluation of the solution h(r, x), or even for investigations of its behaviour. Therefore, we shall compute several approximations for G and h under some simplifying assumptions in the domain of the coefficients of the differential equation. First, we assume the distance x from the well to be large compared to the radius of the well b. Then we obtain the corresponding approximation by letting b + 0 in (25)

420

H. ELIASSON

and J. ELIASSON

and (26). We get A”

10

G1(7.x)

I

=-+e-& E

G(T.x)

w)I(Xy) $(!Z_ = e-‘le

(

ew -+

0

V

E

= e-;‘lcosh-$v V

with v = d/(7” -

lx2)

”

VT

(3 1) E

.

Thus G becomes singular along the first characteristic, which for b = 0 is 7 = X&X, and does not satisfy the boundary condition along the characteristics required by our lemma. This was to be expected by physical reasons, since b + 0 means going to a point source. However, we still get a good approximation for h, since the singularity of G is integrable. From (27), with b = 0, we obtain, by substituting (3 1) for G: t-vz h(t.x)

=

=

=

f(s)e-&f

--m V(t-8)*-d

I

af(t-

fixcosh

cosh$ E

d(t-s)*-~2

ads

U)e-z5CoshU cash (h,&sinh u))du

0

I1i(f(t--1,)

+f(t-f2))e-rcoshYdU.

(32)

Here, t1 and f2 are t, = (icoshu+&*sinhu)x (33) tz= (icoshu-&*sinhu)x. From (32) we obtain the estimate h(t, x) s supf(~)&(x) which shows, since (32) is valid for large x, that the solution h tends exponentially to zero for x + ~0.In fact, (32) is the exact solution of (3), (4) with lim @(ah/ax) = -f(t), as x * 0. (32) is a time integral of the pumping rate function multiplied by a weight function. The weight function has maximum value for u = 0, the large values of u have little effect, as 99 per cent of the weight function area lies within the interval (34)

421

Pumping from elastic artesian aquifers

We have for u = 0 (35)

tl = t.2= Y. From (17) follows that this implies that

where A, and C, are the limiting values when 6.13 0, i.e. the stationary solution. For increasing U, rl increases and r2 decreases until u takes the approximate 1 u=-log,2

1 E

value

for 6 small.

This u-value will always be outside the interval (34), as x is supposed to be large. The convolution will then approximately operate as shown in Fig. 5.

Fig. 5.

We can now put E = 0 in (32), (33); then we have

(36)

The integral (32) now becomes h(r,x)

=~/mj(r-~eU)+f(r-~e--U)exp(-xcoshu)du.

(37)

0

This we transform into m

h(r.x)

=i

1 f(r-7)

exp (-(g+T))$.

0

This integral is easily evaluated on a computer.

Forf=

(38) constant beginning at r = 0,

H. ELIASSON

422

and J. ELIASSON

this becomes

=$[ exp(-($+7))$.

h(f,x)

(39)

The case of zero inflow can be evaluated from this formula by letting 1’ and t’, the length scale, and the time scale tend to infinity in such a way that T .l”zzT.

a.

0

l’*

T

P

s

-=-

2

.

Then Y/T is fixed, but 7 --, 0; the integral is then transformed h(t.x)

into

;ydu

=f

(40)

41

This is the classical Theis formula, the integral oc

I e-U&, u

is called W(X)

s

but - Ei(- x) in most mathematical tables. It tends to m as x ---*0 and goes rapidly to 0 for increasing x. With proper dimensions we have from (4) 2nb$

T

=2rfT=Q f=TO

and furthermore: 2 _-_.......-=-. ? t’

t-

4t

Sr2

II2

4tT

So, finally, we get m

h(t.x)

=

&

I

Cdu. ST2 u ZF

(41)

It has thus been demonstrated how practical formulas may be evaluated from the complete solution of the flow equation. Acknowledgmenrs-This work was supported by the University of Iceland Research Institute and the National Energy Authority of Iceland. The authors also wish to thank professor F. A. Engelund of the Hydraulic Laboratory, Technical university of Denmark for his contributions during preparation of the manuscript. REFERENCE

[I] P. I. POLOUBARINOVA-KOCHINA,

Theory of Ground Water Movement, Princeton (1962). (Received7

October

1971)

423

Pumping from elastic artesian aquifers

R&sum& Cet article examine les probltmes poses par la pompage d’un puits pCnCtrant dans une nappe aquif&e artesienne. Dans le cas oh le taux de pompage est une fonction connue du temps, on peut exprimer la solution de l’tquation differentielle directrice sous la forme dune transformation integrale de la fonction du pompage. Nous calculons l’expression exacte de cette transformation integrale (formules 25,26 et 27) et nous obtenons une bonne approximation de la formule (38): h(x,t)

=iSf(t-r)exp(-($+r)):

pour l’ecart h du niveau de pression de la valeur non perturbee au moment t et a la distance x du puits; dans cette formule, f indique la fonction du taux de pompage. Cette approximation est suffisante si x a une valeur suffisamment importante par comparaison avec le rayon du puits. On remarque Cgalement que ce resultat contient le resultat de Theis et Jacob dans le cadre d’un cas special. ZusammenfassungEs wird das Pumpen von einem Brunnen untersucht, der eine artesische wasserfiihrende Schicht konstanter Dicke durchdringt. Fur den Fall dass die Pumprate eine bekannte Funktion der Zeit ist, wird gefungen, dass die Losung zu der regierenden Differentialgleichung als Integraltransform der Pumpfunktion ausgedriickt werden kann. Wir errechnen den genauen Ausdruck dieses Integraltransforms (Formeln 25. 26.27) und erhalten, als gute Annaherung. die Formel (38): h(x.r)

=iJmf(t-~) exp(-(g+T))+ 0

filr die Abweichung h der Druckhohe vom ungestorten Wert zur Zeit t und beim Abstand x vom Brunnen, wobeifdie Pumpratenfunktion bezeichnet. Diese Annaherung ist gut wennx im Vergleich mit dem Brunnenradius gross ist. Wir zeigen such, dass dieses Resultat die Resultate von Theis und Jacob als Spezialfall einschliesst. Sommario- Si studia il pompaggio di un pozzo the penetra una falda acquifera artesiana di spessore costante. Quando si ha il ritmo di pompaggio the corrisponde a una data e nota funzione di tempo, si trova the la soluzione dell’equazione differenziale regolatrice e esprimibile come una trasformazione integrale della funzione pompante. Computiamo l’espressione esatta di questa trasformazione integrale (formule 25, 26, 27) ed otteniamo con buona approssimazione la formula (3 8): h(x.r)

/:f(t-~)

=;

exp(-($+r))$

per la deviazione h de1 live110 di pressione dal valore indisturbato al tempo t e distanza x da1 pozzo, dove .f denota la funztione de1 ritmo pompante. Tale approssimazione e buona se x e grande in co&onto al raggio del pozzo. Dimostriamo pure the quest0 risultato contiene il risultato di Theis e Jacob come case speciale. A~CT~IICT-

nnacr BaeTCR

ki3y’ieHO

IIOCTO~HHOL B03MOxHO

IIbHOrO

OTKa’iWBaHl(e

BblpaXCKaTb

ll&PZO6pa30BaHHn

A3

~,“eHHe

I$)‘HKUHW

OCHOBHOrO

,LWI

OTKJIOHeHWR

h

YPOBHX

X OT CKBaXGiHbI,

X BeJIUKOM

B CpaBHeHliH

TaT

Ii IiKo6a

Tswca

AaBJIeHUR rAe

f

B

OT HeHapylUeHHOrO PaCXOna.

CKBaXGiHbI.

YaCTHOrO

+YHKL~~SI ,‘paBHeHI

TOYHOl-0

BOnOHOCHbIfi BpeMeaH,

B BAne

BbIpameHHR

oKa3btHHTerpa-

3TOI’O

npeo6-

tca’lecrae xopomeg annpokcnMauwu ($opMyna 38):

I-&JHKUHff

C PanHyCOM

B KaWCTBe

apTe3HaHCKHk

J,H&$e,,eHUUaJIbHOrO KOMIIfOTaLWR

OTKaYUBaHUR.

pa3oaaaus (+opMynbr 25,26,27) ,naer

CTORHWH

IIpOHHKalOI&i

CKBaXWHbI,

B cnyqae xoraa pacxoa ecrb n3aecrnas

M~LIJH~CTI~.

3Ta

nOKaxeM

CJIyqaX.

3Ha’ieHWII

B MOMeHTe

aIIIlpOKCAMaLOffl TaKxe,

‘IT0

BpMeHIi

IlBJllleTCR

pe3j’flbTaT

t ti Ha PaC-

XOpOIUOfi,

COnepEHT

EJIA

A pe3)‘Jlb-