THE MATHEMATICS OF
BASIC PUMPING WELL HYDRAULICS
John R. Hoaglund, III Ph.D.
Table of Contents
1.0 DERIVATION OF THE GOVERNING EQUATIONS FOR GROUNDWATER FLOW 1
1.1 Darcy's law and Related Hydrology 1
1.2 Continuity and the Laplace equation for Steady State flow 1
1.3 The Diffusion Equation for Transient Groundwater Flow 3
1.4 Reduction of the Laplace and diffusion equations to two dimensions: the governing equations for confined and unconfined groundwater flow 3
2.0 DERIVATION OF THE GOVERNING EQUATIONS FOR A PUMPING WELL: TRANSFORMATION OF THE GROUNDWATER FLOW EQUATIONS TO RADIAL COORDINATES 4
2.1 The radial form of the two-dimensional, confined diffusion equation 5
2.2 The radial form of the two-dimensional, unconfined Boussinesq equation 8
3.0 BOUNDARY CONDITIONS FOR THE PUMPING WELL MODELS 9
4.0 SOLUTIONS TO THE STEADY STATE PUMPING WELL MODELS: THE RADIAL LAPLACE EQUATION 9
4.1 Confined case 11
4.2 Unconfined case 13
4.3 Remarks 14
5.0 SOLUTIONS TO THE TRANSIENT PUMPING WELL MODELS 14
5.1 Confined Case: The Theis Solution 14
5.1.1 The Product Solution Strategy 16
5.1.2 The Similarity Solution Strategy 16
5.1.3 The Similarity Solution 18
6.0 TAYLOR SERIES EXPANSION OF THE EXPONENTIAL INTEGRAL AND THE COOPER JACOB APPROXIMATION 25
APPENDIX I THE NEED FOR A FULLY THREE-DIMENSIONAL UNCONFINED STORAGE TERM 27
REFERENCES CITED 29
THE MATHEMATICS OF BASIC PUMPING WELL HYDRAULICS
In hydrology, solutions for a pumping well in confined and unconfined porous media are derived from the diffusion and Boussinesq equations for groundwater flow respectively. These governing equations are second order partial differential equations derived from Darcy's law and continuity. Derivations of the solutions to these equations for steady state and transient flow conditions highlight the underlying simplifying assumptions: pumping a fully penetrating well in either an unconfined or confined aquifer that is isotropic and homogeneous.
Darcy's law relates the specific discharge (groundwater flux) linearly to the gradient of the hydraulic potential. The elevation of the water at any point in the spatial domain is known as hydraulic head. Hydraulic head has been related to fluid potential in a classic paper by Hubbert (Hubbert, 1940). As a spatially independent potential variable, hydraulic head is analogous to temperature in thermal diffusion and voltage in electrical conductance. Darcy's law in the x direction is written:
The proportionality constant, K, is known as the conductivity or permeability (analogous to thermal diffusivity in thermal diffusion and conductance as the reciprocal of resistance in electrical conductivity). It is related to the intrinsic permeability, k, by the following:
Since the groundwater flux can be thought of as a line length of water delivered in a unit time, and since the gradient is dimensionless as a change in head over a change in spatial length, the units of permeability, K, are length per time. Also, if groundwater flux, q, is multiplied by an equipotential surface area (at right angles to the flux), a volume per unit time (discharge Q) results:
If a unit cube is drawn with sides parallel to the x, y, and z axes, conservation of mass requires that any flux into the cube must equal the flux out of the cube. This geometrical consideration is known mathematically as divergence (net outflux) of a vector field and results in the following:
Note that the divergence operator (the second partial derivatives with respects to the independent variables x, y and z) operates on the directional permeability terms as well as the gradient, both of which are functions of the independent variables. This mathematical difficulty can be relaxed with the following condition:
Condition 1. The medium is homogeneous
If one assumes that the medium is homogeneous (note that this refers to the physical definition of homogeneity as opposed to the mathematical definition in differential equations in which terms involving the dependent variable and its derivatives are collectively set equal to zero), the directional permeabilities are not functions of the independent (spatial) variables and can be assumed as constant coefficients:
A second condition simplifies the algebraic manipulation of the directional permeability coefficients:
Condition 2. The medium is isotropic
If one assumes that the medium is isotropic, the permeabilities are not directionally dependent. This implies that Kxx = Kyy = Kzz = K.
The two simplifying assumptions applied together result in the Laplace equation for steady state flow.
Theis (Theis, 1935) defined groundwater storage analogously to capacitance in electrical conductance and heat capacity in thermal diffusion. In three dimensions, the specific storage, Ss, is defined as the volume of water released from a volume of aquifer per unit decline in hydraulic head. The concept of storage enables time dependent changes in head as the unit volume of aquifer obtains, retains, or releases water in storage. Continuity requires that a flux into the volume be balanced by either a flux out of the volume or a time dependent change in hydraulic head dependent upon the storage coefficient. This results in the diffusion equation:
A common simplification of hydrogeological problems is to assume that the vertical gradient and corresponding flow is zero resulting in horizontal, two-dimensional flow. This is usually justified based on the layered nature of aquifers as well as very high horizontal dimension to vertical dimension ratios. Two geometric cases, confined and unconfined, must be considered.
In the confined case, the aquifer with thickness b is confined between impermeable layers above and below with hydraulic head above the top of the aquifer. This is the case for the Theis (Theis, 1935) solution. The coefficients, K and Ss, are multiplied by this thickness due to geometric considerations resulting in T (Transmissivity) and S (Storativity).
The general form of the diffusion equation in two dimensions is
The isotropic homogeneous form in two dimensions is
In the unconfined case, the top of the aquifer is bound by the water table. With horizontal flow assumed, this water level corresponds to the hydraulic head at any point in the saturated aquifer. These conditions are called the Dupuit assumptions. The coefficient K is multiplied by the water level h resulting in an h-dependent, two-dimensional transmissivity. The two dimensional, unconfined storage coefficient is not distinguished from the three dimensional, unconfined storage coefficient.
Replacing the transmissivity coefficient in the two-dimensional diffusion equation results in an equation that describes the two-dimensional flow in an unconfined aquifer. This equation is known as the Boussinesq equation (Bear, 1979). The general Boussinesq equation is
The isotropic homogeneous form is
With the isotropic and homogeneous conditions stated above, a transformation from Cartesian to radial coordinates can be accomplished assuming no angular () dependence. This is done as follows:
Assume radial flow to a point (the origin) where
and h(x,y) = h(r). Then
and
The second partial of h with respect to x is
and similarly the second partial of h with respect to y is
Noting the partials r with respects to x are
and
and similarly the partials r with respects to y are
and
we can substitute into the expressions for the second partial of h with respects to x and y to obtain
Therefore, the diffusion equation becomes:
An alternate derivation of this radial flow equation starting from
Darcy's law is provided by Lohman (Lohman, 1972).
The homogeneous isotropic Boussinesq equation
is expanded using the product rule
Substituting the relations solved above between the first partial derivatives of h with respect to x and y to the first partial derivatives of h with respect to r, and substituting the relations solved above between the second partial derivatives of h with respect to x and y to the second partial derivatives of h with respect to r results in the following
The relation between r and x and y simplifies the equation to
Finally, if the non-linear term, (h/r)(h/r), is assumed to be negligible, as squares of infinitesimals often are, we obtain
1.) An initially flat hydraulic head, h0, is assumed for the entire domain at time t=0. This is expressed mathematically as follows:
h(r,0) = h0
2.) The hydraulic head at infinity is h0 for all time t. This is expressed mathematically as follows:
h(,t) = h0
3.) A constant pumping rate at r=0 for all time t implies that the flux through a cylinder of any radius r is equal to a constant.
For confined conditions, this is expressed mathematically as follows:
For unconfined conditions, this is expressed mathematically as
If no time dependent changes in the hydraulic head (steady state) is assumed or known, both the radial diffusion and radial Boussinesq equations reduce to the radial Laplace equation. The two solutions then differ only by imposition of different statements of the constant pumping rate boundary condition. The radial Laplace equation is
We begin the solution of this ordinary differential equation using reduction of order (Boyce and DiPrima, 1986) by letting
and
Substituting into the radial Laplace equation, we obtain
Separating variables yields a form we can integrate
Completing the integration and exponentiating both sides yields a solution of "v" with an integration constant written in terms of a natural log of a constant.
The integration constant becomes a coefficient and the definition of v is restored
The steps to solving both the confined and unconfined case are similar and include:
1. Applying the pumping boundary condition to evaluate C1,
2. Substituting for C1 in the expression for v which results in a restatement of the pumping boundary condition,
3. Separating variables and integrating to obtain a second coefficient, and
4. Evaluating the second coefficient with the other boundary condition.
Alternatively, one can integrate the pumping boundary condition between limits that include the other boundary condition and the solution corresponding to the radius of interest. This is done in the derivation of the solution for the unconfined case below. The derivation of the confined case uses the steps outlined above to demonstrate how the boundary conditions constrain the integration constants.
# 1. Applying the pumping boundary condition
for a confined aquifer yields
for all r > 0. Substituting the boundary into the expression for v gives
#2 Substituting C1 results in a restatement of the pumping boundary
#3 Separating variables gives a first order differential equation for h
which is solvable for any r > 0. Integrating gives
or
#4 Since h = h0 at r = ,
Therefore the difference:
Since C2 - C2 is 0, the solution is
The pumping boundary for unconfined conditions is
Separation of variables yields the integrable form
Integrating gives
Evaluating the limits yields the solution
An interesting result of the steady state solutions of a pumping well is that they involve logarithmic functions of the independent variable r. It may therefore be more appropriate to refer to the shape of the potentiometric surface around a pumping well as a "logarithm (log) of depression" rather than a "cone of depression." For these solutions to be useful, a practical "finite" radius of zero drawdown must be substituted for the infinite radius of zero drawdown and steady state must be assumed. Since this radius is really an unknown, it is best to determine aquifer parameters using the Theim equation for steady state pump tests involving two observation wells, or from a transient pump test (see transient solution below). Alternatively, if the aquifer parameters have already been determined from either the Theim equation or a transient solution, observations of steady state drawdown at a known radius can be used to determine a practical, steady state radius of influence (zero drawdown) for the pumping well. This radius would be useful for conditions of an initially flat regional potentiometric surface. For conditions involving a regional gradient of the potentiometric surface, it is best to solve for capture zones. For an example using superposition of radial and one-dimensional flow fields, see page 121-123 of Todd (Todd, 1976).
C. V. Theis, with the help of the mathematician C.I. Lubin, introduced the concept of storage in hydrogeology and produced a solution for the confined transient pumping well model which he published in Transactions of the American Geophysical Union in 1935. The full derivation of the solution was not published since Lubin knew of similar solutions in the mathematics and heat transfer literature. The derivation requires knowledge of boundary and initial value problems applied to ordinary and partial differential equations. A concise derivation is given by Li (1972). An expanded derivation assuming knowledge of calculus and following that of Li (1972) is presented here starting with the radial form of the diffusion equation.
The first step is to convert the radial diffusion equation to a similar drawdown form following the chain rule, where
s = h0 - h
noting
Since each term is negative, the radial equation in drawdown form is similar to the radial equation in head form.
The initial and boundary conditions must also be converted:
Head initial condition Equivalent drawdown initial condition
1.) h(r,0) = h0 1.) s(r,0) = 0
Head Boundary conditions Equivalent drawdown boundary conditions
1.) h(,t) = h0 1.) s(,t) = 0
2.)
2.)
The left side of the radial equation in drawdown form can be expressed as a product by noting
The product drawdown radial form is then
This partial differential equation is solved using the drawdown initial and boundary conditions. The first goal is to reduce the partial differential equation to a set of ordinary differential equation(s). This is often accomplished with a technique known as separation of variables. The technique assumes that the final solution can be expressed as a product of two functions, one dependent only on the first independent variable, r, and the other dependent only on the second independent variable, t. This assumed product solution is then differentiated according to the needs of the partial differential equation. The differentials are substituted into the partial differential equation, and two ordinary differential equations, corresponding to each of the two independent-variable functions, are separated and solved independently. This technique results in solutions involving Bessel functions and was not used by Theis.
Alternatively, the partial differential equation can be solved by assuming that the solution can be expressed in terms of a function of a new independent variable that combines the two existing independent variables, r and t. Such a function is called a similarity function and the resulting solution is called a similarity solution. By changing variables, substituting respective partial derivatives of the new similarity function into the partial differential equation, the partial differential equation is reduced to an ordinary one.
Choosing the new function and independent variable requires some mathematical ingenuity. Since the second boundary condition, r (ds/dr), equals Q/2T, we would like this "constant coefficient" in front of the new function so that it cancels conveniently during the change of variables of the second boundary condition. The transformation equation is thereby assumed to be:
The independent variable, "u", is chosen by inspecting the partial differential equation. Changing variables in the partial differential equation requires partial derivative of "s" with respects to the independent variables "r" and "t." Equivalently, these partials involve taking the partials of the similarity function, "f" with respects to its independent variable, "u", which, by the chain rule, requires partials of "u" with respects to the independent variables "r" and "t."
Since "r" and the chain rule partial of "u" with respect to the independent variables "r" is not a function of u, they become "constants" when the second partial derivative of "f" is taken with respects to "u," as required by the chain rule for the second derivative of r ds/dr. The r's and Q/2T's cancel and the partial differential equation becomes:
If we can find the independent variable u as a function of r and t where r, t, S, and T will cancel in the relation
the problem will reduce to ordinary differential equation. One such function is:
Other functions may work also, which underscores a problem with similarity solutions: they are not unique.
We let
and
and change independent variables in the partial differential equation to "u" forming an ordinary differential equation. To do this, we differentiate s in terms of the similarity function taking the necessary chain rule derivatives, and substitute into the right and left sides of the drawdown-radial form of the Lohman partial differential equation.
Right Side
Left Side
Note the first partial of "f" with respects to "u" requires
the partial of "u" with respects to "r:"
Note that r multiplied by the partial of "s" with respects to "r" results in "u" multiplied by partial of "f" with respects to "u:"
Note the second partial of "f" with respects to "u" requires
the partial of "u" with respects to "r:"
Finally for the left side we have
Canceling terms on the left and right side of the equation results
in:
which is equivalent to
or
which is an homogeneous, second order, linear, ordinary differential equation.
The initial and boundary conditions must also undergo the change of variable. Both the initial condition and the first boundary condition,
s(r,0) = 0
and
s(,t) = 0
state that as "u" goes to infinity, f() goes to zero; "u" goes to infinity as "t" goes to zero and/or "r" goes to infinity, and "s" is zero for this "r", "t", and "u," so that the function f(u) is zero.
f() = 0
The second boundary condition
is modified from s = f(u) Q/2T as follows
or
The change of variable is complete and we now have a boundary value problem involving an ordinary differential equation, restated as follows:
with boundaries
f() = 0
and
To solve the ordinary differential equation, we use an integrating factor to reduce the order of the differential equation by turning the left hand side of the equation into a derivative of a product. Integrating factors are discussed in Boyce and Di Prima (1986) of the form
where p(t) is the coefficient function of the independent variable u:
Therefore,
When the ordinary differential equation is multiplied by the integrating factor, we obtain
or
which can be written as the derivative of a product:
Integrating, we obtain
Applying boundary 2 for u = 0,
we evaluate the constant,
We are left with a first order ordinary differential equation
or
To solve, we integrate over the domain from u = 0 to u = which is equivalent to the domain r = 0 to r =
The integral over the domain is known as the zeroth Gamma function (Erd-elyi et. al., 1953-1955). Since the Gamma Function is singular at u = 0, we integrate over the sub-domain of the Incomplete Gamma function (Erd-elyi et. al., 1953-1955) defining the function f(u) at the lower limit.
The integral definition of the Incomplete Gamma function is a form of the exponential integral (Abramowitz and Stegun, 1972). Note that we applied the first boundary condition, f() = 0.
We are left with
Finally we change variables back using f(u) = s 2T / Q
and obtain the Theis solution
Using the Theis solution requires an approximation for the exponential integral. First, the exponential function is approximated by a Taylor series expanded about x = 0 (i.e. the Maclaurin's series).
Each term in the series is divided by the denominator "x." The series is then mathematically assumed to be integrable term by term.
Cooper and Jacob (1946) derived an equation for semilog linear plots of the Theis solution by evaluating limits of the well function integral and truncating the remaining series. According to Todd (1959), the Cooper Jacob well function:
is valid for 0<u 0.01 (Freeze and Cherry, 1979). Delillo (1993) reproduced the approximation by evaluating the well function integral as the sum of two integrals:
The final derivation of the Cooper Jacob form of the Theis equation requires some logarithmic manipulation of the -.5772 term left from evaluating and integrating the Taylor series (Freeze and Cherry, 1979). Noting that -.5772 = -ln(1.78):
Noting that ln(x) = 2.3log(x) (see figure 1) and substituting for the exponential integral in the Theis solution, we obtain the Cooper Jacob equation for semilog (base 10) plots:
The Need for a Fully Three-Dimensional Unconfined Storage Term
Three-dimensional modelling, either mathematical or numerical, of an unconfined aquifer is hindered by the lack of an analytically correct storage term. The field definition and laboratory definition of the unconfined storage term, the specific yield, do not exactly concur. Furthermore, specific yield should be restricted to two-dimensional modelling.
The three-dimensional storage term in use today, the specific storage, was derived from the analysis of aquifer compaction of, and water expansion upon removal from, a confined aquifer which remains fully saturated. It was very easy to derive a two-dimensional, areal, confined aquifer storage term, the storativity (or storage coefficient), by multiplying the specific storage times the saturated thickness, in an analogous fashion with the conversion from hydraulic conductivity to transmissivity. The resulting storativity term was quite useful in problems involving radial flow (or horizontal flow) to a well in a confined aquifer, as in the classic Theis-Lubin solution. This solution, unfortunately, is often misused in problems of three-dimensional flow in confined aquifers.
In unconfined aquifers, the most commonly used term for modelling is the specific yield, also known as the unconfined storativity (Freeze and Cherry, 1979, p. 61). This second name is fitting since it reminds us of the two-dimensional modelling term, storativity, in confined aquifers. In fact, they are defined quite similarly in the field hydrogeologic context. In this context, the storativity is defined as "the volume of water that an aquifer releases from storage per unit surface area per unit decline in the component of hydraulic head normal to that surface" (Freeze and Cherry, 1979, p. 60). Freeze and Cherry also illustrate it as the volume of water released per unit surface area of aquifer per unit decline of the potentiometric surface (Freeze and Cherry, 1979, p. 60). In the same field context, the specific yield term is defined as "the volume of water that an unconfined aquifer releases from storage per unit surface area of aquifer per unit decline in the water table" (Freeze and Cherry, 1979, p. 61). This is also illustrated in Freeze and Cherry (1979, p. 60). In the special case of Dupuit assumptions with strictly horizontal flow (with no change in head with depth) when the water table itself becomes the potentiometric surface, the storativity and specific yield (unconfined storativity) is a two-dimensional modelling term. This leaves us an ambiguity in picking a term for the three-dimensional unconfined case.
The specific storage term can be defined as the volume of water that a unit volume of aquifer releases from storage under a unit decline of head. Since it is volume per volume per unit decline of head, it has the dimensions of L-1. The two-dimensional terms, storativity and unconfined storativity (specific yield), are dimensionless since they denote a volume per area per unit decline in head, or the water table respectively. The laboratory definition of specific yield can be misleading in the field or regional hydrogeologic context. It states that "specific yield of a rock or soil is the ratio of (numerator): the volume of water which the rock or soil, after being saturated, will yield by gravity, to (denominator): the volume of the rock or soil. The definition implies that gravity drainage is complete" (Lohman, et. al., 1972 p. 12; from Meinzer, 1923, p. 28). This term is equal to the porosity minus the specific retention. The volume per volume definition here, however, does not include a reference to a decline in head, water table, or even water level in the laboratory sample. In the field, the volume of water yielded by gravity is a volume generated by taking a unit area of aquifer and lowering the water table, and the volume of rock or soil from which the water was derived may be implied to be the same "aquifer area--drawdown" cube constructed for the water volume assuming an even drop in the water level. We could take the volume of water obtained and divide it by the volume of rock we suspect it came from (this would just be the porosity less specific retention), but we would need to divide this term by the water decline to obtain a fully three-dimensional storage term. Therefore, the laboratory definition is not directly concurrent with the field definition of specific yield which relates volume per unit surface area of aquifer per unit decline in the water table. Furthermore, both the lab and field definitions of specific yield are dimensionless and are thus not fully three-dimensional storage terms, which should have L-1 modelling units.
To solve the units problem, one might be tempted to use the specific storage term. However, the aquifer compaction and water volume expansion phenomena are different from the dewatering phenomenon. Lohman et. al., in their definition of specific storage, reserve the term, "[for use] in problems of three-dimensional, transient flow in a compressible groundwater body" (1972, p. 13). The compressibility term in unconfined aquifers is considered negligible compared to the dewatering term. In two dimensions, the unconfined storage term can be expressed as follows:
be defined.
REFERENCES CITED
Abramowitz, M. and I. Stegun, 1972. Handbook of mathematical functions with formulas, graphs, and mathematical tables. U.S. Government Printing Office, Washington D.C.
Bear, J., 1979. Hydraulics of groundwater. McGraw-Hill Inc. New York. 569 p.
Boyce and DiPrima, 1986. Elementary differential equations and boundary value problems. Wiley and Sons. New York. 654 p.
Cooper, H.H.,Jr., and C.E. Jacob, 1946. "A generalized graphical method for evaluating formation constants and summarizing well field history." Transactions of the American Geophysical Union, Vol. 27, p. 526-534.
Erd-elyi, A., W. Magnus, F. Oberhettinger, and F. Tricomi, 1953-1955. Higher transcendental functions. Based, in part, on notes left by Harry Bateman, and compiled by the staff of the Bateman Manuscript Project, Vol. 2. McGraw-Hill Inc. New York.
Freeze, R.A. and J.A. Cherry, 1979. Groundwater. Prentice-Hall, Inc. Englewood Cliffs, New Jersey. 604 p.
Hubbert, M.K., 1940. "The theory of groundwater motion." Journal of Geology, Vol. 48. p. 785-944.
Li, Wen-Hsiung, 1972. Differential equations of hydraulic transients, dispersion, and groundwater flow; mathematical methods in water resources. Prentice-Hall, Inc. Englewood Cliffs, New Jersey. 316 p.
Lohman, S.W., 1972. Ground-water hydraulics. U.S.G.S. Professional Paper 708. 70 p., 9 plates.
Meinzer. O.E., 1923. "The occurrence of groundwater in the United States, with a discussion of principles." U.S.G.S. Water-Supply Paper 489.
Theis, C.V., 1935. "The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using groundwater storage." Transactions of the American Geophysical Union, Vol. 2. p. 519-524.
Todd, D.K., 1959. Groundwater hydrology. John Wiley and Sons, New York.
Todd, D.K., 1980. Groundwater hydrology, second edition John Wiley and Sons, New York. 535 p.