Государственный комитет Республики Узбекистан

B = C = D = E = 1(8.59)
(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
(8.61)
where W_i,j is a sink term with units [L/T]. W is related to the pumping rate, Q[L³/T] by
(8.62)
Sometimes W is given a more general definition,
(8.63)
where R_i,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 W_i,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 S_{i, j}, (T_x)_{i, j}, and (T_y)_{i, j}, where T_x and T_y 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 ft³/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, Δx_A [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, Δx_M [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),
(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:
(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, Δx²_A; and the capacitance, C, of the capacitors.
   
4. 
   Aquifer coordinates, x_A and y_A (as determined by the spacing. Δx_A); and model coordinates, x_M and y_M (as determined by the spacing, Δx_M).
   
5. 
   Real time, t_A; and model time, t_M.
   

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

6. 
   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, F₁, F₂, F₃, and F₄, such that
(8.66)
(8.67)
(8.68)
(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
(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
(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, F₁, F₂, F₃, and F₄, 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). 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 = V_i/R_i, where V_i is the voltage drop across the resistance box.
Figure 8.31 Excitation-response apparatus for electrical-analog simulation using a resistance-capacitance network.
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,
(8.72)
where