Chapter 8: Groundwater Resource Evaluation

Groundwater Resource Evaluation

In the first seven chapters of this book we have examined the physical and chemical principles that govern groundwater flow and we have investigated the interrelationships that exist between the geological environment, the hydrologic cycle, and natural groundwater flow. In this chapter and the two that follow, we will turn to the interactions between groundwater and man. We will look at the utilization of groundwater as a resource, we will examine its role as an agent for subsurface contamination, and we will assess the part it plays in a variety of geotechnical problems.

8.1 Development of Groundwater Resources

Exploration, Evaluation, and Exploitation

The development of groundwater resources can be viewed as a sequential process with three major phases. First, there is an exploration stage, in which surface and subsurface geological and geophysical techniques are brought to bear on the search for suitable aquifers. Second, there is an evaluation stage that encompasses the measurement of hydrogeologic parameters, the design and analysis of wells, and the calculation of aquifer yields. Third, there is an exploitation, or management phase, which must include consideration of optimal development strategies and an assessment of the interactions between groundwater exploitation and the regional hydrologic system.

It is worth placing these three phases in a historical perspective in North America and Europe, nearly all major aquifers have already been located and are being used to some extent. The era of true exploration for regional aquifers is over. We are now in a period in which detailed evaluation of known aquifers and careful management of known resources will take on greater importance. The layout of this chapter reflects this interpretation of current needs. We will treat aquifer exploration in a single section, and place heavier emphasis on the evaluation and management stages.

Let us assume that we have located an aquifer that has some apparent potential. The scope of groundwater resource evaluation and management studies might best be indicated by the following series of questions:

  1. Where should the wells be located? How many wells are needed? What pumping rates can they sustain?
  2. What will be the effect of the proposed pumping scheme on regional water levels?
  3. What are the long-term yield capabilities of the aquifer?
  4. Will the proposed development have any detrimental influence on other components of the hydrologic cycle?
  5. Are there likely to be any undesirable side effects of development, such as land subsidence or seawater intrusion, that could serve to limit yields?

This chapter is designed to provide the methodology needed to answer questions of this type. The measurement and estimation of hydrogeologic parameters is treated in Sections 8.4 through 8.7. Predictions of drawdown in an aquifer under a proposed pumping scheme can be carried out for simple situations with the analytical methods presented in Section 8.3. More complex hydrogeological environments may require the application of numerical-simulation techniques, as presented in Section 8.8, or electrical-analog techniques, as presented in Section 8.9. Land subsidence is discussed in Section 8.12, and seawater intrusion in Section 8.13.

Well Yield, Aquifer Yield, and Basin Yield

The techniques of groundwater resource evaluation require an understanding of the concept of groundwater yield, and, perhaps surprisingly, this turns out to be a difficult and ambiguous term/to address. The concept is certainly pertinent, in that one of the primary objectives of most groundwater resource studies is the determination of the maximum possible pumping rates that are compatible with the hydrogeologic environment from which the water will be taken. This need for compatibility implies that yields must be viewed in terms of a balance between the benefits of groundwater pumpage and the undesirable changes that will be induced by such pumpage. The most ubiquitous change that results from pumping is lowered water levels, so in the simplest cases groundwater yield can be defined in terms of the maximum rate of pumpage that can be allowed while ensuring that water-level declines are kept within acceptable limits.

This concept of yield can be applied on several scales. If our unit of study is a single well, we can define a well yield; if our unit of study is an aquifer, we can define an aquifer yield; and if our unit of study is a groundwater basin, we can define a basin yield. Well yield can be defined as the maximum pumping rate that can be supplied by a well without lowering the water level in the well below the pump intake. Aquifer yield can be defined as the maximum rate of withdrawal that can be sustained by an aquifer without causing an unacceptable decline in the hydraulic head in the aquifer. Basin yield can be defined as the maximum rate of withdrawal that can be sustained by the complete hydrogeologic system in a groundwater basin without causing unacceptable declines in hydraulic head in the system or causing unacceptable changes to any other component of the hydrologic cycle in the basin. In light of the effects of well interference that are discussed in Section 8.3, it is clear that aquifer yield is highly dependent on the number and spacing of wells tapping an aquifer. If all the wells in a highly developed aquifer pump at a rate equal to their well yield, it is likely that the aquifer yield will be exceeded. In light of the effects of aquitard leakage and aquifer interference that are also discussed in Section 8.3, it is clear that basin yield is highly dependent on the number and spacing of exploited aquifers in a basin. If all the aquifers are pumped at a rate equal to their aquifer yield, it is likely that the basin yield will be exceeded.

These simple concepts should prove useful to the reader in the early sections of this chapter. However, the concept of basin yield deserves reconsideration in greater depth, and this is presented in Section 8.10.

8.2 Exploration for Aquifers

An aquifer is a geological formation that is capable of yielding economic quantities of water to man through wells. It must be porous, permeable, and saturated. While aquifers can take many forms within the wide variety of existing hydrogeological environments, a perusal of the permeability and porosity data of Tables 2.2 and 2.4 and consideration of the discussions of Chapter 4 make it clear that certain geological deposits are of recurring interest as aquifers. Among the most common are unconsolidated sands and gravels of alluvial, glacial, lacustrine, and deltaic origin; sedimentary rocks, especially limestones and dolomites, and sandstones and conglomerates; and porous or fractured volcanic rocks. In most cases, aquifer exploration becomes a search for one or other of these types of geological deposits. The methods of exploration can be grouped under four headings: surface geological, subsurface geological, surface geophysical, and subsurface geophysical.

Surface Geological Methods

The initial steps in a groundwater exploration program are carried out in the office rather than in the field. Much can be learned from an examination of available maps, reports, and data. There are published geologic maps on some scale for almost all of North America; there are published soils maps or surficial geology maps for most areas; and there are published hydrogeological maps for some areas. Geologic maps and reports provide the hydrogeologist with an initial indication of the rock formations in an area, together with their stratigraphic and structural interrelationships. Soils maps or surficial geology maps, together with topographic maps, provide an introduction to the distribution and genesis of the surficial unconsolidated deposits and their associated landforms. Hydrogeologic maps provide a summarized interpretation of the topographic, geologic, hydrogeologic, geochemical, and water resource data available in an area.

Airphoto interpretation is also widely used in groundwater exploration. It is usually possible to prepare maps of landforms, soils, land use, vegetation, and drainage from the airphoto coverage of an area. Each of these environmental properties leads to inferences about the natural groundwater flow systems and/or the presence of potential aquifers. Way (1973) and Mollard (1973) each provide a handbook-style treatment of airphoto-interpretation methods, and both include a large number of interpreted photos, many of which illustrate significant hydrogeological features.

However, even in areas where there is a considerable amount of published information, it is usually necessary to carry out geologic mapping in the field. In view of the importance of unconsolidated sands and gravels as potential aquifers, special attention must be paid to geomorphic landforms and to the distribution of glacial and alluvial deposits. Where sand and gravel deposits are sparse, or where these deposits are shallow and unsaturated, more detailed attention must be paid to the lithology, stratigraphy, and structure of the bedrock formations.

The methods of hydrogeologic mapping outlined in Section 6.1 are useful in determining the scale and depth of natural groundwater flow systems and in mapping the extent of their recharge and discharge areas.

Subsurface Geological Methods

It is seldom sufficient to look only at the surficial manifestations of a hydrogeological environment. It is unlikely that subsurface stratigraphic relationships will be fully revealed without direct subsurface investigation. Once again, the initial step usually involves scanning the available records. Many state and provincial governments now require that geological logs of all water wells be filed in a central bank for the use of other investigators. These data, while varying widely in quality, can often provide the hydrogeologist with considerable information on past successes and failures in a given region.

In most exploration programs, especially those for large-scale industrial or municipal water supplies, it is necessary to carry out test-drilling to better delineate subsurface conditions. Test holes provide the opportunity for geological and geophysical logging and for the coring or sampling of geological materials. Test holes can also be used to obtain water samples for chemical analysis and to indicate the elevation of the water table at a site. Test-drilling programs, together with published geological maps and available well-log records, can be interpreted in terms of the local and regional lithology, stratigraphy, and structure. Their logs can be used to prepare stratigraphic cross sections, geological fence diagrams, isopach maps of overburden thickness or formation thickness, and lithofacies maps. Hydrogeological interpretations might include water-table contours and isopachs of saturated thickness of unconfined aquifers. The results of chemical analyses of groundwater samples, when graphically displayed using the methods of Chapter 7, can provide important evidence on the natural geochemical environment as well as a direct measure of water quality.

Surface Geophysical Methods

There are two regional geophysical techniques that are used to some extent in the exploration for aquifers. These are the seismic refraction method and the electrical resistivity method. The design of geophysical surveys that utilize these approaches, and the interpretation of the resulting geophysical measurements, is a specialized branch of the earth sciences. It is not expected that a groundwater hydrologist become such a specialist, and for this reason our discussion is brief. On the other hand, it is necessary that the hydrogeologist be aware of the power and limitations of the methods. If this brief presentation fails to meet that objective, the reader is directed to a standard geophysics textbook such as Dobrin (1960), or to one of several review papers that deal specifically with geophysical applications in groundwater exploration, such as McDonald and Wantland (1961), Hobson (1967), or Lennox and Carlson (1967).

The seismic refraction method is based on the fact that elastic waves travel through different earth materials at different velocities. The denser the material, the higher the wave velocity. When elastic waves cross a geologic boundary between two formations with different elastic properties, the velocity of wave propagation changes and the wave paths are refracted according to Snell’s law. In seismic exploration, elastic waves are initiated by an energy source, usually a small explosion, at the ground surface. A set of receivers, called geophones, is set up in a line radiating outward from the energy source. Waves initiated at the surface and refracted at the critical angle by a high-velocity layer at depth will reach the more distant geophones more quickly than waves that travel directly through the low velocity surface layer. The time between the shock and the arrival of the elastic wave at a geophone is recorded on a seismograph. A set of seismograph records can be used to derive a graph of arrival time versus distance from shot point to geophone, and this, in turn, with the aid of some simple theory, can be used to calculate layer depths and their seismic velocities.

In groundwater investigations the seismic refraction method has been used to determine such features as the depth to bedrock, the presence of buried bedrock channels, the thickness of surficial fracture zones in crystalline rock, and the areal extent of potential aquifers. The interpretations are most reliable in cases where there is a simple two-layer or three-layer geological configuration in which the layers exhibit a strong contrast in seismic velocity. The velocities of the layers must increase with depth; the method cannot pick up a low-velocity layer (which might well be a porous potential aquifer) that underlies a high-velocity surface layer. The depth of penetration of the seismic method depends on the strength of the energy source. For shallow investigations (say, up to 30 m) hydrogeologists have often employed hammer seismic methods, in which the energy source is simply a hammer blow on a steel plate set on the ground surface.

The electrical resistivity of a geological formation is defined as ρ = RA/L, where R is the resistance to electrical current for a unit block of cross-sectional area A and length L. The resistivity controls the gradient in electrical potential that will be set up in a formation under the influence of an applied current. In a saturated rock or soil, the resistivity is largely dependent on the density and porosity of the material and on the salinity of the saturating fluid. In an electrical resistivity survey an electric current is passed into the ground through a pair of current electrodes and the potential drop is measured across a pair of potential electrodes. The spacing of the electrodes controls the depth of penetration. At each setup an apparent resistivity is calculated on the basis of the measured potential drop, the applied current, and the electrode spacing. Sets of measurements are taken either in the form of lateral profiling or depth profiling. In lateral profiling the electrode spacing is kept constant as electrodes are leapfrogged down a survey line. This method provides areal coverage at a given depth of penetration. It can be used to define aquifer limits or to map areal variations in groundwater salinity. In depth profiling a series of readings is taken at different electrode spacings at a single station. Apparent resistivities are plotted against electrode spacing, and stratigraphic interpretations are made by comparing the resulting curve against published theoretical curves for simple layered geometries. Depth profiling has been widely used to determine the thickness of sand and gravel aquifers that overlie bedrock. It can also be used to locate the saltwater-freshwater interface in coastal aquifers. It is often claimed that the method can “feel” the water table, but this is questionable except in very homogeneous deposits. In urban areas the method is often hampered by the presence of pipes, rails, and wires that interfere with the current fields.

Surface geophysical methods cannot replace test drilling, although by providing data that lead to a more intelligent selection of test-hole drilling, they may lead to a reduction in the amount of drilling required. Stratigraphic interpretations based on seismic or electrical resistivity measurements must be calibrated against test-hole information.

Subsurface Geophysical Methods

There is one geophysical approach that has now become a standard tool in groundwater exploration. This approach involves the logging of wells and test holes by the methods of borehole geophysics. The term encompasses all techniques in which a sensing device is lowered into a hole in order to make a record that can be interpreted in terms of the characteristics of the geologic formations and their contained fluids. The techniques of borehole geophysics were originally developed in the petroleum industry and the standard textbooks on the interpretation of geophysical logs (Pirson, 1963; Wyllie, 1963) emphasize petroleum applications. Fortunately, there are several excellent review articles (Jones and Skibitzke, 1956; Patton and Bennett, 1963; Keys, 1967, 1968) that deal specifically with the application of geophysical logging techniques to groundwater problems.

A complete borehole geophysics program as it is carried out in the petroleum industry usually includes two electric logs (spontaneous potential and resistivity), three radiation logs (natural gamma, neutron, and gamma-gamma), and a caliper log that indicates variations in hole diameter. In hydrogeological applications, emphasis is usually placed on the electric logs.

The simplest electric log is the spontaneous potential (or self-potential) log. It is obtained with the single-point electrode arrangement shown in Figure 8.1 with the current source disconnected. It provides a measure of the naturally occurring potential differences between the surface electrode and the borehole electrode. The origin of these natural electric potentials is not well understood, but they are apparently related to electrochemical interactions that take place between the borehole fluid and the in situ rock-water complex.

Figure 8.1 Single-point electrode arrangement for spontaneous potential and resistivity logging in a borehole.

The second electric log is a resistivity log. There are several electrode arrangements that can be used, but the simplest and the one most widely used in the water well industry is the single-point arrangement shown in Figure 8.1. The potential difference recorded at different depths for a given current strength leads to a log of apparent resistivity versus depth.

The two electric logs can be jointly interpreted in a qualitative sense in terms of the stratigraphic sequence in the hole. Figure 8.2 shows a pair of single-point electric logs measured in a test hole in an unconsolidated sequence of Pleistocene and Upper Cretaceous sediments in Saskatchewan. The geologic descriptions and the geologic log in the center are based on a core-sampling program.

Figure 8.2 Geologic log, electric logs, geologic description, and hydrologic description of a test hole in Saskatchewan (after Christiansen et al., 1971).

The hydrologic description of the potential aquifers at the site is based on a joint interpretation of the geologic and geophysical logs. In most common geological environments, the best water-yielding zones have the highest resistivities. Electric logs often provide the most accurate detail for the selection of well-screen settings.

Dyck et al. (1972) pointed out three disadvantages to single-point electric logs. They do not provide quantitative values of formation resistivity; they are affected by hole diameter and borehole fluid resistivity; and they have only a shallow radius of investigation. To emphasize the first point, the resistivity log on Figure 8.2 records simply the resistance measured between the two electrodes rather than an apparent resistivity. Multiple-point electric logs are more versatile. They can be used for quantitative calculations of the resistivity of formation rocks and their enclosed fluids. These calculations lie beyond the scope of this presentation. Campbell and Lehr (1973) provide a good summary of the techniques. Dyck et al. (1972) provide some sample calculations in the context of a groundwater exploration program.

Keys (1967, 1968) has suggested that radiation logs, especially the natural gamma log, may have applications to groundwater hydrology. A logging suite that might be considered complete for hydrogeological purposes would include a driller’s log (including drilling rate), a geologic log, a spontaneous potential log, a resistivity log, a natural gamma log, and a caliper log.

Drilling and Installation of Wells and Piezometers

The drilling of piezometers and wells, and their design, construction, and maintenance, is a specialized technology that rests only in part on scientific and engineering principles. There are many books (Briggs and Fiedler, 1966; Gibson and Singer, 1971; Campbell and Lehr, 1973; U.S. Environmental Protection Agency, 1973a, 1976) that provide a comprehensive treatment of water well technology. In addition, Walton (1970) presents material on the technical aspects of groundwater hydrology, and his text includes many case histories of water well installations and evaluations. Reeve (1965), Hvorslev (1951), Campbell and Lehr (1973), and Kruseman and de Ridder (1970) discuss methods of piezometer construction and installation. In this text we will limit ourselves to a brief overview of these admittedly important practical matters. Most of what follows is drawn from Campbell and Lehr (1973).

Water wells are usually classified on the basis of their method of construction. Wells may be dug by hand, driven or jetted in the form of well points, bored by an earth auger, or drilled by a drilling rig. The selection of the method of construction hinges on such questions as the purpose of the well, the hydrogeological environment, the quantity of water required, the depth and diameter envisaged, and economic factors. Dug, bored, jetted and driven wells are limited to shallow depths, unconsolidated deposits, and relatively small yields. For deeper, more productive wells in unconsolidated deposits, and for all wells in rock, drilling is the only feasible approach.

There are three main types of drilling equipment: cable tool, rotary, and reverse rotary. The cable tool drills by lifting and dropping a string of tools suspended on a cable. The bit at the bottom of the tool string rotates a few degrees between each stroke so that the cutting face of the bit strikes a different area of the hole bottom with each stroke. Drilling is periodically interrupted to bail out the cuttings. With medium- to high-capacity rigs, 40- to 60-cm-diameter holes can be drilled to depths of several hundred meters and smaller diameter holes to greater depths. The cable-tool approach is successful over a wide range of geological materials, but it is not capable of drilling as quickly or as deeply as rotary methods. With the conventional rotary method, drilling fluid is forced down the inside of a rapidly rotating drill stem and out through openings in the bit. The drilling fluid flows back to the surface, carrying the drill cuttings with it, by way of the annulus formed between the outside of the drill pipe and the hole wall. In a reverse rotary system, the direction of circulation is reversed. Reverse rotary is particularly well suited to drilling large-diameter holes in soft, unconsolidated formations.

The conventional rotary rig is generally considered to be the fastest, most convenient, and least expensive system to operate, especially in unconsolidated deposits. Penetration rates for rotary rigs depend on such mechanical factors as the weight, type, diameter, and condition of the bit, and its speed of rotation; the circulation rate of the drilling fluid and its properties; and the physical characteristics of the geological formation. In rock formations, drillability (defined as depth of penetration per revolution) is directly related to the compressive strength of the rock.

The direct rotary method is heavily dependent on its hydraulic circulation system. The most widely used drilling fluid is a suspension of bentonitic clay in water, known as drilling mud. During drilling, the mud coats the hole wall and in so doing contributes to the hole stability and prevents losses of the drilling fluid to permeable formations. When even heavy drilling mud cannot prevent the caving of hole walls, well casing must be emplaced as drilling proceeds. Caving, lost circulation, and conditions associated with the encounter of flowing artesian water constitute the most common drilling problems.

The design of a deep-cased well in an unconsolidated aquifer must include consideration of the surface housing, the casing, the pumping equipment, and the intake. Of these, it is the intake that is most often of primary concern to groundwater hydrologists. In the first half of this century it was quite common to provide access for the water to the well by a set of perforations or hand-sawn slots in the casing. It is now recognized that well yields can be significantly increased by the use of well screens. The size of the intake slots in a properly designed well screen is related to the grain-size distribution of the aquifer. Development of a screened well by pumping, surging, or backwashing draws the lines out of the aquifer, through the well screen, and up to the surface. By removing the lines from the formation in the vicinity of the well, a natural gravel pack is created around the screen that increases the efficiency of the intake. In some cases, an artificial gravel pack is emplaced to improve intake properties. Figure 8.3 shows several typical designs for wells in consolidated and unconsolidated formations.

Figure 8.3 Typical well designs for consolidated and unconsolidated formations.

The productivity of a well is often expressed in terms of the specific capacity, Cs, which is defined as Cs = Qhw, where Q is the pumping rate and Δhw is the drawdown in the well. In this equation, Δhw = Δh + ΔhL, where Δh is the drawdown in hydraulic head in the aquifer at the well screen boundary, and ΔhL is the well loss created by the turbulent flow of water through the screen and into the pump intake. Δh is calculated from the standard well-hydraulics equations developed in Section 8.3. ΔhL can be estimated by methods outlined in Walton (1970) and Campbell and Lehr (1973). In general, ΔhL \ll Δh.

8.3 The Response of ldeal Aquifers to Pumping

The exploitation of a groundwater basin leads to water-level declines that serve to limit yields. One of the primary goals of groundwater resource evaluation must therefore be the prediction of hydraulic-head drawdowns in aquifers under proposed pumping schemes. In this section, the theoretical response of idealized aquifers to pumping will be examined. We will investigate several types of aquifer configuration, but in each case the geometry will be sufficiently regular and the boundary conditions sufficiently simple to allow the development of an analytical solution to the boundary-value problem that represents the case at hand. These solutions, together with solutions to more complex boundary-value problems that describe less ideal conditions, constitute the foundation of the study of well hydraulics. This section provides an introduction to the topic, but the material covered is far from all-inclusive. There is a massive literature in the field and the committed reader is directed to Walton’s (1970) comprehensive treatment, to Hantush’s (1964) monograph, or to the excellent handbooks of Ferris et al. (1962) and Kruseman and de Ridder (1970).

Radial Flow to a Well

The theoretical analyses are based on an understanding of the physics of flow toward a well during pumping. All the necessary concepts have been introduced in Chapter 2. The distinction between confined and unconfined aquifers was explained there, as was the relation between the general concept of hydraulic head in a three-dimensional geologic system and the specific concept of the potentiometric surface on a two-dimensional, horizontal, confined aquifer. Definitions were presented for the fundamental hydrogeologic parameters: hydraulic conductivity, porosity, and compressibility; and for the derived aquifer parameters: transmissivity and storativity. It was explained there that pumping induces horizontal hydraulic gradients toward a well, and as a result hydraulic heads are decreased in the aquifer around a well during pumping. What is required now is that we take these fundamental concepts, put them into the form of a boundary-value problem that represents flow to a well in an aquifer, and examine the theoretical response.

At this point it is worth recalling from Section 2.10 that the definition of storativity invokes a one-dimensional concept of aquifer compressibility. The ∝ in Eq. (2.63) is the aquifer compressibility in the vertical direction. The analyses that follow in effect assume that changes in effective stress induced by aquifer pumping are much larger in the vertical direction than in the horizontal.

The concept of aquifer storage inherent in the storativity term also implies an instantaneous release of water from any elemental volume in the system as the head drops in that element.

Let us begin our analysis with the simplest possible aquifer configuration. Consider an aquifer that is (1) horizontal, (2) confined between impermeable formations on top and bottom, (3) infinite in horizontal extent, (4) of constant thickness, and (5) homogeneous and isotropic with respect to its hydrogeological parameters.

For the purposes of our initial analysis, let us further limit our ideal system as follows: (1) there is only a single pumping well in the aquifer, (2) the pumping rate is constant with time, (3) the well diameter is infinitesimally small, (4) the well penetrates the entire aquifer, and (5) the hydraulic head in the aquifer prior to pumping is uniform throughout the aquifer.

The partial differential equation that describes saturated flow in two horizontal dimensions in a confined aquifer with transmissivity T and storativity S was developed in Section 2.11 as Eq. (2.77):

\frac{\partial^2h}{\partial r^2} + \frac{1}{r}\frac{\partial h}{\partial r} = \frac{S}{T}\frac{\partial h}{\partial t}

Since it is clear that hydraulic-head drawdowns around a well will possess radial symmetry in our ideal system, it is advantageous to convert Eq. (2.77) into radial coordinates. The conversion is accomplished through the relation r = \sqrt{x^2 + y^2} and the equation of flow becomes (Jacob, 1950)

\frac{\partial^2h}{\partial r^2} + \frac{1}{r} \frac{\partial h}{r \partial \hspace{1mm} r} (8.1)

The mathematical region of flow, as illustrated in the plan view of Figure 8.4, is a horizontal one-dimensional line through the aquifer, from r = 0 at the well to r = ∞ at the infinite extremity.

Figure 8.4 Radial flow to a well in a horizontal confined aquifer.

The initial condition is

h(r, 0) = h0 for all r (8.2)

where h0 is the constant initial hydraulic head.

The boundary conditions assume no drawdown in hydraulic head at the infinite boundary:

h(∞, t) = h0 for all t (8.3)

And a constant pumping rate Q[L3/T] at the well:

\text{lim}_{r\rightarrow0} \left( r\frac{\partial h}{\partial r}\right)^n = \frac{Q}{2 \pi T} \hspace{1cm} \text{for} \hspace{1mm} t > 0 (8.4)

Condition (8.4) is the result of a straightforward application of Darcy’s law at the well face.

The solution h(r, t) describes the hydraulic head field at any radial distance r at any time after the start of pumping. For reasons that should be clear from a perusal of Figure 8.4, solutions are often presented in terms of the drawdown in head h0h(r, t).

The Theis Solution

Theis (1935), in what must be considered one of the fundamental breakthroughs in the development of hydrologic methodology, utilized an analogy to heat-flow theory to arrive at an analytical solution to Eq. (8.1) subject to the initial and boundary conditions of Eqs. (8.2) through (8.4). His solution, written in terms of the drawdown, is

h_0 - h(r, t) = \frac{Q}{4 \pi T} \int^\infty_u \frac {e^{-u}du}{u} (8.5)


u = \frac{r^2}{4 \pi t} (8.6)

The integral in Eq. (8.5) is well known in mathematics. It is called the exponential integral and tables of values are widely available. For the specific definition of u given by Eq. (8.6), the integral is known as the well function, W(u). With this notation, Eq. (8.5) becomes

h_0 - h = \frac{Q}{4Tt}W(u) (8.7)

Table 8.1 provides values of W(u) versus u, and Figure 8.5(a) shows the relationship W(u) versus 1/u graphically. This curve is commonly called the Theis curve.

If the aquifer properties, T and S, and the pumping rate, Q, are known, it is possible to predict the drawdown in hydraulic head in a confined aquifer at any distance r from a well at any time t after the start of pumping. It is simply necessary to calculate u from Eq. (8.6), look up the value of W(u) on Table 8.1, and calculate h0h from Eq. (8.7).

Table 8.1 Values of W(u) for Various Values of u

u 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0
× 1 0.219 0.049 0.013 0.0038 0.0011 0.00036 0.00012 0.000038 0.000012
× 10–1 1.82 1.22 0.91 0.70 0.56 0.45 0.37 0.31 0.26
× 10–2 4.04 3.35 2.96 2.68 2.47 2.30 2.15 2.03 1.92
× 10–3 6.33 5.64 5.23 4.95 4.73 4.54 4.39 4.26 4.14
× 10–4 8.63 7.94 7.53 7.25 7.02 6.84 6.69 6.55 6.44
× 10–5 10.94 10.24 9.84 9.55 9.33 9.14 8.99 8.86 8.74
× 10–6 13.24 12.55 12.14 11.85 11.63 11.45 11.29 11.16 11.04
× 10–7 15.54 14.85 14.44 14.15 13.93 13.75 13.60 13.46 13.34
× 10–8 17.84 17.15 16.74 16.46 16.23 16.05 15.90 15.76 15.65
× 10–9 20.15 19.45 19.05 18.76 18.54 18.35 18.20 18.07 17.95
× 10–10 22.45 21.76 21.35 21.06 20.84 20.66 20.50 20.37 20.25
× 10–11 24.75 24.06 23.65 23.36 23.14 22.96 22.81 22.67 22.55
× 10–12 27.05 26.36 25.96 25.67 25.44 25.26 25.11 24.97 24.86
× 10–13 29.36 28.66 28.26 27.97 27.75 27.56 27.41 27.28 27.16
× 10–14 31.66 30.97 30.56 30.27 30.05 29.87 29.71 29.58 29.46
× 10–15 33.96 33.27 32.86 32.58 32.35 32.17 32.02 31.88 31.76
SOURCE: Wenzel, 1942.

Figure 8.5(b) shows a calculated plot of h0h versus t for the specific set of parameters noted on the figure. A set of field measurements of drawdown versus time measured in a piezometer that is set in an ideal confined aquifer with these properties would show this type of record.

The shape of the function h0h versus t, when plotted on log-log paper as in Figure 8.5(b), has the same form as the plot of W(u) versus 1/u shown in Figure 8.5(a). This is a direct consequence of the relations embodied in Eqs. (8.6) and (8.7), where it can be seen that h0h and W(u), and t and 1/u, are related to one another through a constant term.

It is also possible to calculate values of h0h at various values of r at a given time t. Such a calculation leads to a plot of the cone of depression (or drawdown cone) in the potentiometric surface around a pumping well. Figure 8.4 provides a schematic example. The steepening of the slope of the cone near the well is reflected in the solution, Eq. (8.7). The physical explanation is clear if one carries out the simple flow-net construction shown in the plan view of Figure 8.4 and then carries the hydraulic head values down onto the section.

For a given aquifer the cone of depression increases in depth and extent with increasing time. Drawdown at any point at a given time is directly proportional to the pumping rate and inversely proportional to aquifer transmissivity and aquifer storativity. As shown in Figure 8.6, aquifers of low transmissivity develop tight, deep drawdown cones, whereas aquifers of high transmissivity develop shallow cones of wide extent. Transmissivity exerts a greater influence on drawdown than does storativity.

In that geologic configurations are seldom as ideal as that outlined above, the time-drawdown response of aquifers under pumpage often deviates from the Theis solution shown in Figure 8.5.

Figure 8.5 (a) Theoretical curve of W(u) versus 1/u. (b) Calculated curve of h0h versus t.

Theis solution shown in Figure 8.5. We will now turn to some of the theoretical response curves that arise in less ideal situations. Specifically, we will look at (1) leaky aquifers, (2) unconfined aquifers, (3) multiple-well systems, (4) stepped pumping rates, (5) bounded aquifers, and (6) partially penetrating wells.

Figure 8.6 Comparison of drawdown cones at a given time for aquifers of (a) low transmissivity; (b) high transmissivity; (c) low storativity; (d) high storativity.

Leaky Aquifers

The assumption inherent in the Theis solution that geologic formations overlying and underlying a confined aquifer are completely impermeable is seldom satisfied. Even when production wells are screened only in a single aquifer, it is quite usual for the aquifer to receive a significant inflow from adjacent beds. Such an aquifer is called a leaky aquifer, although in reality it is the aquitard that is leaky. The aquifer is often just one part of a multiple-aquifer system in which a succession of aquifers are separated by intervening low-permeability aquitards. For the purposes of this section, however, it is sufficient for us to consider the three-layer case shown in Figure 8.7. Two aquifers of thickness b1 and b2 and horizontal hydraulic conductivities K1 and K2 are separated by an aquitard of thickness b’ and vertical hydraulic conductivity K’. The specific storage values in the aquifers are SS1 and SS2, while that in the aquitard is S’S.

Since a rigorous approach to flow in multiple-aquifer systems involves boundary conditions that make the problem intractable analytically, it has been customary to simplify the mathematics by assuming that flow is essentially horizontal in the aquifers and vertical in the aquitards. Neuman and Witherspoon (1969a) report that the errors introduced by this assumption are less than 5% when the conductivities of the aquifers are more than 2 orders of magnitude greater than that of the aquitard.

Figure 8.7 Schematic diagram of a two-aquifer “leaky” system. Recall that T = Kb and S = S1b.

The development of leaky-aquifer theory has taken place in two distinct sets of papers. The first, by Hantush and Jacob (1955) and Hantush (1956, 1960), provided the original differentiation between the Theis response and that for leaky aquifers. The second, by Neuman and Witherspoon (1969a, 1969b, 1972) evaluated the significance of the assumptions inherent in the earlier work and provided more generalized solutions.

The analytical solution of Hantush and Jacob (1955) can be couched in the same form as the Theis solution [Eq. (8.7)] but with a more complicated well function. In fact, Hantush and Jacob developed two analytical solutions, one valid only for small t and one valid only for large t, and then interpolated between the two solutions to obtain the complete response curve. Their solution is presented in terms of dimensionless parameter, r/B, defined by the relation

\frac{r}{B} = r \sqrt{\frac{K'}{K_1b_1b'}} (8.8)

In analogy with Eq. (8.7), we can write their solution as

h_0 - h = \frac{Q}{4Tt}W(u,r/B) (8.9)

where W(r/B) is known as the leaky well function.

Hantush (1956) tabulated the values of W(r/B). Figure 8.8 is a plot of this function against 1/u. If the aquitard is impermeable, then K’ = 0, and from Eq. (8.8), r/B = 0. In this case, as shown graphically in Figure 8.8, the Hantush-Jacob solution reduces to the Theis solution.

If T1 (= K1b1) and S1 (= S1b1) are known for the aquifer and K’ and b’ are known for the aquitard, then the drawdown in the hydraulic head in the pumped aquifer for any pumpage Q at any radial distance r at any time t can be calculated from Eq. (8.9), after first calculating u for the pumped aquifer from Eq. (8.6), r/B from Eq. (8.8), and W(u, r/B) from Figure 8.9.

Figure 8.8 Theoretical curves of W(u,r/B) versus 1/u for a leaky aquifer (after Walton, 1960).
Figure 8.9 Theoretical curves of W(u,r/B11, r/B21, β11, β21) versus 1/u for a leaky two-aquifer system (after Neuman and Witherspoon, 1969a).

The original Hantush and Jacob (1955) solution was developed on the basis of two very restrictive assumptions. They assumed that the hydraulic head in the unpumped aquifer remains constant during the removal of water from the pumped aquifer and that the rate of leakage into the pumped aquifer is proportional to the hydraulic gradient a across the leaky aquitard. The first assumption implies that the unpumped aquifer has an unlimited capacity to provide water for delivery through the aquitard to the pumped aquifer. The second assumption completely ignores the effects of the storage capacity of the aquitard on the transient solution (i.e., it is assumed that S’S = 0).

In a later paper, Hantush (1960) presented a modified solution in which consideration was given to the effects of storage in the aquitard. More recently, Neuman and Witherspoon (1969a, 1969b) presented a complete solution that includes consideration of both release of water from storage in the aquitard and head drawdowns in the unpumped aquifer. Their solutions require the calculation of four dimensionless parameters, which, with reference to Figure 8.7, are defined as follows:

\frac{r}{B_{11}} = r \sqrt{\frac{K'}{K_1b_1b'}}

\frac{r}{B_{21}} = r \sqrt{\frac{K'}{K_2b_2b'}}

\beta_{11} = \frac{r}{4b_1} \sqrt{\frac{K'S'_S}{K_1S_{S_1}}}

\beta_{21} = \frac{r}{4b_2} \sqrt{\frac{K'S'_S}{K_2S_{S_2}}}

Neuman and Witherspoon’s solutions provide the drawdown in both aquifers as a function of radial distance from the well, and in the aquitard as a function of both radial distance and elevation above the base of the aquitard. Their solutions can be described in a schematic sense by the relation

h_0 = h(r,z,t) = \frac{Q}{4\pi T}W(u, r/B_{11}, r/B_{21}, \beta_{11}, \beta_{21}) (8.11)

Tabulation of this well function would require many pages of tables, but an indication of the nature of the solutions can be seen from Figure 8.9, which presents the theoretical response curves for the pumped aquifer, and at three elevations in the aquitard, for a specific set of r/B and β values. The Theis solution is shown on the diagram for comparative purposes.

Because of its simplicity, and despite the inherent dangers of using a simple model for a complex system, the r/B solution embodied in Figure 8.8 is widely used for the prediction of drawdowns in leaky-aquifer systems. Figure 8.10 shows an h0h versus t plot for a specific case as calculated from Eq. (8.9) with the aid of Figure 8.8.

Figure 8.10 Calculated curve of h0h versus t for a leaky aquifer, based on Hantush-Jacob theory.

The drawdown reaches a constant level after about 5 × 103 seconds. From this point on, the r/B solution indicates that steady-state conditions hold throughout the system, with the infinite storage capacity assumed to exist in the upper aquifer feeding water through the aquitard toward the well. If the overlaying aquitard were impermeable rather than leaky, the response would follow the dotted line. As one would expect, drawdowns in leaky aquifers are less than those in non-leaky aquifiers, as there is now an additional source of water over and above that which can be supplied by the aquifer itself. Predictions based on the Theis equation therefore provide a conservative estimate for leaky systems; that is, they over-predict the drawdown, or, put another way, are unlikely to reach the values predicted by the Theis equation for a given pumping scheme in a multiaquifer system.

Unconfined Aquifers

When water is pumped from a confined aquifer, the pumpage induces hydraulic gradients toward the well that create drawdowns in the potentiometric surface. The water produced by the well arises from two mechanisms: expansion of the water in the aquifer under reduced fluid pressures, and compaction of the aquifer under increased effective stresses (Section 2.10). There is no dewatering of the geologic system. The flow system in the aquifer during pumping involves only horizontal gradients toward the well; there are no vertical components of flow. When water is pumped from an unconfined aquifer, on the other hand, the hydraulic gradients that are induced by the pumpage create a drawdown cone in the water table itself and there are vertical components of flow (Figure 8.11). The water produced by the well arises from the two mechanisms responsible for confined delivery plus the actual dewatering of the unconfined aquifer.

There are essentially three approaches that can be used to predict the growth of unconfined drawdown cones in time and space. The first, which might be termed the complete analysis, recognizes that the unconfined well-hydraulics problem (Figure 8.11) involves a saturated-unsaturated flow system in which water-table drawdowns are accompanied by changes in the unsaturated moisture contents above the water table (such as those shown in Figure 2.23).

Figure 8.11 Radial flow to a well in an unconfined aquifer.

The complete analysis requires the solution of a boundary-value problem that includes both the saturated and unsaturated zones. An analytical solution for this complete case was presented by Kroszynski and Dagan (1975) and several numerical mathematical models have been prepared (Taylor and Luthin, 1969; Cooley, 1971; Brutsaert et al., 1971). The general conclusion of these studies is that the position of the water table during pumpage is not substantially affected by the nature of the unsaturated flow above the water table. In other words, while it is conceptually more appealing to carry out a complete saturated-unsaturated analysis, there is little practical advantage to be gained, and since unsaturated soil properties are extremely difficult to measure in situ, the complete analysis is seldom used.

The second approach, which is by far the simplest, is to use the same equation as for a confined aquifer [Eq. (8.7)] but with the argument of the well function [Eq. (8.6)] defined in terms of the specific yield Sy, rather than the storativity S. The transmissivity T must be defined as T = Kb, where b is the initial saturated thickness. Jacob (1950) has shown that this approach leads to predicted drawdowns that are very nearly correct as long as the drawdown is small in comparison with the saturated thickness. The method in effect relies on the Dupuit assumptions (Section 5.5) and fails when vertical gradients become significant.

The third approach, and the one most widely used in practice, is based on the concept of delayed water-table response. This approach was pioneered by Boulton (1954, 1955, 1963) and has been significantly advanced by Neuman (1972, 1973b, 1975a). It can be observed that water-level drawdowns in piezometers adjacent to pumping wells in unconfined aquifers tend to decline at a slower rate than that predicted by the Theis solution. In fact, there are three distinct segments that can be recognized in time-drawdown curves under water-table conditions. During the first segment, which covers only a short period after the start of pumping, an unconfined aquifer reacts in the same way as does a confined aquifer. Water is released instantaneously from storage by the compaction of the aquifer and by the expansion of the water. During the second segment, the effects of gravity drainage are felt. There is a decrease in the slope of the time-drawdown curve relative to the Theis curve because the water delivered to the well by the dewatering that accompanies the falling water table is greater than that which would be delivered by an equal decline in a confined potentiometric surface. In the third segment, which occurs at later times, time-drawdown data once again tend to conform to a Theis-type curve.

Boulton (1963) produced a semiempirical mathematical solution that reproduces all three segments of the time-drawdown curve in an unconfined aquifer. His solution, although useful in practice, required the definition of an empirical delay index that was not related clearly to any physical phenomenon. In recent years there has been a considerable amount of research (Neuman, 1972; Streltsova, 1972; Gambolati, 1976) directed at uncovering the physical processes responsible for delayed response in unconfined aquifers. It is now clear that the delay index is not an aquifer constant, as Boulton had originally assumed. It is related to the vertical components of flow that are induced in the flow system and it is apparently a function of the radius r and perhaps the time t.

The solution of Neuman (1972, 1973b, 1975a) also reproduces all three segments of the time-drawdown curve and it does not require the definition of any empirical constants. Neuman’s method recognizes the existence of vertical new components, and the general solution for the drawdown, h0h, is a function of both r and z, as defined in Figure 8.11. His general solution can be reduced to one that is a function of r alone if an average drawdown is considered. His complex analytical solution can be represented in simplified form as

h_0 - h = \frac{Q}{4\pi T}W(u_A, u_B, \eta) (8.12)

where W(uA, uB, η) is known as the unconfined well function and η = r2/b2, Figure 8.12 is a plot of this function for various values of η. The type A curves that grow out of the left-hand Theis curve of Figure 8.12, and that are followed at early time, are given by

h_0 - h = \frac{Q}{4\pi T}W(u_A, \eta) (8.13)


u_A = \frac{r^2S}{4Tt}

and S is the elastic storativity responsible for the instantaneous release of water to the well. The type B curves that are asymptotic to the right-hand Theis curve of Figure 8.12, and that are followed at later time, are given by

h_0 - h = \frac{Q}{4\pi T}W(u_B, \eta) (8.14)


u_B = \frac{r^2S_y}{4Tt}

and Sy is the specific yield that is responsible or the delayed release of water to the well.

Figure 8.12 Theoretical curves of W(uA, uB, η) versus 1/uA and 1/uB for an unconfined aquifer (after Neuman, 1975a).

For an anisotropic aquifer with horizontal hydraulic conductivity Kr and vertical hydraulic conductivity Kz, the parameter η is given by

\eta = \frac{r^2K_z}{b^2K_r} (8.15)

If the aquifer is isotropic, Kz = Kr, and η = r2/b2. The transmissivity T is defined as T = Krb. Equations (8.12) through (8.15) are only valid if Sy \gg S and h0h \ll b.

The prediction of the average drawdown at any radial distance r from a pumping well at any time t can be obtained from Eqs. (8.13) through (8.15) given Q, S, Sy, Kr, Kz, and b.

Multiple Well Systems, Stepped Pumping Rates, Well Recovery, and Partial Penetration

The drawdown in hydraulic head at any point in a confined aquifer in which more than one well is pumping is equal to the sum of the drawdowns that would arise from each of the wells independently. Figure 8.13 schematically displays the drawdown h0h at a point B situated between two pumping wells with pumping rates, Q1 = Q2. If Q1Q2, the symmetry of the diagram about the plane AA’ would be lost but the principles remain the same.

Figure 8.13 Drawdown in the potentiometric surface of a confined aquifer being pumped by two wells with Q1 = Q2.

For a system of n wells pumping at rates Q1, Q2, . . . Qn, the arithmetic summation of the Theis solutions leads to the following predictive equation for the drawdown at a point whose radial distance from each well is given by r1, r2, . . . rn,

h_0 - h = \frac{Q_1}{4 \pi T}W(u_1) + \frac{Q_2}{4 \pi T}W(u_2) + ... + \frac{Q_n}{4 \pi T}W(u_n) (8.16)


u_i = \frac{r_2S}{4Tt_t} \hspace{1cm} i = 1, 2, ..., n

and tt, is the time since pumping started at the well whose discharge is Qi.

The summation of component drawdowns outlined above is an application of the principle of superposition of solutions. This approach is valid because the equation of flow [Eq. (8.1)] for transient flow in a confined aquifer is linear (i.e., there are no cross terms of the form ∂h/∂r · ∂h/∂t. Another application of the principle of superposition is in the case of a single well that is pumped at an initial rate Q0 and then increased to the rates Q1, Q2, . . . Qm in a stepwise fashion by the additions ΔQ1, ΔQ2, . . . ΔQm. Drawdown at a radial distance r from the pumping well is given by

h_0 - h = \frac{Q_0}{4\pi T}W(u_0) + \frac{\Delta Q_1}{4\pi T}W(u_1) + ... + \frac{\Delta Q_n}{4\pi T}W(u_m) (8.17)


u_j = \frac{r^2S}{4T_j} \hspace{1cm} j = 0, 1, 2, ..., m

and tj, is the time since the start of the pumping rate Qj.

A third application of the superposition principle is in the recovery of a well after pumping has stopped. If t is the time since the start of pumping and t’
is the time since shutdown, then the drawdown at a radial distance r from the well is given by

h_0 - h = \frac{Q}{4\pi T}[W(u_1)-W(u_2)] (8.18)


u_1 = \frac{r^2S}{4Tt} \hspace {1cm} u_2 = \frac{r^2S}{4Tt'}

Figure 8.14 schematically displays the drawdowns that occur during the pumping period and the residual drawdowns that remain during the recovery period.

Figure 8.14 Schematic diagram of the recovery in hydraulic head in an aquifer after pumping is stopped.

It is not always possible, or necessarily desirable, to design a well that fully penetrates the aquifer under development. This is particularly true for unconfined aquifers, but may also be the case for thick confined aquifers. Even for wells that are fully penetrating, screens may be set over only a portion of the aquifer thickness.

Partial penetration creates vertical flow gradients in the vicinity of the well that render the predictive solutions developed for full penetration inaccurate. Hantush (1962) presented adaptations to the Theis solution for partially penetrating wells, and Hantush (1964) reviewed these solutions for both confined and leaky-confined aquifers. Dagan (1967), Kipp (1973), and Neuman (1974) considered the effects of partial penetration in unconfined aquifers.

Bounded Aquifers

When a confined aquifer is bounded on one side by a straight-line impermeable boundary, drawdowns due to pumping will be greater near the boundary [Figure 8.15(a)] than those that would be predicted on the basis of the Theis equation for an aquifer of infinite areal extent. In order to predict head drawdowns in such systems, the method of images, which is widely used in heat-flow theory, has been adapted for application in the groundwater milieu (Ferris et al., 1962).

Figure 8.15 (a) Drawdown in the potentiometric surface of a confined aquifer bounded by an impermeable boundary; (b) equivalent system of infinite extent; (c) plan view.

With this approach, the real bounded system is replaced for the purposes of analysis by an imaginary system of infinite areal extent [Figure 8.15(b)]. In this system there are two wells pumping: the real well on the left and an image well on the right. The image well pumps at a rate, Q, equal to the real well and is located at an equal distance, x1, from the boundary. If we sum the two component drawdowns in the infinite system (in identical fashion to the two-well case shown in Figure 8.13), it becomes clear that this pumping geometry creates an imaginary impermeable boundary (i.e., a boundary across which there is no flow) in the infinite system at the exact position of the real impermeable boundary in the bounded system. With reference to Figure 8.15(c), the drawdown in an aquifer bounded by an impermeable boundary is given by

h_0 - h = \frac{Q}{4 \pi T}[W(u_r) - W(u_i)] (8.19)


u_r = \frac{r^2_rS}{4Tt} \hspace{1cm} and \hspace{1cm} u_i=\frac{r^2_iS}{4Tt}

One can use the same approach to predict the decreased drawdowns that occur in a confined aquifer in the vicinity of a constant-head boundary, such as would be produced by the slightly unrealistic case of a fully penetrating stream [Figure 8.15(d)]. For this case, the imaginary infinite system [Figure 8.15(e)] includes the discharging real well and a recharging image well. The summation of the cone of depression from the pumping well and the cone of impression from the recharge well leads to an expression for the drawdown in an aquifer bounded by a constant head boundary:

h_0 - h = \frac{Q}{4 \pi T}[W(u_r) - W(u_i)] (8.20)

where ur and ui are as defined in connection with Eq. (8.19).

It is possible to use the image well approach to provide predictions of drawdown in systems with more than one boundary. Ferris et al. (1962) discuss several geometric configurations. One of the more realistic (Figure 8.16) applies to a pumping well in a confined alluvial aquifer in a more-or-less straight river valley. For this case, the imaginary infinite system must include the real pumping well R, an image well I1 equidistant from the left-hand impermeable boundary, and an image well I2 equidistant from the right-hand impermeable boundary. These image wells themselves give birth to the need for further image wells. For example, I3 reflects the effect of I2 across the left-hand boundary, and I4 reflects the effect of I1 across the right-hand boundary. The result is a sequence of imaginary pumping wells stretching to infinity in each direction. The drawdown at point P in Figure 8.16 is the sum of the effects of this infinite array of wells. In practice, image wells need only be added until the most remote pair produces a negligible effect on water-level response (Bostock, 1971).

Figure 8.16 Image-well system for pumpage from a confined aquifer in a river valley bounded by impermeable boundaries.

The Response of Ideal Aquitards

The most common geological occurrence of exploitable confined aquifers is in sedimentary systems of interbedded aquifers and aquitards. In many cases the aquitards are much thicker than the aquifers and although their permeabilities are low, their storage capacities can be very high. In the very early pumping history of a production well, most of the water comes from the depressurization of the aquifer in which the well is completed. As time proceeds the leakage properties of the aquitards are brought into play and at later times the majority of the water being produced by the well is aquitard leakage. In many aquifer-aquitard systems, the aquitards provide the water and the aquifers transmit it to the wells. It is thus of considerable interest to be able to predict the response of aquitards as well as aquifers.

In the earlier discussion of leaky aquifers, two theories were introduced: the Hantush-Jacob theory, which utilizes the W(u, r/B) curves of Figure 8.8, and the Neuman-Witherspoon theory, which utilizes the W(u, r/B11, r/B21, β11, β21) curves of Figure 8.9. In that the Hantush-Jacob theory does not include the storage properties of the aquitard, it is not suitable for the prediction of aquitard response The Neuman-Witherpoon solution, in the form of Eq. (8.11) can be used to predict the hydraulic head h(r, z, t) at any elevation z in the aquitard (Figure 8.7) at any time t, at any radial distance r, from the well. In many cases, however, it may be quite satisfactory to use a simpler approach. If the hydraulic conductivity of the aquitards is at least 2 orders of magnitude less than the hydraulic conductivity in the aquifers, it can be assumed that flow in the aquifers is horizontal and leakage in the aquitards is vertical. If one can predict, or has measurements of, h(r, t) at some point in an aquifer, one an often predict the hydraulic head h(z, t) at an overlying point in the aquitard by the application of a one-dimensional flow theory, developed by Karl Terzaghi, the founder of modern soil mechanics.

Consider an aquitard of thickness b’ (Figure 8.17) sandwiched between two producing aquifers. If the initial condition is a constant hydraulic head h = h0 in the aquitard, and if the drawdowns in hydraulic head in the adjacent aquifers can be represented as an instantaneous step function Δh the system can be represented by the following one-dimensional boundary-value problem.

Figure 8.17 Response of an ideal aquitard to a step drawdown in head in the two adjacent aquifers.

From Eq. (2.76), the one-dimensional form of the flow equation is

\frac{\partial^2h}{\partial z^2} = \frac{\rho g(\alpha ' + n'\beta)}{K'} \frac{\partial h}{\partial t} (8.21)

where the primed parameters are the aquitard properties. The initial condition is

h(z,0) = h_0

and the boundary conditions are

h(0,t) = h_0 - \Delta h

h(b', t) = h_0 - \Delta h

Terzaghi (1925) provided an analytical solution to this boundary-value problem. He noted that for clays n’β
\ll α’ in Eq. (8.21). He grouped the remaining aquitard parameters into a single parameter cv, known as the coefficient of consolidation and defined as

c_v = \frac{K'}{\rho g \alpha '} (8.22)

He further defined the dimensionless time factor, Tf, as

T_f = \frac{4c_vt}{(b')^2} (8.23)

Given the aquitard parameter cv and the geometric parameter b’, one can calculate Tf, for any time t.

Figure 8.17 is a graphical presentation of Terzaghi’s solution h(z, Tf). It allows the prediction of the hydraulic head at any elevation z at any time t in an aquitard sandwiched between two producing aquifers, as long as the drop in hydraulic head Δh can be estimated in the aquifers. It is also possible to interpret this solution for an aquitard that drains to only one aquifer. For example, if the lower boundary of the aquitard on the inset to Figure 8.17 is impermeable, only the upper half of the curves shown in the figure are used for the prediction of h(z, t). The z = 0 line passes through the center of the figure, and the parameters cy and Tf are defined as above. Wolff (1970) has described a case history that utilizes the concepts of one-dimensional aquitard response.

Predictions of aquitard response, and the inverse application of this theory to determine aquitard parameters, as discussed in Section 8.6, are also important in assessing contaminant migration (Chapter 9) and land subsidence (Section 8.12).

The Real World

Each of the analytical solutions presented in this section describes the response to pumping in a very idealized representation of actual aquifer configurations. In the real world, aquifers are heterogeneous and anisotropic; they usually vary in thickness; and they certainly do not extend to infinity. Where they are bounded, it is not by straight-line boundaries that provide perfect confinement. In the real world, aquifers are created by complex geologic processes that lead to irregular stratigraphy, interfingering of strata, and pinchouts and trendouts of both aquifers and aquitards. The predictions that can be carried out with the analytical expressions presented in this section must be viewed as best estimates. They have greater worth the more closely the actual hydrogeological environment approaches the idealized configuration.

In general, well-hydraulics equations are most applicable when the unit of study is a well or well held. They are less applicable on a larger scale, where the unit of study is an entire aquifer or a complete groundwater basin. Short-term yields around wells are very dependent on aquifer properties and well-field geometry, both of which are emphasized in the well-hydraulics equations. Long-term yields on an aquifer scale are more often controlled by the nature of the boundaries. Aquifer studies on the larger scale are usually carried out with the aid of models based on numerical simulation or electric-analog techniques. These approaches are discussed in Sections 8.8 and 8.9.

The predictive formulas developed in this section and the simulation techniques described in later sections allow one to calculate the drawdowns in hydraulic head that will occur in an aquifer in response to groundwater development through wells. They require as input either the three basic hydrogeological parameters: hydraulic conductivity, K, porosity, n, and compressibility, α; or the two derived aquifer parameters: transmissivity, T, and storativity, S. There is a wide variety of techniques that can be used to measure these parameters. In the next section, we will discuss laboratory tests; in Section 8.5, piezometer tests; and in Section 8.6, pumping tests. In Section 8.7, we will examine some estimation techniques, and in Section 8.8, the determination of aquifer parameters by inverse simulation. The formulas presented in this section are the basis for the pumping-test approach that is described in Section 8.6.

8.4 Measurement of Parameters: Laboratory Tests

The laboratory tests described in this section can be considered as providing point values of the basic hydrogeologic parameters. They are carried out on small samples that are collected during test-drilling programs or during the mapping of surficial deposits. If the samples are undisturbed core samples, the measured values should be representative of the in situ point values. For sands and gravels, even disturbed samples may yield useful values. We will describe testing methods for the determination of hydraulic conductivity, porosity, and compressibility in the saturated state; and we will provide references for the determination of the characteristic curves relating moisture content, pressure head, and hydraulic conductivity in the unsaturated state. We will emphasize principles; for a more complete description of each testing apparatus and more detailed directions on laboratory procedures, the reader is directed to the soil-testing manual by Lambe (1951), the permeability handbook of the American Society of Testing Materials (1967), or the pertinent articles in the compendium of soil analysis methods edited by Black (1965). Our discussions relate more to soils than to rocks, but the principles of measurement are the same. The rock mechanics text by Jaeger (1972) discusses rock-testing procedures.

Hydraulic Conductivity

The hydraulic conductivity, K, was defined in Section 2.1, and its relationship to the permeability, k, was explored in Section 2.3. The saturated hydraulic conductivity of a soil sample can be measured with two types of laboratory apparatus. The first type, known as a constant-head permeameter, is shown in Figure 8.18(a); the second type, a falling-head permeameter, is shown in Figure 8.18(b).

Figure 8.18 (a) Constant-head permeater; (b) falling-head permeater (after Todd, 1959).

In a constant-head test, a soil sample of length L and cross-sectional area A is enclosed between two porous plates in a cylindrical tube, and a constant-head differential H is set up across the sample. A simple application of Darcy’s law leads to the expression

K = \frac{QL}{AH} (8.24)

where Q is the steady volumetric discharge through the system. It is important that no air become entrapped in the system, and for this reason it is wise to use deaired water. If disturbed samples are being tested in the permeameter, they should be carefully saturated from below as they are emplaced.

In a falling-head test [Figure 8.18(b)], the head, as measured in a tube of cross-sectional area a, is allowed to fall from H0 to H1 during time t. The hydraulic conductivity is calculated from

K = \frac{aL}{At} \text{ln} \left( \frac{H_0}{H_1} \right) (8.25)

This equation can be derived (Todd, 1959) from the simple boundary-value problem that describes one-dimensional transient flow across the soil sample. In order that the head decline be easily measurable in a finite time period, it is necessary to choose the standpipe diameter with regard to the soil being tested. Lambe (1951) suggests that for a coarse sand a standpipe whose diameter is approximately equal to that of the permeameter is usually satisfactory, whereas a fine silt may necessitate a standpipe whose diameter is one-tenth the permeameter diameter. Lambe also suggests that the point \sqrt{H_0 H_1} be marked on the standpipe. If the time required for the head decline from H0 to \sqrt{H_0 H_1} is not equal to that for the decline from \sqrt{H_0 H_1} to H1, the test has not functioned correctly and a check should be made for leaks or entrapped air.

Klute (1965a) notes that the constant-head system is best suited to samples with conductivities greater than 0.01 cm/min while the falling-head system is best suited to samples with lower conductivity. He also notes that elaborate, painstaking measurements are not generally required for conductivity determinations on field samples. The variability among samples is usually large enough that precise determination of the conductivity of a given sample is not warranted.

For clayey materials the hydraulic conductivity is commonly determined from a consolidation test, which is described in the subsection on compressibility below.


In principle, the porosity, n, as defined in Section 2.5, would be most easily measured by saturating a sample, measuring its volume, VT, weighing it and then oven drying it to constant weight at 105°C. The weight of water removed could be converted to a volume, knowing the density of water. This volume is equivalent to the volume of the void space, Vv; and the porosity could be calculated from n = Vv/VT.

In practice, it is quite difficult to exactly and completely saturate a sample of given volume. It is more usual (Vomocil, 1965) to make use of the relationship

n = 1 - \frac{\rho_b}{\rho_s} (8.26)

which can be developed by simple arithmetic manipulation of the basic definition of porosity. In Eq. (8.26), ρb is the bulk mass density of the sample and ρs is the particle mass density. The bulk density is the oven-dried mass of the sample divided by its held volume. The particle density is the oven-dried mass divided by the volume of the solid particles, as determined by a water-displacement test. In cases where great accuracy is not required, ρs = 2.65 g/cm3 can be assumed for most mineral soils.


The compressibility of a porous medium was defined in Section 2.9 with the aid of Figure 2.19. It is a measure of the relative volumetric reduction that will take place in a soil under an increased effective stress. Compressibility is measured in a consolidation apparatus of the kind commonly used by soils engineers. In this test, a soil sample is placed in a loading cell of the type shown schematically in Figure 2.19(a). A load L is applied to the cell, creating a stress σ, where σ = L/A, A being the cross-sectional area of the sample. If the soil sample is saturated and the fluid pressure on the boundaries of the sample is atmospheric (i.e., the sample is free-draining), the effective stress, σe, which leads to consolidation of the sample, is equal to the applied stress, σ.

The reduction in sample thickness, b, is measured after equilibrium is achieved at each of several loading increments, and the results are converted into a graph of void ratio, e, versus effective stress, σe, as shown in Figure 2.19(b). The compressibility, α, is determined from the slope of such a plot by

\alpha = \frac{de/(1 + e_0)}{d \sigma_e} (8.27)

where e0 is the initial void ratio prior to loading. As noted in Section 2.9, α is a function of the applied stress and it is dependent on the previous loading history.

Lambe (1951) describes the details of the testing procedure. The most common loading method is a lever system on which weights of known magnitude are hung. There are two types of loading cell in common use. In the fixed-ring container [Figure 8.19(a)], all the sample movement relative to the container is downward. In the floating-ring container [Figure 8.19(b)], compression occurs toward the middle from both top and bottom. In the floating-ring container, the effect of friction between the container wall and the soil specimen is smaller than in the fixed-ring container. In practice, it is difficult to determine the magnitude of the friction in any case, and because its effect is thought to be minor, it is normally neglected. Cohesionless sands are usually tested as disturbed samples. Cohesive clays must be carefully trimmed to fit the consolidometer ring.

Figure 8.19 (a) Fixed-ring consolidometer; (b) floating-ring consolidometer (after Lambe, 1951).

In soil mechanics terminology, the slope of the eσe curve is called the coefficient of compressibility, av. The relationship between av and α is easily seen to be

a_v = \frac{-de}{d\sigma_e} = (1+e_0)\alpha (8.28)

More commonly, soils engineers plot the void ratio, e, against the logarithm of σe. When plotted in this manner, there is usually a significant portion of the curve that is a straight line. The slope of this line is called the compression index, Cc, where

C_c = \frac{-de}{d(\log \sigma_e)} (8.29)

In most civil engineering applications the rate of consolidation is just as important as the amount of consolidation. This rate is dependent both on the compressibility, α, and the hydraulic conductivity, K. As noted in connection with Eq. (8.22), soils engineers utilize a grouped parameter known as the coefficient of consolidation, Cv, with is defined as

V_v = \frac{K}{\rho g \alpha} (8.30)

At each loading level in a consolidation test, the sample undergoes a transient drainage process (fast for sands, slow for clays) that controls the rate of consolidation of the sample. If the rate of decline in sample thickness is recorded for each loading increment, such measurements can be used in the manner described by Lambe (1951) to determine the coefficient of consolidation, Cv, and the hydraulic conductivity, K, of the soil.

In Section 8.12, we will further examine the mechanism of one-dimensional consolidation in connection with the analysis of land subsidence.

Unsaturated Characteristic Curves

The characteristic curves, K(ψ) and theta(ψ), that relate the moisture content, θ, and the hydraulic conductivity, K, to the pressure head, ψ, in unsaturated soils were described in Section 2.6. Figure 2.13 provided a visual example of the hysteretic relationships that are commonly observed. The methods used for the laboratory determination of these curves have been developed exclusively by soil scientists. It is not within the scope of this text to outline the wide variety of sophisticated laboratory instrumentation that is available. Rather, the reader is directed to the soil science literature, in particular to the review articles by L. A. Richards (1965), Klute (1965b), Klute (1965c), and Bouwer and Jackson (1974).

8.5 Measurement of Parameters: Piezometer Tests

It is possible to determine in situ hydraulic conductivity values by means of tests carried out in a single piezometer. We will look at two such tests, one suitable for point piezometers that are open only over a short interval at their base, and one suitable for screened or slotted piezometers that are open over the entire thickness of a confined aquifer. Both tests are initiated by causing an instantaneous change in the water level in a piezometer through a sudden introduction or removal of a known volume of water. The recovery of the water level with time is then observed. When water is removed, the tests are often called bail tests; when it is added, they are known as slug tests. It is also possible to create the same effect by suddenly introducing or removing a solid cylinder of known volume.

The method of interpreting the water level versus time date that arise from bail tests or slug tests depends on which of the two test configurations is felt to be most representative. The method of Hvorslev (l951) is for a point piezometer, while that of Cooper et al. (l967) is for a confined aquifer. We will now describe each in turn.

The simplest interpretation of piezometer-recovery data is that of Hvorslev (1951). His initial analysis assumed a homogeneous, isotropic, infinite medium in which both soil and water are incompressible. With reference to the bail test of Figure 8.20(a), Hvorslev reasoned that the rate of inflow, q, at the piezometer tip at any time t is proportional to the hydraulic conductivity, K, of the soil and to the unrecovered head difference, H – h, so that

q(t) = \pi r^2 \frac{dh}{dt} = FK(H-h) (8.31)

where F is a factor that depends on the shape and dimensions of the piezometer intake. If q = q0 at t = 0, it is clear that q(t) will decrease asymptotically toward zero as time goes on.

Figure 8.20 Hvorslev piezometer test. (a) Geometry; (b) method of analysis.

Hvorslev defined the basic time lag, T0, as t = 0, is

T_0 = \frac{\pi r^2}{FK} (8.32)

When this parameter is substituted in Eq. (8.31), the solution to the resulting ordinary differential equation, with the initial condition, h = H0 at t = 0, is

\frac{H-h}{H-H_0} = e^{-t/T_0} (8.33)

A plot of field recovery data, H – h versus t, should therefore show an exponential decline in recovery rate with time. If, as shown on Figure 8.20(b), the recovery is normalized to HH0 and plotted on a logarithmic scale, a straight-line plot results. Note that for Hh/HH0 = 0.37, ln(Hh/HH0) = –1, and from Eq. (8.33), T0 = t. The basic time lag, T0, can be defined by this relation; or if a more physical definition is desired, it am be seen, by multiplying both top and bottom of Eq. (8.32) by HH0, that T0 is the time that would be required for the complete equalization of the head difference if the original rate of inflow were maintained. That is, T0 = V/q0 where V is the volume of water removed or added.

To interpret a set of field recovery data, the data are plotted in the form of Figure 8.20(b). The value of T0 is measured graphically, and K is determined from Eq. (8.32). For a piezometer intake of length L and radius R [Figure 8.20(a)], with L/R > 8, Hvorslev (1951) has evaluated the shape factor, F. The resulting expression for K is

K = \frac{r^2 \text{ln} (L/R)}{2LT_0} (8.34)

Hvorslev also presents formulas for anisotropic conditions and for a wide variety of shape factors that treat such cases as a piezometer open only at its basal cross section and a piezometer that just encounters a permeable formation underlying an impermeable one. Cedergren (1967) also lists these formulas.

In the field or agricultural hydrology, several in situ techniques, similar in principle to the Hvorslev method but differing in detail, have been developed for the measurement of saturated hydraulic conductivity. Boersma (1965) and Bouwer and Jackson (1974) review those methods that involve auger holes and piezometers.

For bail tests of slug tests run in piezometers that are open over the entire thickness of a confined aquifer, Cooper et al. (1967) and Papadopoulos et al. (1973) have evolved a test-interpretation procedure. Their analysis is subject to the same assumptions as the Theis solution for pumpage from a confined aquifer. Contrary to the Hvorslev method of analysis, it includes consideration of both formation and water compressibilities. It utilizes a curve-matching procedure to determine the aquifer coefficients T and S. The hydraulic conductivity K can then be determined on the basis of the relation, K = T/b. Like the Theis solution, the method is based on the solution to a boundary-value problem that involves the transient equation of groundwater flow, Eq. (2.77). The mathematics will not be described here.

For the bail-test geometry shown in Figure 8.21(a), the method involves the preparation of a plot of recovery data in the form Hh/HH0 versus t. The plot is prepared on semilogarithmic paper with the reverse format to that of the Hvorslev test; the Hh/HH0 scale is linear, while the t scale is logarithmic. The field curve is then superimposed on the type curves shown in Figure 8.21(b).

Figure 8.21 Piezometer test in a confined aquifer. (a) Geometry; (b) type curves (after Papadopoulos et al., 1973).

With the axes coincident, the data plot is translated horizontally into a position where the data best fit one of the type curves. A matchpoint is chosen (or rather, a vertical axis is matched) and values of t and W are read of the horizontal scales at the matched axis of the field plot and the type plot, respectively. For ease of calculation it is common to choose a matched axis at W = 1.0. The transmissivity T is then given by

T = \frac{Wr^2}{t} (8.35)

where the parameters are expressed in any consistent set of units.

In principle, the storativity, S, can be determined from the a value of the matched curve and the expression shown on Figure 8.21(b). In practice, since the slopes of the various a lines are very similar, the determination of S by this method is unreliable.

The main limitation on slug tests and hail tests is that they are heavily dependent on a high-quality piezometer intake. If the wellpoint or screen is corroded or clogged, measured values may be highly inaccurate. On the other hand, if a piezometer is developed by surging or backwashing prior to testing, the measured values may reflect the increased conductivities in the artificially induced gravel pack around the intake.

It is also possible to determine hydraulic conductivity in a piezometer or single well by the introduction of a tracer into the well bore. The tracer concentration decreases with time under the influence of the natural hydraulic gradient that exists in the vicinity of the well. This approach is known as the borehole dilution method, and it is described more fully in Section 9.4.

8.6 Measurement of Parameters: Pumping Tests

In this section, a method of parameter measurement that is specifically suited to the determination of transmissivity and storativity in confined and unconfined aquifers will be described. Whereas laboratory tests provide values of the hydrogeological parameters, and piezometer tests provide in situ values representative of a small volume of porous media in the immediate vicinity of a piezometer tip, pumping tests provide in situ measurements that are averaged over a large aquifer volume.

The determination of T and S from a pumping test involves a direct application of the formulas developed in Section 8.3. There, it was shown that for a given pumping rate, if T and S are known, it is possible to calculate the time rate of drawdown, h0h versus t at any point in an aquifer. Since this response depends solely an the values of T and S, it should be possible to take measurements of h0h versus t at some observational point in an aquifer and work backward through the equations to determine the values ofTand S.

The usual course of events during the initial exploitation of an aquifer involves (1) the drilling of a test well with one or more observational piezometers, (2) a short-term pumping test to determine the values of T and S, and (3) application of the predictive formulas of Section 8.3, using the T and S values determined in the pumping test, to design a production well or wells that will fulfill the pumpage requirements of the project without lending to excessive long-term drawdowns. The question of what constitutes an “excessive” drawdown and how drawdowns and well yields are related to groundwater recharge rates and the natural hydrologic cycle are discussed in Section 8.10.

Let us now examine the methodology of pumping-test interpretation in more detail. There are two methods that are in common usage for calculating aquifer coefficients from time-drawdown data. Both approaches are graphical. The first involves curve matching on a log-log plot (the Theis method), and the second involves interpretations with a semilog plot (the Jacob method).

Log-Log Type-Curve Matching

Let us first consider data taken from an aquifer in which the geometry approaches that of the idealized Theis configuration. As was explained in connection with Figure 8.5, the time-drawdown response in an observational piezometer in such an aquifer will always have the shape of the Theis curve, regardless of the values of T and S in the aquifer. However, for high T a measurable drawdown will reach the observation point faster than for low T, and the drawdown data will begin to march up the Theis curve sooner. Theis (1935) suggested the following graphical procedure to exploit this curve-matching property:

  1. Plot the function W(u) versus 1/u on log-log paper. (Such a plot of dimensionless theoretical response is known as a type curve.)
  2. Plot the measured time-drawdown values, h0h versus t, on log-log paper of the same size and scale as the W(u) versus 1/u curve.
  3. Superimpose the field curve on the type curve keeping the coordinate axes parallel. Adjust the curves until most of the observed data points fall on the type curve.
  4. Select and arbitrary match point and read off the paired values of W(u), 1/u, h0h, and t at the match point. Calculate u from 1/u.
  5. Using these values, together with the pumping rate Q and the radial distance r from well to piezometer, calculate T from the relationship

T = \frac{QW(u)}{4\pi (h_0-h)} (8.36)

  1. Calculate S from the relationship

S = \frac{4yTt}{r^2} (8.37)

Equations (8.36) and (8.37) follow directly from Eqs. (8.7) and (8.6). They are valid for any consistent system of units. Some authors prefer to present the equations in the form

T = \frac{AQW(u)}{h_0 - h} (8.38)

S = \frac{uTt}{Br^2} (8.39)

where the coefficients A and B are dependent on the units used for the various parameters. For SI units, with h0h and r measured in meters, t in seconds, Q in m3/s, and T in m2/s, A = 0.08 and B = 0.25. For the inconsistent set of practical units widely used in North America, with h0h and r measured in feet, t in days, Q in US. gal/min, and T in U.S. gal/day/ft, A = 114.6 and B = 1.87. For Q and T in terms of Imperial gallons, A remains unchanged and B = 1.56.

Figure 8.22 illustrates the curve-matching procedure and calculations for a set of field data. The alert reader will recognize these data as being identical to the calculated data originally presented in Figure 8.5(b). It would probably be intuitively clearer if the match point were taken at some point on the coincident portions of the superimposed curves. However, a few quick calculations should convince doubters that it is equally valid to take the matchpoint anywhere on the overlapping fields once they have been fixed in their correct relative positions. For ease of calculation, the matchpoint is often taken at W(u) = 1.0, u = 1.0.

Figure 8.22 Determination of T and S from h0h versus t data using the log-log curve-matching procedure and the W(u) versus 1/u-type curve.

The log-log curve-matching technique can also be used for leaky aquifers (Walton, 1962) and unconfined aquifers (Prickett, 1965; Neuman, 1975a). Figure 8.23 provides a comparative review of the geometry of these systems and the types of h0h versus t data that should be expected in an observational piezometer in each case. Sometimes time-drawdown data unexpectedly display one of these forms, thus indicating a geological configuration that has gone unrecognized during the exploration stage of aquifer evaluation.

Figure 8.23 Comparison of log-log h0h versus t data for ideal, leaky, unconfined, and bounded systems.

For leaky aquifers the time-drawdown data can be matched against the leaky type curves of Figure 8.8. The r/B value of the matched curve, together with the matchpoint values of W(u, r/B), u, h0h and t can be substituted into Eqs. (8.6), (8.8), and (8.9) to yield the aquifer coefficients T and S. Because the development of the r/B solutions does not include consideration of aquitard storativity, an r/B curve matching approach is not suitable for the determination of the aquitard conductivity K’. As noted in the earlier subsection on aquitard response, there are many aquifer-aquitard configurations where the leakage properties of the aquitards are more important in determining long-term aquifer yields than the aquifer parameters themselves. In such cases it is necessary to design a pumping-test configuration with observational piezometers that bottom in the aquitards as well as in the aquifers. One can then use the pumping-test procedure outlined by Neuman and Witherspoon (1972), which utilizes their more general leaky-aquifer solution embodied in Eqs. (8.6), (8.10), and (8.11). They present a ratio method that obviates the necessity of matching field data to type curves as complex as those of Figure 8.9. The method only requires matching against the Theis curve, and calculations are relatively easy to carry out.

As an alternative approach (Wolff, 1970), one can simply read off a Tf value from Figure 8.17 given a hydraulic head value h measured in an aquitard piezometer at elevation z at time t. Knowing the aquitard thickness, b’, one can solve Eq. (8.23) for cv. If an α value can be estimated, Eq. (8.22) can be solved for K’.

For unconfined aquifers the time-drawdown data should be matched against the unconfined type curves of Figure 8.12. The η value of the matched curve, together with the match-point values of W(uA, uB, η), uA, uB, h0h and t can be substituted into Eqs. (8.13) through (8.15) to yield the aquifer coefficients T, S, and Sy. Moench and Prickett (1972) discuss the interpretation of data at sites where lowered water levels cause a conversion from confined to unconfined conditions.

Figure 8.23(d) shows the type of log-log response that would be expected in the vicinity of an impermeable or constant-head boundary. However, bounded systems are more easily analyzed with the semilog approach that will now be described.

Semiology Plots

The semilog method of pump-test interpretation rests on the fact that the exponential integral, W(u), in Eqs. (8.5) and (8.7) can be represented in infinite series. The Theis solution then becomes

h_0 - h = \frac{Q}{4 \pi T} \left( -0.5772 - \text{ln } u + u - \frac{u^2}{2 \cdot 2!} +  \frac{u^3}{3 \cdot 3!} + ...\right) (8.40)

Cooper and Jacob (1946) noted that for small u the sum of the series beyond ln u becomes negligible, so that

h_0 - h = \frac{Q}{4\pi T}(-0.5772 - \text{ln } u) (8.41)

Substituting Eq. (8.6) for u, and noting that ln u = 2.3 log u, that –ln u = ln u, and that ln 1.78 = 0.5772, Eq. (8.41) becomes

h_0 - h = \frac{2.3Q}{4\pi T} \log \frac{2.25Tt}{r^2S} (8.42)

Since Q, r, T, and S are constants, it is clear that h0h versus log t should plot as a straight line.

Figure 8.24(a) shows the time-drawdown data of Figure 8.22 plotted on a semilog graph. If Δh is the drawdown for one log cycle of time and t0 is the time intercept where the drawdown line intercepts the zero drawdown axis, it follows from further manipulation with Eq. (8.42) that the values of T and S, in consistent units, are given by

T = \frac{2.3Q}{4\pi \Delta T} (8.43)

S = \frac{2.25Tt_0}{r^2} (8.44)

As with the log-log methods, these equations can be reshaped as

T = \frac{CQ}{\Delta h} (8.45)

S = \frac{DTt_0}{r^2} (8.46)

where C and D are coefficients that depend on the units used. For Δh and r in meters, t in seconds, Q in m3/s, and T in m2/s, C = 0.18 and D = 2.25. For Δh and r in feet, t in days, Q in U.S. gal/min, and T in U.S. gal/day/ft, C = 264 and D = 0.3. For and T in terms of Imperial gallons, C = 264 and D = 0.36.

Figure 8.24 (a) Determination of T and S from h0h versus t data using the semilog method; (b) semilog plot in the vicinity of an impermeable boundary.

Todd (1959) states that the semilog method is valid for u < 0.01. Examination of the definition of u [Eq. (8.6)] shows that this condition is most likely to be satisfied for piezometers at small r and large t.

The semilog method is very well suited to the analysis of bounded confined aquifers. As we have seen, the influence of a boundary is equivalent to that of a recharging or discharging image well. For the case of an impermeable boundary, for example, the effect of the additional imaginary pumping well is to double the slope of the h0h versus log t plot [Figure 8.24(b)]. The aquifer coefficients S and T should be calculated from Eqs. (8.43) and (8.44) on the earliest limb of the plot (before the influence of the boundary is felt). The time, t1, at which the break in slope takes place can be used together with Eqs. (8.19) to calculate ri, the distance from piezometer to image well [Figure 8.15(c)]. It takes records from three piezometers to unequivocally locate the position of the boundary if it is not known from geological evidence.

Advantages and Disadvantages of Pumping Tests

The determination of aquifer constants through pumping tests has become a standard step in the evaluation of groundwater resource potential. In practice, there is much art to successful pump testing and the interested reader is directed to Kruseman and de Ridder (1970) and Stallman (1971) for detailed advice on the design of pumping-test geometries, and to Walton’s (1970) many case histories.

The advantages of the method are probably self-evident. A pumping test provides in situ parameter values, and these values are, in effect, averaged over a large and representative aquifer volume. One obtains information on both conductivity (through the relation K = T/b) and storage properties from a single test. In aquifer-aquitard systems it is possible to obtain information on the very important leakage properties of the system if observations are made in the aquitards as well as the aquifers.

There are two disadvantages, one scientific and one practical. The scientific limitation relates to the nonuniqueness of pumping-test interpretation. A perusal of Figure 8.23(b), (c), and (d) indicates the similarity in time-drawdown response that can arise from leaky, unconfined, and bounded systems. Unless there is very clear geologic evidence to direct groundwater hydrologists in their interpretation, there will be difficulties in providing a unique prediction of the effects of any proposed pumping scheme. The fact that a theoretical curve can be matched by pumping test data in no way proves that the aquifer fits the assumptions on which the curve is based.

The practical disadvantage of the method lies in its expense. The installation of test wells and observational piezometers to obtain aquifer coefficients is probably only justified in cases where exploitation of the aquifer by wells at the test site is contemplated. In most such cases, the test well can be utilized as a production well in the subsequent pumping program. In geotechnical applications, in contamination studies, in regional flow-net analysis, or in any flow-net approach that requires hydraulic conductivity data but is not involved with well development, the use of the pumping-test approach is usually inappropriate. It is our opinion that the method is widely overused. Piezometer tests are simpler and cheaper, and they can provide adequate data in many cases where pumping tests are not justified.

8.7 Estimation of Saturated Hydraulic Conductivity

It has long been recognized that hydraulic conductivity is related to the grain-size distribution of granular porous media. In the early stages of aquifer exploration or in regional studies where direct permeability data are sparse, this interrelationship can prove useful for the estimation of conductivity values. In this section, we will examine estimation techniques based on grain-size analyses and porosity determinations. These types of data are often widely available in geological reports, agricultural soil surveys, or reports of soil mechanics testing at engineering sites.

The determination of a relation between conductivity and soil texture requires the choice of a representative grain-size diameter. A simple and apparently durable empirical relation, due to Hazen in the latter part of the last century, relies on the effective grain size, d10, and predicts a power-law relation with K:

K = Ad^2_{10} (8.47)

The d10 value can be taken directly from a grain-size gradation curve as determined by sieve analysis. It is the grain-size diameter at which 10% by weight of the soil particles are finer and 90% are coarser. For K in cm/s and d10 in mm, the coefficient A in Eq. (8.47) is equal to 1.0. Hazen’s approximation was originally determined for uniformly graded sands, but it can provide rough but useful estimates for most soils in the fine sand to gravel range.

Textural determination of hydraulic conductivity becomes more powerful when some measure of the spread of the gradation curve is taken into account. When this is done, the median grain size, d50, is usually taken as the representative diameter. Masch and Denny (1966) recommend plotting the gradation curve [Figure 8.25(a)] using Krumbein’s φ units, where φ = –log2d, d being the grain-size diameter (in mm). As a measure of spread, they use the inclusive standard deviation, σ1, where

\sigma_1 = \frac{d_{16}-d_{84}}{4} + \frac{d_5 - d_{95}}{6.6} (8.48)

For the example shown in Figure 8.25(a), d50 = 2.0 and σ1 = 0.8. The curves shown in Figure 8.25(b) were developed experimentally in the laboratory on prepared samples of unconsolidated sand. From them, one can determine K, knowing d50, and σ1.

Figure 8.25 Determination of saturated hydraulic conductivity from grain-size gradation curves for unconsolidated sands (after Masch and Denny, 1966).

For a fluid of density, ρ, and viscosity, μ, we have seen in Section 2.3 [Eq. (2.26)] that the hydraulic conductivity of a porous medium consisting of uniform spherical grains of diameter, d, is given by

K = \left(\frac{\rho g}{\mu}\right)Cd^2 (8.49)

For a nonuniform soil, we might expect the d in Eq. (8.49) to become dm, where dm is some representative grain size, and we would expect the coefficient C to be dependent on the shape and packing of the soil grains. The fact that the porosity, n, represents an integrated measure of the packing arrangement has led many investigators to carry out experimental studies of the relationship between C and n. The best known of the resulting predictive equations for hydraulic conductivity is the Kozeny-Carmen equation (Bear, 1972), which takes the form

K = \left(\frac{\rho g}{\mu}\right)\left[ \frac{n^3}{(1-n)^2} \right]\left( \frac{d^2_m}{180} \right) (8.50)

In most formulas of this type, the porosity term is identical to the central element of Eq. (8.50), but the grain-size term can take many forms. For example, the Fair-Hatch equation, as reported by Todd (1959), take the form

K = \left(\frac{\rho g}{\mu}\right)\left[ \frac{n^3}{(1-n)^2} \right]\left[ \frac{1}{m \left( \frac{\theta}{100}\sum \frac{P}{d_m} \right)^2} \right] (8.51)

where m is a packing factor, found experimentally to be about 5; θ is a sand shape factor, varying from 6.0 for spherical grains to 7.7 for angular grains; P is the percentage of sand held between adjacent sieves; and dm is the geometric mean of the rated sizes of adjacent sieves.

Both Eqs. (8.50) and (8.51) are dimensionally correct. They are suitable for application with any consistent set of units.

8.8 Prediction of Aquifer Yield by Numerical Simulation

The analytical methods that were presented in Section 8.3 for the prediction of drawdown in multiple-well systems are not sophisticated enough to handle the heterogeneous aquifers of irregular shape that are often encountered in the field. The analysis and prediction of aquifer performance in such situations is normally carried out by numerical simulation on a digital computer.

There are two basic approaches: those that involve a finite-difference formulation, and those that involve a finite-element formulation. We will look at finite-difference methods in moderate detail, but our treatment of finite-element methods will be very brief.

Finite-Difference Methods

As with the steady-state finite-difference methods that were described in Section 5.3, transient simulation requires a discretization of the continuum that makes up the region of flow. Consider a two-dimensional, horizontal, confined aquifer of constant thickness, b; and let it be discretized into a finite number of blocks, each with its own hydrogeologic properties, and each having a node at the center at which the hydraulic head is defined for the entire block. As shown in Figure 8.26(a), some of these blocks may be the site of pumping wells that are removing water from the aquifer.

Let us now examine the flow regime in one of the interior nodal blocks and its four surrounding neighbors. The equation of continuity for transient, saturated flow states that the net rate of flow into any nodal block must be equal to the time rate of change of storage within the nodal block. With reference to Figure 8.26(b), and following the developments of Section 2.11, we have

Q_{15} + Q_{25} + Q_{35} + Q_{45} = S_{S_5} \Delta x \Delta y b \frac{\partial h_5}{\partial t} (8.52)

where SS5 is the specific storage of nodal block 5. From Darcy’s law,

Q_{15} = K_{15} \frac{h_1 - h_5}{\Delta y}\Delta x \hspace{1mm} b (8.53)

where K15 is a representative hydraulic conductivity between nodes 1 and 5. Similar expressions can be written for Q25, Q35, and Q45.

Figure 8.26 Discretization of two-dimensional, horizontal, confined aquifer.

Let us first consider the case of a homogeneous, isotropic medium for which K15 = K25 = K35 = K45 = K and SS1 = SS2 = SS3 = SS4 = SS. If we arbitrarily select a square nodal grid with Δx = Δy, and note that T = Kb and S = SSb, substitution of expression such as that of Eq. (8.53) into Eq. (8.52) leads to

T(h_1 + h_2 + h_3 + h_4 - 4h_5) = S \Delta x^2 \frac{\partial h_5}{\partial t} (8.54)

The time derivative on the right-hand side can be approximated by

\frac{\partial h_5}{\partial t} = \frac{h_5(t) - h_5(t - \Delta t)}{\Delta t} (8.55)

where Δt is the time step that is used to discretize the numerical model in a time-wise sense. If we now convert to the ijk notation indicated on Figure 8.26(c), where the subscript (i, j) refers to the nodal position and the superscript k = 0, 1, 2, . . . indicates the time step, we have

h^k_{i, j-1} + h^k_{i+1,j} + h^k_{i-1, j} + h^k_{i, j+1} - 4h^k_{i,j} = \frac{S \Delta x^2}{T \Delta t}(h^k_{i,j}-h^{k-1}_{i,j}) (8.56)

In a more general form,

Ah^k_{i,j} = B^k_{i,j-1}+Ch^k_{i+1,j}+Dh^k_{i,j+1}+Eh^k_{i-1,j} + F (8.57)


A = \frac{S \Delta x^2}{T \Delta t} + 4 (8.58)

B = C = D = E = 1 (8.59)

F = \frac{S \Delta x^2}{T \Delta t} \cdot h^{k-1}_{i,j} (8.60)

Equation (8.57) is the finite-difference equation for an internal node (i, j) in a homogeneous, isotropic, confined aquifer. Each of the parameters S, T, Δx, and Δt that appear in the definitions of the coefficients are known, as is the value of the hydraulic head, h(i, j), at the previous time step, k – 1. In a similar fashion, it is possible to develop finite-difference equations for boundary nodes and corner nodes, and for nodes from which pumping takes place. In each case, the finite-difference equation is similar in form to Eq. (8.57), but the expressions for the coefficients will differ. For boundary nodes, some of the coefficients will be zero. For an internal pumping node, the coefficients A, B, C, D, and E are as given in Eqs. (8.58) and (8.59), but

F = \frac{\Delta x^2}{T}\left( \frac{S}{\Delta t} \cdot h^{k-1}_{i,j} + W_{i,j} \right) (8.61)

where Wi,j is a sink term with units [L/T]. W is related to the pumping rate, Q[L3/T] by

W_{i,j} = \frac{Q_{i,j}}{\Delta x^2} (8.62)

Sometimes W is given a more general definition,

W_{i,j} = \frac{Q_{i,j}}{\Delta x^2} - R_{i,j} (8.63)

where Ri,j is a source term with units [L/T] that represents vertical leakage into the aquifer from overlying aquitards. In this case Eq. (8.61) is used for all nodes in the system and Wi,j is specified for every node. It will be negative for nodes accepting leakage and positive for nodes undergoing pumping.

It is possible to develop Eq. (8.57) in a more rigorous way, starting with the partial differential equation that describes transient flow in a horizontal confined aquifer. In Appendix IX, the rigorous approach is used to determine the values for the coefficients A, B, C, D, E and F, in the general finite-difference equation for an internal node in a heterogeneous, anisotropic aquifer. In such a system each node (i, j) may be assigned its own specific values of Si, j, (Tx)i, j, (Ty)i, j, where Tx and Ty are the principal components of the transmissivity tensor in the x and y coordinate directions. The derivation of Appendix IX is carried out for a rectangular nodal grid in which Δx ≠ Δy. A further sophistication, which is not considered there, would allow an irregular nodal grid in which the Δx and Δy values are themselves a function of nodal position. Irregular nodal spacings are often required in the vicinity of pumping wells where hydraulic gradients tend to be large. The concepts that underlie the development of these more complex finite-difference formulations is identical to that which led to Eq. (8.57). The more complex the finite-difference equations embodied in the computer program, the more versatile is that program as a numerical simulator of aquifer performance.

It is possible, then, to develop a finite-difference equation, at some degree of sophistication, for every node in the nodal grid. If there are N nodes, there are N finite-difference equations. At each time step, there are also N unknowns: namely, the N values of h(i, j) at the N nodes. At each time step, we have N linear, algebraic equations in N unknowns. This set of equations must be solved simultaneously at each time step, starting from a set of initial conditions wherein h(i, j) is known for all (i, j), and proceeding through the time steps, k = 1, 2, . . . . Many methods are available for the solution of the system of equations, and numerical aquifer models are often classified on the basis of the approach that is used. For example, the method of successive over relaxation that was described in Section 5.3 for the numerical simulation of steady-state flow nets is equally applicable to the system of equations that arises at each time step of a transient aquifer model. More commonly a method known as the alternating-direction implicit procedure is used. Remson et al. (1971) and Pinder and Gray (1977) provide a systematic and detailed presentation of these various methods as they pertain to aquifer simulation. Advanced mathematical treatment of the methods is available in the textbook by Forsythe and Wasow (1960). The original development of most numerical-simulation techniques took place in the petroleum engineering field, where the primary application is in the simulation of oil-reservoir behavior. Pinder and Bredehoeft (1968) adapted the powerful alternating-direction implicit procedure to the needs of groundwater hydrologists.

There are two aquifer-simulation programs that have been completely documented and widely applied in North America. One is the U.S. Geological Survey model, which is an outgrowth of Pinder and Bredehoeft’s original work. Trescott et al. (1976) provide an updated manual for the most recent version of the computer program. The other is the Illinois State Water Survey model, which is fully documented by Prickett and Lonnquist (1971). Bredehoeft and Pinder (1970) have also shown how a sequence of two-dimensional aquifer models can be coupled together to form a quasi-three-dimensional model of an aquifer-aquitard system.

As a practical example, we will consider the analysis carried out by Pinder and Bredehoeft (1968) for an aquifer at Musquoduboit Harbour, Nova Scotia. The aquifer there is a glaciofluvial deposit of limited areal extent. Figure 8.27(a) shows the initial estimate of the areal distribution of transmissivity for the aquifer as determined from the rather sparse hydrogeological data that were available. Simulations with this transmissivity matrix failed to reproduce the drawdown patterns observed during a pumping test that was run near the center of the aquifer. The aquifer parameters were then adjusted and readjusted over several computer runs until a reasonable duplication was achieved between the measured time-drawdown data and the results of the digital model. Additional test-well logs tended to support the adjusted parameters at several points. The final transmissivity distribution is shown in Figure 8.27(b). The model was then put into prediction mode; Figure 8.27(c) is a plot of the predicted drawdown pattern 206.65 days after the start of exploitation by a proposed production well pumping at a rate of Q = 0.963 ft3/s.

Render (1971, 1972) and Huntoon (1974) provide additional case histories of interest.

Finite-Element Methods

The finite-element method, first noted in Section 5.3 in connection with the simulation of steady-state flow nets, can also be used for the simulation of transient aquifer performance. As in the finite-difference approach, the finite-element approach leads to a set of N algebraic equations in N unknowns at each time step, where the N unknowns are the values of the hydraulic heads at a set of nodal points distributed through the aquifer. The fundamental difference lies in the nature of the nodal grid. The finite-element method allows the design of an irregular mesh that can be hand-tailored to any specific application. The number of nodes can often be significantly reduced from the number required for a finite-difference simulation. The finite-element approach also has some advantages in the way it treats boundary conditions and in the simulation of anisotropic media.

The development of the finite-element equations for each node requires an understanding of both partial differential equations and the calculus of variations. Remson, Hornberger, and Molz (1971) provide an introductory treatment of the method as it applies to aquifer simulation. Pinder and Gray (1977) provide an advanced treatment. Zienkiewicz (1967) and Desai and Abel (1972) are the most widely quoted general reference texts. The finite-element method was introduced into the groundwater literature by Javandel and Witherspoon (1969). Pinder and Frind (1972) were among the first to utilize the method for the prediction of regional aquifer performance. Gupta and Tanji (1976) have reported an application of a three-dimensional finite-element model for the simulation of flow in an aquifer-aquitard system in the Sutter Basin, California.

Model Calibration and the Inverse Problem

If measurements of aquifer transmissivity and storativity were available at every nodal position in an aquifer-simulation model, the prediction of drawdown patterns would be a very straightforward matter. In practice, the data base on which models must be designed is often very sparse, and it is almost always necessary to calibrate the model against historical records of pumping rates and drawdown patterns. The parameter adjustment procedure that was described in connection with Figure 8.27 represents the calibration phase of the modeling procedure for that particular example. In general, a model should be calibrated against one period of the historical record, then verified against another period of record. The application of a simulation model for a particular aquifer then becomes a three-step process of calibration, verification, and prediction.

Figure 8.27 Numerical simulation of aquifer performance at Musquoduboit Harbour, Nova Scotia (after Pinder and Bredehoeft, 1968).

Figure 8.28 is a flowchart that clarifies the steps involved in the repetitive trial-and-error approach to calibration. Parameter correction may be carried out on the basis of purely empirical criteria or with a performance analyzer that embodies formal optimization procedures. The contribution by Neuman (1973a) includes a good review and a lengthy reference list. The role of subjective information in establishing the constraints for optimization was treated by Lovell et al. (1972). Gates and Kisiel (1974) considered the question of the worth of additional data. They analyzed the trade-off between the cost of additional measurements and the value they have in improving the calibration of the model.

Figure 8.28 Flowchart of the trial-and-error calibration process (after Neuman, 1973a).

The term calibration usually refers to the trial-and-error adjustment of aquifer parameters as outlined in Figure 8.28. This approach involves the repetitive application of the aquifer model in its usual mode. In each simulation the boundary-value problem is set up in the usual way with the transmissivity, T(x, y), storativity, S(x, y), leakage, R(x, y, t), and pumpage, Q(x, y, t), known, and the hydraulic head, h(x, y, t), unknown. It is possible to carry out the calibration process more directly by utilizing an aquifer simulation model in the inverse mode. In this case only a single application of the model is required, but the model must be set up as an inverse boundary-value problem where h(x, y, t) and Q(x, y, t) are known and T(x, y), S(x, y) and R(x, y, t) are unknown. When posed in this fashion, the calibration process is known as the inverse problem.

In much of the literature, the term parameter identification is used to encompass all facets of the problem at hand. What we have called calibration is often called the indirect approach to the parameter identification problem, and what we have called the inverse problem is called the direct approach.

The solution of the inverse formulation is not, in general, unique. In the first place there may be too many unknowns; and in the second place, h(x, y, t) and Q(x, y, t) are not known for all (x, y). In practice, pumpage takes place at a finite number of points, and the historical records of head are available at only a finite number of points. Even if R(x, y, t) is assumed constant or known, the problem remains ill-posed mathematically. Emsellem and de Marsily (1971) have shown, however, that the problem can be made tractable by using a “flatness criteria” that limits the allowable spatial variations in T and S. The mathematics of their approach is not simple, but their paper remains the classic discussion of the inverse problem. Newman (1973a, 1975b) suggests using available measurements of T and S to impose constraints on the structure of T(x, y) and S(x, y) distributions. The contributions of Yeh (1975) and Sagar (1975) include reviews of more recent developments.

There is another approach to inverse simulation that is simpler in concept but apparently open to question as to its validity (Neuman, 1975b). It is based on the assumption of steady-state conditions in the flow system. As first recognized by Stallman (1956), the steady-state hydraulic head pattern, h(x, y, z) in a three-dimensional system can be interpreted inversely in terms of the hydraulic conductivity distribution, K(x, y, z). In a two-dimensional, unpumped aquifer, h(x, y) can be used to determine T(x, y). Nelson (1968) showed that the necessary condition for the existence and uniqueness of a solution to the steady-state inverse problem is that, in addition to the hydraulic heads, the hydraulic conductivity or transmissivity must be known along a surface crossed by all streamlines in the system. Frind and Pinder (1973) have pointed out that, since transmissivity and flux are related by Darcy’s law, this criterion can be stated alternatively in terms of the flux that crosses a surface. If water is being removed from an aquifer at a steady pumping rate, the surface to which Nelson refers occurs around the circumference of the well and the well discharge alone provides a sufficient boundary condition for a unique solution. Frind and Pinder (1973) utilized a finite-element model to solve the steady-state inverse problem. Research is continuing on the question of what errors are introduced into the inverse solution when a steady-state approach is used for model calibration for an aquifer that has undergone a transient historical development.

8.9 Prediction of Aquifer Yield by Analog Simulation

Numerical simulation of aquifer performance requires a moderately large computer and relatively sophisticated programming expertise. Electric-analog simulation provides an alternative approach that circumvents these requirements at the expense of a certain degree of versatility.

Analogy Between Electrical Flow and Groundwater Flow

The principles underlying the physical and mathematical analogy between electrical flow and groundwater flow were introduced in Section 5.2. The application under discussion was the simulation of steady-state flow nets in two-dimensional vertical cross sections. One of the methods described there utilized a resistance-network analog that was capable of handling heterogeneous systems of irregular shape, In this section, we will pursue analog methods further, by considering the application of two-dimensional resistance-capacitance networks for the prediction of transient hydraulic-head declines in heterogeneous, confined aquifers of irregular shape.

Consider a horizontal confined aquifer of thickness b. If it is overlaid with a square grid of spacing, ΔxA [as in Figure 8.26(a)], any small homogeneous portion of the discretized aquifer [Figure 8.29(a)] can be modeled by a scaled-down array of electrical resistors and capacitors on a square grid of spacing, ΔxM [Figure 8.29(b)]. The analogy between electrical flow in the resistance-capacitance network and groundwater flow in the horizontal confined aquifer can be revealed by examining the finite-difference form of the equations of flow for each system. For ground-water flow, from Eq. (8.54),

T(h_1 + h_2 + h_3 + h_4 - 4h_5) = S\Delta x^2_A \frac{\partial h_5}{\partial t_A} (8.64)

Figure 8.29 Small homogeneous portion of discretized aquifer and analogous resistor-capacitor network (after Prickett, 1975).

For the electrical circuit, from Kirchhoff’s laws:

\frac{1}{R}(V_1 + V_2 + V_3 + V_4 - 4V_5) = C\frac{\partial V_5}{\partial t_M} (8.65)

Comparison of Eqs. (8.64) and (8.65) leads to the analogous quantities:

  1. Hydraulic head, h; and voltage, V.
  2. Transmissivity, T; and the reciprocal of the resistance, R, of the resistors.
  3. The product of the storativity, S, times the nodal block area, Δx2A; and the capacitance, C, of the capacitors.
  4. Aquifer coordinates, xA and yA (as determined by the spacing. ΔxA); and model coordinates, xM and yM (as determined by the spacing, ΔxM).
  5. Real time, tA; and model time, tM.

In addition, if pumpage is considered, there is an analogy between:

  1. Pumping rate, Q, at a well; and current strength, I, at an electrical source.

Resistance-Capacitance Network

The network of resistors and capacitors that constitutes the analog model is usually mounted on a Masonite pegboard perforated with holes on approximately 1-inch centers. There are four resistors and one capacitor connected to each terminal. The resistor network is often mounted on the front of the board, and the capacitor network, with each capacitor connected to a common ground, on the back. The boundary of the network is designed in a stepwise fashion to approximate the shape of the actual boundary of the aquifer.

The design of the components of the analog requires the choice of a set of scale factors, F1, F2, F3, and F4, such that

F_1 = \frac{h}{V} (8.66)

F_2 = \frac{\Delta x_A}{\Delta x_M} (8.67)

F_3 = \frac{t_A}{t_M} (8.68)

F_4 = \frac{Q}{I} (8.69)

Heterogeneous and transversely anisotropic aquifers can be simulated by choosing resistors and capacitors that match the transmissivity and storativity at each point in the aquifer. Comparison of the hydraulic flow through an aquifer section and the electrical flow through an analogous resistor [Figure 8.30(a)] leads to the relation

R = \frac{F_4}{F_1T} (8.70)

Comparison of the storage in an aquifer section and the electrical capacitance of an analogous capacitor [Figure 8.30(b)] leads to the relation

C = \frac{F_1S\Delta x^2_A}{F_4F_3} (8.71)

Figure 8.30 Aquifer nodal block and analogous (a) resistor and (b) capacitor (after Prickett, 1975).

The resistors and capacitors that make up the network are chosen on the basis of Eqs. (8.70) and (8.71). The scale factors, F1, F2, F3, and F4, must be selected in such a way that (1) the resistors and capacitors fall within the range of inexpensive, commercially available components; (2) the size of the model is practical; and (3) the response times of the model are within the range of available excitation-response equipment.

Figure 8.31 is a schematic diagram that shows the arrangement of excitation-response apparatus necessary for electric-analog simulation using a resistance-capacitance network. The pulse generator, in tandem with a waveform generator, produces a rectangular pulse of specific duration and amplitude. This input pulse is displayed on channel 1 of a dual-channel oscilloscope as it is fed through a resistance box to the specific terminal of the resistance-capacitance network that represents the pumped well. The second channel on the oscilloscope is used to display the time-voltage response obtained by probing various observation points in the network. The input pulse is analogous to a step-function increase in pumping rate; the time-voltage graph is analogous to a time-drawdown record at an observational piezometer. The numerical value of the head drawdown is calculated from the voltage drawdown by Eq. (8.66).

Figure 8.31 Excitation-response apparatus for electrical-analog simulation using a resistance-capacitance network.

The time at which any specific drawdown applies is given by Eq. (8.68). Any pumping rate, Q, may be simulated by setting the current strength, I, in Eq. (8.69). This is done by controlling the resistance, R_i, of the resistance box in Figure 8.31. The current strength is given by I = Vi/Ri, where Vi is the voltage drop across the resistance box.

Walton (1970) and Prickett (1975) provide detailed coverage of the electric-analog approach to aquifer simulation. Most groundwater treatments owe much to the general discussion of analog simulation by Karplus (1958). Results of analog simulation are usually presented in the form of maps of predicted water-level drawdowns similar to that shown in Figure 8.27(c). Patten (1965), Moore and Wood (1967), Spieker (1968), and Render (1971) provide case histories that document the application of analog simulation to specific aquifers.

Comparison of Analog and Digital Simulation

Prickett and Lonnquist (1968) have discussed the advantages, disadvantages, and similarities between analog and digital techniques of aquifer simulation. They note that both methods use the same basic field data, and the same method of assigning hydrogeologic properties to a discretized representation of the aquifer. Analog simulation requires knowledge of specialized electronic equipment; digital simulation requires expertise in computer programming. Digital simulation is more flexible in its ability to handle irregular boundaries and pumping schemes that vary through time and space. It is also better suited to efficient data readout and display.

The physical construction involved in the preparation of a resistance-capacitance network is both the strength and the weakness of the analog method. The fact that the variables of the system under study are represented by analogous physical quantities and pieces of equipment is extremely valuable for the purposes of teaching or display, but the cost in time is large. The network, once built, describes only one specific aquifer. In digital modeling, on the other hand, once a general computer program has been prepared, data decks representing a wide variety of aquifers and aquifer conditions can be run with the same program. The effort involved in designing and keypunching a new data deck is much less than that involved in designing and building a new resistance-capacitance network. This flexibility is equally important during the calibration phase of aquifer simulation.

The advantages of digital simulation weigh heavily in its favor, and with the advent of easy accessibility to large computers, the method is rapidly becoming the standard tool for aquifer management. However, analog simulation will undoubtedly continue to play a role for some time, especially in developing countries where computer capacities are not yet large.

8.10 Basin Yield

Safe Yield and Optimal Yield of a Groundwater Basin

Groundwater yield is best viewed in the context of the full three-dimensional hydrogeologic system that constitutes a groundwater basin. On this scale of study we can turn to the well-established concept of safe yield or to the more rigorous concept of optimal yield.

Todd (1959) defines the safe yield of a groundwater basin as the amount of water that can be withdrawn from it annually without producing an undesired result. Any withdrawal in excess of safe yield is an overdraft. Domenico (1972) and Kazmann (1972) review the evolution of the term. Domenico notes that the “undesired results” mentioned in the definition are now recognized to include not only the depletion of the groundwater reserves, but also the intrusion of water of undesirable quality, the contravention of existing water rights, and the deterioration of the economic advantages of pumping. One might also include excessive depletion of streamflow by induced infiltration and land subsidence.

Although the concept of safe yield has been widely used in groundwater resource evaluation, there has always been widespread dissatisfaction with it (Thomas, 1951; Kazmann, 1956). Most suggestions for improvement have encouraged consideration of the yield concept in a socioeconomic sense within the overall framework of optimization theory. Domenico (1972) reviews the development of this approach, citing the contributions of Bear and Levin (1967), Buras (1966), Burt (1967), Domenico et al. (1968), and others. From an optimization viewpoint, groundwater has value only by virtue of its use, and the optimal yield must be determined by the selection of the optimal groundwater management scheme from a set of possible alternative schemes. The optimal scheme is the one that best meets a set of economic and/or social objectives associated with the uses to which the water is to be put. In some cases and at some points in time, consideration of the present and future costs and benefits may lead to optimal yields that involve mining groundwater, perhaps even to depletion. In other situations, optimal yields may reflect the need for complete conservation. Most often, the optimal groundwater development lies somewhere between these extremes.

The graphical and mathematical methods of optimization, as they relate to groundwater development, are reviewed by Domenico (1972).

Transient Hydrologic Budgets and Basin Yield

In Section 6.2 we examined the role of the average annual groundwater recharge, R, as a component in the steady-state hydrologic budget for a watershed. The value of R was determined from a quantitative interpretation of the steady-state, regional, groundwater flow net. Some authors have suggested that the safe yield of a groundwater basin be defined as the annual extraction of water that does not exceed the average annual groundwater recharge. This concept is not correct. As pointed out by Bredehoeft and Young (1970), major groundwater development may significantly change the recharge-discharge regime as a function of time. Clearly, the basin yield depends both on the manner in which the effects of withdrawal are transmitted through the aquifers and on the changes in rates of groundwater recharge and discharge induced by the withdrawals. In the form of a transient hydrologic budget for the saturated portion of a groundwater basin,

Q(t) = R(t) - D(t) + \frac{dS}{dt} (8.72)


Q(t) = total rate of groundwater withdrawal

R(t) = total rate of groundwater recharge to the basin

D(t) = total rate of groundwater discharge from the basin

dS/dt = rate of change of storage in the saturated zone of the basin.

Freeze (1971a) examined the response of R(t) and D(t) to an increase in Q(t) in a hypothetical basin in a humid climate where water tables are near the surface. The response was simulated with the aid of a three-dimensional transient analysis of a complete saturated-unsaturated system such as that of Figure 6.10 with a pumping well added. Figure 8.32 is a schematic representation of his findings. The diagrams show the time-dependent changes that might be expected in the various terms of Eq. (8.72) under increased pumpage. Let us first look at the case shown in Figure 8.32(a), in which withdrawals increase with time but do not become excessive. The initial condition at time t0 is a steady-state flow system in which the recharge, R0, equals the discharge, D0. At times t1, t2, t3, and t4, new wells begin to tap the system and the pumping rate Q undergoes a set of stepped increases. Each increase is initially balanced by a change in storage, which in an unconfined aquifer takes the form of an immediate water-table decline. At the same time, the basin strives to set up a new equilibrium under conditions of increased recharge, R.

Figure 8.32 Schematic diagram of transient relationships between recharge rates, discharge rates, and withdrawal rates (after Freeze, 1971a).

The unsaturated zone will now be induced to deliver greater flow rates to the water table under the influence of higher gradients in the saturated zone. Concurrently, the increased pumpage may lead to decreased discharge rates, D. In Figure 8.32(a), after time t4, all natural discharge ceases and the discharge curve rises above the horizontal axis, implying the presence of induced recharge from a stream that had previously been accepting its baseflow component from the groundwater system. At time t5, the withdrawal Q is being fed by the recharge, R, and the induced recharge, D; and there has been a significant decline in the water table. Note that the recharge rate attains a maximum between t3 and t4. At this rate, the groundwater body is accepting all the infiltration that is available from the unsaturated zone under the lowered water-table conditions.

In Figure 8.32(a), steady-state equilibrium conditions are reached prior to each new increase in withdrawal rate. Figure 8.32(b) shows the same sequence of events under conditions of continuously increasing groundwater development over several years. This diagram also shows that if pumping rates are allowed to increase indefinitely, an unstable situation may arise where the declining water table reaches a depth below which the maximum rate of groundwater recharge R can no longer be sustained. After this point in time the same annual precipitation rate no longer provides the same percentage of infiltration to the water table. Evapotranspiration during soil-moisture-redistribution periods now takes more of the infiltrated rainfall before it has a chance to percolate down to the groundwater zone. At t4 in Figure 8.32(b), the water table reaches a depth below which no stable recharge rate can be maintained. At t5 the maximum available rate of induced recharge is attained. From time t5 on, it is impossible for the basin to supply increased rates of withdrawal. The only source lies in an increased rate of change of storage that manifests itself in rapidly declining water tables. Pumping rates can no longer be maintained at their original levels. Freeze (1971a) defines the value of Q at which instability occurs as the maximum stable basin yield. To develop a basin to its limit of stability would, of course, be foolhardy. One dry year might cause an irrecoverable water-table drop. Production rates must allow for a factor of safety and must therefore be somewhat less than the maximum stable basin yield.

The discussion above emphasizes once again the important interrelationships between groundwater flow and surface runoff. If a groundwater basin were developed up to its maximum yield, the potential yields of surface-water components of the hydrologic cycle in the basin would be reduced. It is now widely recognized that optimal development of the water resources of a watershed depend on the conjunctive use of surface water and groundwater. The subject has provided a fertile field for the application of optimization techniques (Maddock, 1974; Yu and Haimes, 1974). Young and Bredehoeft (1972) describe the application of digital computer simulations of the type described in Section 8.8 to the solution of management problems involving conjunctive groundwater and surface-water systems.

8.11 Artificial Recharge and Induced Infiltration

In recent years, particularly in the more populated areas of North America where water resource development has approached or exceeded available yield, there has been considerable effort placed on the management of water resource systems. Optimal development usually involves the conjunctive use of groundwater and surface water and the reclamation and reuse of some portion of the available water resources. In many cases, it involves the importation of surface water from areas of plenty to areas of scarcity, or the conservation of surface water in times of plenty for use in times of scarcity. These two approaches require storage facilities, and there is often advantage to storing water underground where evaporation losses are minimized. Underground storage may also serve to replenish groundwater resources in areas of overdraft.

Any process by which man fosters the transfer of surface water into the groundwater system can be classified as artificial recharge. The most common method involves infiltration from spreading basins into high-permeability, unconfined, alluvial aquifers. In many cases, the spreading basins are formed by the construction of dikes in natural channels. The recharge process involves the growth of a groundwater mound beneath the spreading basin. The areal extent of the mound and its rate of growth depend on the size and shape of the recharging basin, the duration and rate of recharge, the stratigraphic configuration of subsurface formations, and the saturated and unsaturated hydraulic properties of the geologic materials. Figure 8.33 shows two simple hydrogeological environments and the type of groundwater mound that would be produced in each case beneath a circular spreading basin. In Figure 8.33(a), recharge takes place into a horizontal unconfined aquifer bounded at the base by an impermeable formation. In Figure 8.33(b), recharge takes place through a less-permeable formation toward a high-permeability layer at depth.

Figure 8.33 Growth of a groundwater mound beneath a circular recharge basin.

Both cases have been the subject of a large number of predictive analyses, not only for circular spreading basins but also for rectangular basins and for recharge from an infinitely long strip. The latter case, with boundary conditions like those shown in Figure 8.33(b), also has application to canal and river seepage. It has been studied in this context by Bouwer (1965), Jeppson (1968), and Jeppson and Nelson (1970). The case shown in Figure 8.33(a), which also has application to the development of mounds beneath waste disposal ponds and sanitary landfills, has been studied in even greater detail. Hantush (1967) provides an analytical solution for the prediction of h(r, t), given the initial water-table height, h0, the diameter of the spreading basin, a, the recharge rate, R, and the hydraulic conductivity and specific yield, K and Sy, of the unconfined aquifer. His solution is limited to homogeneous, isotropic aquifers and a recharge rate that is constant in time and space. In addition, the solution is limited to a water-table rise that is less than or equal to 50% of the initial depth of saturation, h0. This requirement implies that R \ll K. Bouwer (1962) utilized an electric-analog model to analyze the same problem, and Marine (1975a, 1975b) produced a numerical simulation. All three of these analyses have two additional limitations. First, they neglect unsaturated flow by assuming that the recharge pulse traverses the unsaturated zone vertically and reaches the water table unaffected by soil moisture-conditions above the water table. Second, they utilize the Dupuit-Forchheimer theory of unconfined flow (Section 5.5) which neglects any vertical flow gradients that develop in the saturated zone in the vicinity of the mound. Numerical simulations carried out on the complete saturated-unsaturated system using the approaches of Rubin (1968), Jeppson and Nelson (1970), and Freeze (1971a) would provide a more accurate approach to the problem, but at the expense of added complexity in the calculations.

Practical research on spreading basins has shown that the niceties of predictive analysis are seldom reflected in the real world. Even if water levels in spreading ponds are kept relatively constant, the recharge rate almost invariably declines with time as a result of the buildup of silt and clay on the basin floor and the growth of microbial organisms that clog the soil pores. In addition, air entrapment between the wetting front and the water table retards recharge rates. Todd (1959) notes that alternating wet and dry periods generally furnish a greater total recharge than does continuous spreading. Drying kills the microbial growths, and tilling and scraping of the basin floor during dry periods reopens the soil pores.

There are several excellent case histories that provide an account of specific projects involving artificial recharge from spreading basins. Seaburn (1970) describes hydrologic studies carried out at two of the more than 2000 recharge basins that are used on Long Island, east of New York City, to provide artificial recharge of storm runoff from residential and industrial areas. Bianchi and Haskell (1966, 1968) describe the piezometric monitoring of a complete recharge cycle of mound growth and dissipation. They report relatively good agreement between the field data and analytical predictions based on Dupuit-Forchheimer theory. They note, however, that the anomalous water-level rises that accompany air entrapment (Section 6.8) often make it difficult to accurately monitor the growth of the groundwater mound.

While water spreading is the most ubiquitous form of artificial recharge, it is limited to locations with favorable geologic conditions at the surface. There have also been some attempts made to recharge deeper formations by means of injection wells. Todd (1959) provides several case histories involving such diverse applications as the disposal of storm-runoff water, the recirculation of air-conditioning water, and the buildup of a freshwater barrier to prevent further intrusion of seawater into a confined aquifer. Most of the more recent research on deep-well injection has centered on utilization of the method for the disposal of industrial wastewater and tertiary-treated municipal wastewater (Chapter 9) rather than for the replenishment of groundwater resources.

The oldest and most widely used method of conjunctive use of surface water and groundwater is based on the concept of induced infiltration. If a well produces water from alluvial sands and gravels that are in hydraulic connection with a stream, the stream will act as a constant-head line source in the manner noted in Figures 8.15(d) and 8.23(d). When a new well starts to pump in such a situation the pumped water is initially derived from the groundwater zone, but once the cone of depression reaches the stream, the source of some of the pumped water will be streamflow that is induced into the groundwater body under the influence of the gradients set up by the well. In due course, steady-state conditions will be reached, after which time the cone of depression and the drawdowns within it remain constant. Under the steady flow system that develops at such times, the source of all the pumped groundwater is streamflow. One of the primary advantages of induced infiltration schemes over direct surface-water utilization lies in the chemical and biological purification afforded by the passage of stream water through the alluvial deposits.

8.12 Land Subsidence

In recent years it has become apparent that the extensive exploitation of groundwater resources in this century has brought with it an undesired environmental side effect. At many localities in the world, groundwater pumpage from unconsolidated aquifer-aquitard systems has been accompanied by significant land subsidence. Poland and Davis (1969) and Poland (1972) provide descriptive summaries of all the well-documented cases of major land subsidence caused by the withdrawal of fluids. They present several case histories where subsidence has been associated with oil and gas production, together with a large number of cases that involve groundwater pumpage. There are three cases—the Wilmington oil field in Long Beach, California, and the groundwater overdrafts in Mexico City, Mexico, and in the San Joaquin valley, California—that have led to rates of subsidence of the land surface of almost 1 m every 3 years over the 35-year period 1935–1970. In the San Joaquin valley, where groundwater pumpage for irrigation purposes is to blame, there are three separate areas with significant subsidence problems. Taken together, there is a total area of 11,000 km2 that has subsided more than 0.3 m. At Long Beach, where the subsiding region is adjacent to the ocean, subsidence has resulted in repeated flooding of the harbor area. Failure of surface structures, buckling of pipe lines, and rupturing of oil-well casing have been reported. Remedial costs up to 1962 exceeded $100 million.

Mechanism of Land Subsidence

The depositional environments at the various subsidence sites are varied, but there is one feature that is common to all the groundwater-induced sites. In each case there is a thick sequence of unconsolidated or poorly consolidated sediments forming an interbedded aquifer-aquitard system. Pumpage is from sand and gravel aquifers, but a large percentage of the section consists of high-compressibility clays. In earlier chapters we learned that groundwater pumpage is accompanied by vertical leakage from the adjacent aquitards. It should come as no surprise to find that the process of aquitard drainage leads to compaction*Following Poland and Davis (1969), we are using the term “compaction” in its geological sense. In engineering jargon the term is often reserved for the increase in soil density achieved through the use of rollers, vibrators or other heavy machinery. of the aquitards just as the process of aquifer drainage leads to compaction of the aquifers. There are two fundamental differences, however: (1) since the compressibility of clay is 1–2 orders of magnitude greater than the compressibility of sand, the total potential compaction of an aquitard is much greater than that for an aquifer; and (2) since the hydraulic conductivity of clay may be several orders of magnitude less than the hydraulic conductivity of sand, the drainage process, and hence the compaction process, is much slower in aquitards than in aquifers.

Consider the vertical cross section shown in Figure 8.34. A well pumping at a rate Q is fed by two aquifers separated by an aquitard of thickness b.

Figure 8.34 One-dimensional consolidation of an aquitard.

Let us assume that the geometry is radially symmetric and that the transmissivities in the two aquifers are identical. The time-dependent reductions in hydraulic head in the aquifers (which could be predicted from leaky-aquifer theory) will be identical at points A and B. We wish to look at the hydraulic-head reductions in the aquitard along the line AB under the influence of the head reductions in the aquifers at A and B. If hA(t) and h_B(t) are approximated by step functions with a step Δh (Figure 8.34), the aquitard drainage process can be viewed as the one-dimensional, transient boundary-value problem described in Section 8.3 and presented as Eq. (8.21). The initial condition is h = h0 all along AB, and the boundary conditions are h = h0 – Δh at A and at B for all t > 0. A solution to this boundary-value problem was obtained by Terzaghi (1925) in the form of an analytical expression for h(z, t). An accurate graphical presentation of his solution appears as Figure 8.17. The central diagram on the right-hand side of Figure 8.34 is a schematic plot of his solution; it shows the time-dependent decline in hydraulic head at times t0, t1 . . . , t ∞ along the line AB. To obtain quantitative results for a particular case, one must know the thickness b’, the vertical hydraulic conductivity K’, the vertical compressibility α’, and the porosity n’ of the aquitard, together with the head reduction Δh on the boundaries.

In soil mechanics the compaction process associated with the drainage of a clay layer is known as consolidation. Geotechnical engineers have long recognized that for most clays α \gg , so the latter term is usually omitted from Eq. (8.21). The remaining parameters are often grouped into a single parameter cv, defined by

C_v = \frac{K'}{\rho g \alpha '} (8.73)

The hydraulic head h(z, t) can be calculated from Figure 8.17 with the aid of Eq. (8.23) given cv, Δh, and b.

In order to calculate the compaction of the aquitard given the hydraulic head declines at each point on AB as a function of time, it is necessary to recall the effective stress law: σT = σe + p. For σT = constant, e = -dp. In the aquitard, the head reduction at any point z between the times t1 and t2 (Figure 8.34) is dh = h1(z, t1) – h2(z, t2). This head drop creates a fluid pressure reduction: dp = ρg dψ = ρg d(hz) = pg dh, and the fluid pressure reduction is reflected by an increase in the effective stress e = -dp. It is the change in effective stress, acting through the aquitard compressibility α’, that causes the aquitard compaction Δb’. To calculate Δb’ along AB between the times t1 and t2, it is necessary to divide the aquitard into m slices. Then, from Eq. (2.54),

\Delta b'_{t_1 - t_2} = b' \sum^m_{i=1} \rho g \alpha ' dh_i (8.74)

where dhi is the average head decline in the ith slice.

For a multiaquifer system with several pumping wells, the land subsidence as a function of time is the summation of all the aquitard and aquifer compactions. A complete treatment of consolidation theory appears in most soil mechanics texts (Terzaghi and Peck, 1967; Scott, 1963). Domenico and Mifflin (1965) were the first to apply these solutions to cases of land subsidence.

It is reasonable to ask whether land subsidence can be arrested by injecting groundwater back into the system. In principle this should increase the hydraulic heads in the aquifers, drive water back into the aquitards, and cause an expansion of both aquifer and aquitard. In practice, this approach is not particularly effective because aquitard compressibilities in expansion have only about one-tenth the value they have in compression. The most successful documented injection scheme is the one undertaken at the Wilmington oil field in Long Beach, California (Poland and Davis, 1969). Repressuring of the oil reservoir was initiated in 1958 and by 1963 there had been a modest rebound in a portion of the subsiding region and the rates of subsidence were reduced elsewhere.

Field Measurement of Land Subsidence

If there are any doubts about the aquitard-compaction theory of land subsidence, they should be laid to rest by an examination of the results of the U.S. Geological Survey subsidence research group during the last decade. They have carried out field studies in several subsiding areas in California, and their measurements pro- vide indisputable confirmation of the interrelationships between hydraulic head declines, aquitard compaction, and land subsidence.

Figure 8.35 is a contoured map, based on geodetic measurements, of the land subsidence in the Santa Clara valley during the period 1934–1960.

Figure 8.35 Land subsidence in feet, 1934–1960, Santa Clara valley, California (after Poland and Davis, 1969).

Subsidence is confined to the area underlain by unconsolidated deposits of alluvial and shallow-marine origin. The centers of subsidence coincide with the centers of major pumping, and the historical development of the subsidence coincides with the period of settlement in the valley and with the increased utilization of groundwater.

Quantitative confirmation of the theory is provided by results of the type shown in Figure 8.36. An ingeniously simple compaction-recorder installation [Figure 8.36(a)] produces a graph of the time-dependent growth of the total compaction of all material between the land surface and the bottom of the hole.

Figure 8.36 (a) Compaction-recorder installation; (b) compaction measurement site near Sunnyvale, California; (c) measured compactions, land subsidence, and hydraulic head variations at the Sunnyvale site, 1960–1962 (after Poland and Davis, 1969).

Near Sunnyvale in the Santa Clara valley, three compaction recorders were established at different depths in the confined aquifer system that exists there [Figure 8.36(b)]. Figure 8.36(c) shows the compaction records together with the total land subsidence as measured at a nearby benchmark, and the hydraulic head for the 250- to 300-m-depth range as measured in an observation well at the measurement site. Decreasing hydraulic heads are accompanied by compaction. Increasing hydraulic heads are accompanied by reductions in the rate of compaction, but there is no evidence of rebound. At this site “the land subsidence is demonstrated to be equal to the compaction of the water-bearing deposits within the depth tapped by water wells, and the decline in artesian head is proved to be the sole cause of the subsidence” (Poland and Davis, 1969, p. 259).

Riley (1969) noted that data of the type shown on Figure 8.36(c) can be viewed as the result of a large-scale field consolidation test. If the reductions in aquitard volume reflected by the land subsidence are plotted against the changes in effective stress created by the hydraulic-head declines, it is often possible to calculate the average compressibility and the average vertical hydraulic conductivity of the aquitards. Helm (1975, 1976) has carried these concepts forward in his numerical models of land subsidence in California.

It is also possible to develop predictive simulation models that can relate possible pumping patterns in an aquifer-aquitard system to the subsidence rates that will result. Gambolati and Freeze (1973) designed a two-step mathematical model for this purpose. In the first step (the hydrologic model), the regional hydraulic-head drawdowns are calculated in an idealized two-dimensional vertical cross section in radial coordinates, using a model that is a boundary-value problem based on the equation of transient groundwater flow. Solutions are obtained with a numerical finite-element technique. In the second step of the modeling procedure (the subsidence model), the hydraulic head declines determined with the hydrologic model for the various aquifers are used as time-dependent boundary conditions in a set of one-dimensional vertical consolidation models applied to a more refined geologic representation of each aquitard. Gambolati et al. (1974a, 1974b) applied the model to subsidence predictions for Venice, Italy. Recent measurements summarized by Carbognin et al. (1976) verify the model’s validity.

8.13 Seawater Intrusion

When groundwater is pumped from aquifers that are in hydraulic connection with the sea, the gradients that are set up may induce a flow of salt water from the sea toward the well. This migration of salt water into freshwater aquifers under the influence of groundwater development is known as seawater intrusion.

As a first step toward understanding the nature of the processes involved, it is necessary to examine the nature of the saltwater-freshwater interface in coastal aquifers under natural conditions. The earliest analyses were carried out independently by two European scientists (Ghyben, 1888; Herzberg, 1901) around the turn of the century. Their analysis assumed simple hydrostatic conditions in a homogeneous, unconfined coastal aquifer. They showed [Figure 8.37(a)] that the interface separating salt water of density ρs and fresh water of density ρf must project into the aquifer at an angle α < 90°.

Figure 8.37 Saltwater-freshwater interface in an unconfined coastal aquifer (a) under hydrostatic conditions; (b) under conditions of steady-state seaward flow (after Hubbert, 1940).

Under hydrostatic conditions, the weight of a unit column of fresh water extending from the water table to the interface is balanced by a unit column of salt water extending from sea level to the same depth as the point on the interface. With reference to Figure 8.37(a), we have

\rho_s gz_s = \rho_f g(z_s + z_w) (8.75)


z_s = \frac{\rho_f}{\rho_s - \rho_f}z_w (8.76)

For ρf = 1.0 and ρs = 1.025,

z_s = 40z_w (8.77)

Equation (8.77) is often called the Ghyben-Herzberg relation.

If we specify a change in the water-table elevation of Δzw, then from Eq. (8.77), Δzs = 40 Δzw. If the water table in an unconfined coastal aquifer is lowered 1 m, the saltwater interface will rise 40 m.

In most real situations, the Ghyben-Herzberg relation underestimates the depth to the saltwater interface. Where freshwater flow to the sea takes place, the hydrostatic assumptions of the Ghyben-Herzberg analysis are not satisfied. A more realistic picture was provided by Hubbert (1940) in the form of Figure 8.37(b) for steady-state outflow to the sea. The exact position of the interface can be determined for any given water-table configuration by graphical flow-net construction, noting the relationships shown on Figure 8.37(b) for the intersection of equipotential lines on the water table and on the interface.

The concepts outlined in Figure 8.37 do not reflect reality in yet another way. Both the hydrostatic analysis and the steady-state analysis assume that the interface separating fresh water and salt water in a coastal aquifer is a sharp boundary. In reality, there tends to be a mixing of salt water and fresh water in a zone of diffusion around the interface. The size of the zone is controlled by the dispersive characteristics of the geologic strata. Where this zone is narrow, the methods of solution for a sharp interface may provide a satisfactory prediction of the fresh-water flow pattern, but an extensive zone of diffusion can alter the flow pattern and the position of the interface, and must be taken into account. Henry (1960) was the first to present a mathematical solution for the steady-state case that includes consideration of dispersion. Cooper et al. (1964) provide a summary of the various analytical solutions.

Seawater intrusion can be induced in both unconfined and confined aquifers. Figure 8.38(a) provides a schematic representation of the saltwater wedge that would exist in a confined aquifer under conditions of natural steady-state outflow. Initiation of pumping [Figure 8.38(b)] sets up a transient flow pattern that leads to declines in the potentiometric surface on the confined aquifer and inland migration of the saltwater interface. Pinder and Cooper (1970) presented a numerical mathematical method for the calculation of the transient position of the saltwater front in a confined aquifer. Their solution includes consideration of dispersion.

Figure 8.38 (a) Saltwater-freshwater interface in a confined coastal aquifer under conditions of steady-state seaward flow; (b) seawater intrusion due to pumping.

One of the most intensively studied coastal aquifers in North America is the Biscayne aquifer of southeastern Florida (Kohout, 1960a, 1960b). It is an unconfined aquifer of limestone and calcareous sandstone extending to an average depth of 30m below sea level. Field data indicate that the saltwater front undergoes transient changes in position under the influence of seasonal recharge patterns and the resulting water-table fluctuations. Lee and Cheng (1974) and Segol and Pinder (1976) have simulated transient conditions in the Biscayne aquifer with finite-element numerical models. Both the field evidence and the numerical modeling confirm the necessity of considering dispersion in the steady-state and transient analyses. The nature of dispersion in groundwater flow will be considered more fully in Chapter 9 in the context of groundwater contamination.

Todd (1959) summarizes five methods that have been considered for controlling seawater intrusion: (1) reduction or rearrangement of the pattern of groundwater pumping, (2) artificial recharge of the intruded aquifer from spreading basins or recharge wells, (3) development of a pumping trough adjacent to the coast by means of a line of pumping wells parallel to the coastline, (4) development of a freshwater ridge adjacent to the coast by means of a line of recharge wells parallel to the coastline, and (5) construction of an artificial subsurface barrier. Of these five alternatives, only the first has been proven effective and economic. Both Todd (1959) and Kazmann (1972) describe the application of the freshwater-ridge concept in the Silverado aquifer, an unconsolidated, confined, sand-and-gravel aquifer in the Los Angeles coastal basin of California. Kazmann concludes that the project was technically successful, but he notes that the economics of the project remain a subject of debate.

*Following Poland and Davis (1969), we are using the term “compaction” in its geological sense. In engineering jargon the term is often reserved for the increase in soil density achieved through the use of rollers, vibrators or other heavy machinery.

Suggested Readings

BOUWER, H., and R. D. JACKSON. 1974. Determining soil properties. Drainage for Agriculture, ed. J. van Schilfgaarde. American Society of Agronomy, Madison, Wis., pp. 611–672.

COOPER, H. H. JR., F. A. KOHOUT, H. R. HENRY, and R. E. GLOVER. 1964. Sea water in coastal aquifers. U.S. Geol. Surv. Water-Supply Paper 1613C, 84 pp.

FERRIS, J. G., D. B. KNOWLES, R. H. BROWNE, and R. W. STALLMAN. 1962. Theory of aquifer tests. U.S. Geol. Surv. Water-Supply Paper I536E.

HANTUSH, M. S. 1964. Hydraulics of wells. Adv. Hydrosci., 1, pp. 281–432.

KRUSEMAN, G. P., and N. A. DE RIDDER. 1970. Analysis and evaluation of pumping test data. Intern. Inst. for Land Reclamation and Improvement Bull. 11, Wageningen, The Netherlands.

NEUMAN, S. P., and P. A. WITHERSPOON. 1969. Applicability of current theories of flow in leaky aquifers. Water Resources Res., 5, pp. 817–829.

POLAND, J. F., and G. H. DAVIS. 1969. Land subsidence due to withdrawal of fluids. Geol. Soc. Amer. Rev. Eng. Geol., 2, pp. 187–269.

PRICKETT, T. A. 1975. Modeling techniques for groundwater evaluation. Adv. Hydrosci., 11, pp. 46–66, 91–116.

REMSON, I., G. M. HORNBERGER, and F. J. MOLZ. 1971. Numerical Methods in Subsurface Hydrology. Wiley Interscience, New York, pp. 56–122.

STALLMAN, R. W. 1971. Aquifer-test design, observation and data analysis. Techniques of Water Resources Investigations of the U.S. Geological Survey, Chapter B1. Government Printing Office, Washington, D.C.

YOUNG, R. A., and J. D. BREDEHOEFT. 1972. Digital computer simulation for solving management problems of conjunctive groundwater and surface-water systems. Water Resources Res., 8, pp. 533–556.