- Email: [email protected]
Preliminary results from a methodology study in the field of groundwater resources management
152
Mathematics and Computers in Simulation XXIV (1982) 152-160 North-Holland Publishing Company
PRELIMINARY RESULTS FROM A METHODOLOGY GROUNDWATER RESOURCES MANAGEMENT
M. BENEDINI, G. GIULIANO,
B. SIRANGELO
STUDY IN THE FIELD OF
and S. TROISI
Water Research Institute, National Research Council, Rome, Italy
Increasing water demand, conflicts among utilizations, water-body pollution have given a considerable impulse to the application of systems analysis techniques also to the field of groundwater resources management. In line with this, some mathematical simulation and optimization models taking into account the various aspects of management have been prepared. The first of these models, known as IDROSIM, is actually a package including various flow models developed both for the problem of parameter identification (distribution of transmissivity and porosity parameters, etc.) and for the so called direct problem. The second model, called IASO, deals with modifications in water quality as a function of pollution from urban, industrial or agricultural sources. The third model, called IDROPROG, sets out to apply linear programming techniques to the aquifer in order to optimize water resource utilization. Finally, problems of combined management of both groundwater and surface water are investigated by means of a fourth model, dubbed IDRINTEGRO. All these models are currently being tested on a coastal aquifer extending over an area of 600 km* on the Ionian side of the Puglia region of Italy. This paper makes a comprehensive survey of the principles on which the models are based, and provides further details on the IDROSIM simulation model.
1. INTRODUCTION There is a growing trend in groundwater studies to seek more efficient ways of utilizing aquifers. Many different aspects need to be considered in selecting operative criteria for planning the use of ground water resources. "Nonphysical" aspects such as the economic, social and regulatory implications may also be taken into account. The systems approach to water resources management has been applied to the field of surface water for a long time now 1 and world-wide research now under way tends to suggest that similar methods could be profitably extended to the field of groundwater. However, before going into the application of such techniques to groundwater,some peculiarities of aquifer systems need to be discussed. The physical properties of the water bearing material, i.e. its permeability, porosity, storage, etc., differ at virtually every point of the aquifer. Each of these properties affects the aquifer's hydraulic behaviour with relevant points in many management problems such as
0378-4754/82/0000-0000/$02.75
those connected with pollutant dispersion and propagation of impulses in the piezometric surface. As a result, these characteristics are extremely important in decisions as to how to use the groundwater resources. The same tools for describing aquifer behavior both quantitatively and qualitatively are required also in the analysis of the above mentioned non-physical contingencies. For example, by representing how an aquifer is influenced by pumping it is possible to evaluate the amount of water that can be drawn from an aquifer; as well as its quality.This in turn makes it possible to draw up the cost estimates needed to evaluate the economics of using that particular resource. 2. NEW METHODS FOR GROUNDWATER RESOURCES MANAGEMENT With the objective of developing a new methodology suitable for groundwater resources management in the most general context, four mathematical models have been worked out.
0 1982 IMACS/North-Holland
M. Benedini et al. /Methodology
study of groundwater
2.1 The first of these,dubbed IDROSIM,is actually a package incorporating models previously developed for solving both "parameter identification problem" and the change of piezometric surface over a certain period of time, under different initial and boundary conditions. Classical studies show that flow field in a non homogeneous, isotropic aquifer can be expressed mathematically by a second order partial differential equation, which, for a confined aquifer bearing in mind Darcy's Law and matrix compressibility features, can be written as:
where:K S h q
= = = =
permeability (LT-1) specific storage coefficient (L-l) piezometric head (L) unit flow rate drawn off from the acquifer or fed in at time t (LT-1)
A similar expression also holds true for an unconfined aquifer, providing the storage coefficient S can be replaced by the effective porosity, E., and variations of head are small with respect to the saturated thickness of the aquifer. 2.2 The second model named IASO deals with the dispersion of a pollutant in an aquifer (in a two-dimensional case) according to the following fundamental expression, the assumption being made that the pollutant does not react chemically with the solid matrix:
resources management
153
Starting from an assessment of the various phenomenological aspects /3/, computer programs are now being developed for the integration of equation (2) on the basis of a groundwater velocity map produced by the IDROSIM model. 2.3 The third model i.e. IDROPROG, embraces mathematical prograsiming procedures capable of determining that piezometric surface pattern which may be considered as "optimal" in relation to some pre-determined targets of groundwater withdrawal. As is known, formulating a mathematical programming model entails maximizing a function of the problem variables, i.e., the so-called "objective function", respecting suitable boundary conditions for the same variables. These conditions are expressed by other functions called "constraints". Of all the programming methods available, linear programming appears to be one of the most useful. In this method both the objective function and the constraints are expressed as linear functions of the variables. In Darcy conditions /4/, considering the aquifer as being discretized into a certain number of grid elements, a head value "h" can be attributed to the i-th grid element. Any sort of limitations on the use of groundwater can then be reduced to the constraints applied to the head, hi2
h. 1
with h. determined "exogenously" depending on the patticular conditions set for the i-th grid element. Sometimes the limitation may also be imposed on the gradient, in which case the constraint can be expressed as
where:c = pollutant concentration D = dispersion coefficient V x = seepage velocity components V i Y This model takes into account the way a pollutant disperses into an aquifer with reference to the determination of the dispersion coefficient and the groundwater velocity.
h
i
-hk"
Iik
for every pair of grid elements, i and k, with a given Iik. Under stationary flow conditions, hydraulic head can be expressed in the form of n hi = hi0 -
2
i diaij
1 where h. is the head corresponding to the initial con'8. itions while Q is the constant flow i
M. Benedini et al. /Methodology
154
study of groundwater
rate at well i. The "a..II term is an "influence coefficient, denoting ia t e effect caused on the i-th grid element by the water draw-off at the j-th grid element. Establishing
resources management
The corresponding costs of Q and QI, C and C can be assumed to be linear Punctions 0!?Q. I This means that the cost for District i can be expressed as:
that
pi
= cpi + CIi
b. = h. -l-Y 1 10 io with a predetermined h. , b. expresses the piezometric level drop ini?he :-th grid element due to aquifer tapping. Thence the final form of the linear constraints is given by
Qi aij 5 r
bi
i
Considering flow rate Q. as "production" and b. as being the upper limit, the above expressioni shows how withdrawal rates should be determined. The objective function may be expressed as: max
Consequently, considering all the districts in the area, the following total cost function will be:
The "manager's" problem will then be to determine the "best" combination between groundwater extracted inside the area and of the exogenous water so as to minimize overall cost function 0 taking into account the boundaries of Q and QI. P Once the minimization will lead to:
District 2
1 -yl=
has been completed, this
Q&Qpmaxl+ qlQI,
'i
subject to the abovementioned
constraints.
2.4 The fourth model, i.e. IDRINTEGRO, deals with problems of groundwater management in the framework of a multi-objective programing scheme. In order to clarify the matter let us suppose that an area, which is under a regional authority's jurisdiction for water resources management, may be divided up into n water districts. For each district monthly demand requirements may be defined that are inclusive of agricultural, urban, industrial and recreational uses. Separately, by means of simulation and programming models the "manager" can decide on the maximum amount of groundwater to be withdrawn in every district. He can also determine the maximum withdrawal rate (Q max) for the whole acquifer falling under his gurisdiction. The deficit, i.e. the difference between the amount of groundwater available monthly and the demand for the different uses, will have to be compensated month-by-month by bringing in water from external sources (QI), e.g. long distance supply, desalination, waste water treatment and so forth.
CX Q P2 pmax2+
District 2 -Cp,=
=
. . . District n -y
n
=
... o( Q n2 pmaxn +
72912
... %lQIn
It is thus possible to indicate how much water must be pumped from the aquifer month-by-month and district-by-district and how much must be imported from outside districts with adequate groundwater supplies. Furthermore, this system enables the "manager" to fix taxation levels for the volume of water required by the users of every district. Our attention will now be focussed on the mathematical model which is more relevant to the subjects of this Working Conference, i.e. the IDROSIM model.
3. THE IDROSIM MODEL This model is designed to simulate piezometric head variations over a certain period of time for a coastal aquifer subject to exploitation and seawater intrusion. IDROSIM is actually a cluster of models which
M. Benedini et al. /Methodology
study of groundwater resources management
perform different calculations. The entire model development may be broken up into two phases, (Fig. 1). Phase 1 - Preparation of a tape containing bulky input data and - Processing of the basic data fed in (data includes aquifer pattern, boundary conditions, transmissivity, storage coefficient and undisturbed piezometric levels). Phase 2 - Access to the tape by other programs controlled by card inputs with final printout. Phase 1 is further sub-divided units:
into several
a) input of aquifer pattern data via Program DATACFG,
r-‘--'--
Fig. 1 - IDROSIM-Program
System Structure.
155
b) input of Transmissivity (T) distribution via the ITMC Identification Program or via Direct Assignment Program ASSIGN T, c) input of Storage Coefficient (S) distribution via the IMMAG Identification Program or via Direct Assignment Program ASSIGN S, d) calculation of ,an Optimal Relaxation Factor (characteristic of every aquifer) allowing the number of iterations in calculating a piezometric surface for pre-set conditions to be minimized. This is achieved using the FINDWE Program, e) calculation of piezometric surface for undisturbed conditions via Program SOLVEHO. Phase 2, on the other hand, envisages the following units: a) elaboration of piezometric head patterns of the aquifer under well pumping rates either for a steady rate of seawater intrusion,via
-. -.-_
-.
(s= card input; c= print output; n= tape output and input).
156
M. Benedini et al. /Methodology
study of groundwater
Program PIEZO, or a varying rate of seawater intrusion, via Program PINDER, b) calculation of Aquifer Influence Coefficient via Program CD 1, of tidal effects on the aquifer C) calculation using either a unidimensional approximation method, via Program DH 1, or a two-dimensional approximation method, via Program DH 2. Of the first phase programs, ITMC is designed to determine transmissivity values over a rectangular area on the assumption that piezometric head data are available for some points within that area. The ITMC Program originates from Equation (1) and assumes that the right hand side is negligible, thus giving:
(3)
resources management
The IMMAG Program sets out to solve the problem of identifying Storage Coefficient S (or effective porosity) in a non-homogeneous aquifer starting from Equation (l), Transmissivity T (ITMCor ASSIGN T - derived) and the in-aquifer tide propagation being known. Identifying S therefore consists of two distinct phases,namely working out piezometric head oscillations on the basis of field measurements, establishing the minimum of the error function between calculated and measured patterns. Here, too, a part of the program extrapolates the results obtained for a part of the aquifer to the adjacent areas. ASSIGN S program uses the same extrapolation method as the ASSIGN T program wherever storage coefficient data are obtained from field data. Before describing the other programs, Equation (1) for steady flow conditions may conveniently be recalled:
From the above equation the function of T can be obtained by the "Method of Characteristics", whenever the function h(x,y) is known and is continuous with its first and second partial derivatives for all points of the aquifer. However, the piezometric head is usually known only in very few points so it becomes necessary to interpolate discrete data in order to obtain an analytical form of h. The "Cubic Natural Spline" was chosen as interpolating function as it has the above mentioned analytical properties (continuity of function and of first and second derivatives) and at the same time gives a good . approximation of both the function h and its derivatives. Consequently, identification of transmissivity values consists of two quite distinct aspects, namely working out the interpolating Spline func tion and determining transmissivity values along characteristic lines. In other words, the problem consists in (a) determining how Parameter T is distributed over a certain area, (b) picking out some physical plausibility constraints for a part of the aquifer /7/, (c) extrapolating these data for the remaining part.
(1’)
The solution approach, based on finite differences, transforms the problem of solving Equation (1') into a problem of solving a system of linear algebraic equations of the general form: A.x = b where: A= x= b=
coefficients matrix unknown vector of heads known vector related to boundary conditions
The i-th equation of System (5) can be expressed by writing a flow rate continuity equation for the i-th grid element of the discretized aquifer. Referring to Fig. 2, the following
__i_____’ i __
w
‘______~__ I N __
c
E
__
__
s These extrapolation methods are also used for the ASSIGN T Program wherever there is a distribution of direct transmissivity data for the several parts of the aquifer.
(5)
__i._.._ @
I _____~_
Fig. '2- Discretizing
Grid
M. Benedini et al. /Methodology
study of groundwater
resources management
1.57
The PINDER Program may be considered as an offholds:
T~~(h~-h~) + TCW(hC-hW)+ TCE(hC-hE)+TCS(hC-hS)= =
(6)
qC
where; T
CN
T CW T CE T
CS
= (TN + TC)/2 = (TW + TC)/Z = (TE + TC)/2
shoot of the PIEZO Program as it may deal with fresh and seawater flows. Most authors have tackled the problem from an analytical point of view, while PINDER and COOPER /lo/ propose a different approach which, although somewhat elaborate, seems sufficiently all-embracing. The problem can be reduced to a system of three equations, the first of which describes fresh water flow, the second concerns saline water flow and the third expresses the relationship between the interface and the fresh and saline water piezometric heads. After some manipulations, the following final equations can be obtained:
= (TS + T&/Z A - For saline water
Reordering, we have,
=2q
C
Moving along the aquifer grid row-by-row, n equations may Be obtained where n indicates the total number of grid elements. Whenever one of the grid elements falls on the boundary, the corresponding term will be expressed as a known quantity.
d%(X-Xi
;
y-Yj
\
)
(8)
/
B - For fresh water
To have a sufficiently detailed numerical solution, discretization in real cases calls for a large system of linear equations (5), often in the order of thousands. Classical direct methodology (elimination, triangular decomposition, etc.) was discarded due to features of coefficient matrix, the possibility of being limited to approximate solutions, and some operative considerations to implement the algorithm on a computer, Iterative calculation methods (Jacobi, Gauss-Seidel, S.O.R., etc.) were preferred. The IDROSEM model utilizes the S.O.R. iterative method and so, in conjunction with the FINDW Program, it can determine the value of the Optimal Relaxation Factor, w , for the aquifer being studied. The Program System can then go on to solve Equation (5) using the SOLVEHO Program where no pumping occurs and the PIEZO Program where it does. The latter program also enables seawater intrusion to be calculated if the sea water conditions can be assumed to be stationary.
(9) where
I/“’ = I
htd’cl
,T_
b’“’ = ),“I+
zi
b’d’
As usual, A C the water water,E the dary of the
is the sea le~el,B the ground level., table,D the piezometric of saline interface and F an impervious bounsaline water field fl,
M. Benedini et al. /Methodology
158
study of groundwater resources management
and 6'= R = b = h = K = S = p =
transmission in a simplified unidimensional approach /8/ or in a more realistic two-dimensio&approach /7/.
withdrawal flow rate per unit area recharging flow rate per unit area Dirac function piezometric head permeability storage coefficient density
Finally, Program CD 1 enables calculating the "influence coefficients" that establish to what extent one pumping well is influenced by another one operating inside the same aquifer. As mentioned in the preceding paragraph, this program acts as a link between the IDROSIM and IDROPROG models with the aim of defining the "influence coefficients" which establish the constraints of the linear programming model.
(9)
6
= porosity
Subscripts x, y and z represent components in their respective directions. Subscripts s and d indicate saline and fresh water, respectively (Fig. 3). The PINDER Program discretizes Equations (8) and (9), and, once integrated by means of the finite difference method, it provides a picture of the piezometric surfaces of both fresh and saline water as well as of the interface in a coastal aquifer at different utilization rates flf. If, on the other hand, an overall picture is required of how a coastal aquifer is affected at its borders, for example, to what extent groundwater and sea are connected, then Programs DH 1 and DH 2 are used to solve the problem of tidal
4. APPLYING NEW METHODS TO GROUNDWATER RESOURCES MANAGEMENT The purpose of this paper is restricted to the presentation of methodologies and is not concerned with any numerical results of their applications. The latter will be dealt with in separate papers in preparation. The described methods are applied o a coastal 5 aquifer extending over some 600 km on the Ionian side (Fig. 4 of Italy's Salento peninsula which was chosen as a case study. The aquifer in question floats on a body of continent-invading seawater. Aquifer thickness tends progressively to zero in proximity of the
I -:_: -
_.-__--____---_ __ ~-_---__-_._- - __ --
_-_---
:.,_
. . . .
Fig. 3 - Fresh water/saline
.
..
water interface in a coastal aquifer.
-
-
-_:
__
M. Benedini et al. /Methodology
GU
study of groundwater
L 0
159
resources management
F
F
TARANTO
i:
Fig. 4 - Location of the investigated aquifer.
coast line where the permeability of the outcropping rocks allows a free seaward outflow of the groundwater. Studies already carried out /9/ on the dynamics of the diffusion zone between the aquifer and underlying seawater provide the necessary data for establishing the aquifer pattern /ll/. The area chosen is bounded upstream by the water divide separating the flow of groundwater towards the Ionian or the Adriatic Seas. This divide is not, however, permanently fixed but is subject to lateral shifts of up to about one kilometer due to the difference in level between the two seas, differences that at some times last even for long periods, To the northwest, the boundary ?s tdentified by the appearance of another bordering hydrogeolo-
gical system (the Murgia system), i.e. by the sudden variation of certain significant features such as piezometric gradients and heads, water flow rate, etc. To the southeast however, a conventional boundary has been adopted taking into account the flow of groundwater entering the zone from the bordering area. In a joint effort with Bari University, Applied Geology and Geotechnics Institute, several datataking campaigns have been carried out on the aquifer to obtain the necessary inputs for the four models dealt with in this paper. The datataking consisted basically in obtaining: - piezometric heads in the different selected water points - both undisturbed and subject to pumping or tide effects,
160
M. Benedini et al. /Methodology
study of groundwater
- local concentration of inorganic chemicals and specific water quality indicators, by means of on site sampling and specific laboratory analysis. - concentration of bacteria, e,g, total coliforms, fecal coliforms and fecal streptococchi. The aquifer was divided up into a grid of 16,000 elements, each measuring 200m x 2OOm, on which pumping and discharge points were also located. ACKNOWLEDGEMENTS The authors particularly wish to thank the technical assistants B. Soria and A. De Fazio for their contribution in the computer work, as well as Mr. R, Jacovelli and the people at SVIM S.p. A. for their fruitful help in preparing the computer programs.
REFERENCES /l/ QUADERNO IRSA 30 - Un esperrmento metodologico di PSanificazione e Gestione delle risorse idriche in un grande bacino idrografico, Roma, 1976, 232 pages. (A Methodology Experiment in the Planning and Managing of the Water Resources of a Large Hydrographic Basic) /2/ CASTELLANO, L. - TROISI, S,: Alcune considerazioni sulla dispersione di un inquinante in una falda sotterranea. Rapport0 tecni‘co IRSA R/79, 1980. /3/ BENEDINI, M. & TROISI, S. - Procedimenti nuo vi per lo studio della gestione delle risorse idriche sotterranee. Geol. Appl. e Idrogeol. Vol. XII, Parte 1, pgs. 79-109, 1977. (New Procedures for the Study of Groundwater Resources Management). 141 SCHWARZ, J. -I Linear Models for Groundwater Management. Jour. of Hydro. 28 (1976), pgs. 377-397. /5/ YU, W & HAINES, Y.Y. - Multi-level Optimization for Conjunctive Use of Groundwater, Water Resources Research, Vol. 10, 11, 4 August, 1974, pgs. 625-635. /6/ CAROTENUTO, L.; DI PILLO, G.; RAICONI, C. & TROISI, S. - Modelli matematici e identificazione parametrica di falde acquifere. Quaderno IRSA 41 pgs.: 96, 1979. (Mathematical Models and the Identification of Aquifer Parameters). /7/ SIRANGELO, B. & TROISI, S. - Studio metodologico di nuove tecniche di gestione delle risorse idriche sotterranee - Applicazione ad
resources management
una falda reale, Quaderno IRSA 50, 1980. (Methodology Study for New Groundwater Management Techniques - A Case Study). /8/ MAGRI, G. & TROISI, S. - Sulla influenza delle fluttuazioni di specchi d'acqua sui livelli delle falde costiere - Applicazione allo studio della circolazione idrica sotterranea nella penisola salentina, Geol. APP~. e Idrogeol., vol. IV, pgs. 25-40,1966. (On Row Free Water Surface Body Fluctuations Influence Coastal Aquifer Levels - Study of Groundwater Circulation in the Salento Peninsula). /9/ COTECCHIA, V. - Studi e ricerche sulle acque sotterranee e sull'intrusione marina in Puglie (Penisola Salentina), 345 pages, Quademo IRSA 20, 1977. (Studies and Research into Groundwater and Seawater Intrusion in Apulia (the Salento Peninsula, Italy). /lO/PINDER, G.F. & COOPER, H.H. - A Numerical Technique for Calculating the Transient Position of the Saltwater Front, Water Resources Research, 6 (3), 1970. /~~/PINDER, G.F. - A Galerkin- Finite Element Simulation of Groundwater Contamination of Long Island, Ne York, Water Resources Research, 9 (6), 1973.
Mathematics and Computers in Simulation XXIV (1982) 152-160 North-Holland Publishing Company
PRELIMINARY RESULTS FROM A METHODOLOGY GROUNDWATER RESOURCES MANAGEMENT
M. BENEDINI, G. GIULIANO,
B. SIRANGELO
STUDY IN THE FIELD OF
and S. TROISI
Water Research Institute, National Research Council, Rome, Italy
Increasing water demand, conflicts among utilizations, water-body pollution have given a considerable impulse to the application of systems analysis techniques also to the field of groundwater resources management. In line with this, some mathematical simulation and optimization models taking into account the various aspects of management have been prepared. The first of these models, known as IDROSIM, is actually a package including various flow models developed both for the problem of parameter identification (distribution of transmissivity and porosity parameters, etc.) and for the so called direct problem. The second model, called IASO, deals with modifications in water quality as a function of pollution from urban, industrial or agricultural sources. The third model, called IDROPROG, sets out to apply linear programming techniques to the aquifer in order to optimize water resource utilization. Finally, problems of combined management of both groundwater and surface water are investigated by means of a fourth model, dubbed IDRINTEGRO. All these models are currently being tested on a coastal aquifer extending over an area of 600 km* on the Ionian side of the Puglia region of Italy. This paper makes a comprehensive survey of the principles on which the models are based, and provides further details on the IDROSIM simulation model.
1. INTRODUCTION There is a growing trend in groundwater studies to seek more efficient ways of utilizing aquifers. Many different aspects need to be considered in selecting operative criteria for planning the use of ground water resources. "Nonphysical" aspects such as the economic, social and regulatory implications may also be taken into account. The systems approach to water resources management has been applied to the field of surface water for a long time now 1 and world-wide research now under way tends to suggest that similar methods could be profitably extended to the field of groundwater. However, before going into the application of such techniques to groundwater,some peculiarities of aquifer systems need to be discussed. The physical properties of the water bearing material, i.e. its permeability, porosity, storage, etc., differ at virtually every point of the aquifer. Each of these properties affects the aquifer's hydraulic behaviour with relevant points in many management problems such as
0378-4754/82/0000-0000/$02.75
those connected with pollutant dispersion and propagation of impulses in the piezometric surface. As a result, these characteristics are extremely important in decisions as to how to use the groundwater resources. The same tools for describing aquifer behavior both quantitatively and qualitatively are required also in the analysis of the above mentioned non-physical contingencies. For example, by representing how an aquifer is influenced by pumping it is possible to evaluate the amount of water that can be drawn from an aquifer; as well as its quality.This in turn makes it possible to draw up the cost estimates needed to evaluate the economics of using that particular resource. 2. NEW METHODS FOR GROUNDWATER RESOURCES MANAGEMENT With the objective of developing a new methodology suitable for groundwater resources management in the most general context, four mathematical models have been worked out.
0 1982 IMACS/North-Holland
M. Benedini et al. /Methodology
study of groundwater
2.1 The first of these,dubbed IDROSIM,is actually a package incorporating models previously developed for solving both "parameter identification problem" and the change of piezometric surface over a certain period of time, under different initial and boundary conditions. Classical studies show that flow field in a non homogeneous, isotropic aquifer can be expressed mathematically by a second order partial differential equation, which, for a confined aquifer bearing in mind Darcy's Law and matrix compressibility features, can be written as:
where:K S h q
= = = =
permeability (LT-1) specific storage coefficient (L-l) piezometric head (L) unit flow rate drawn off from the acquifer or fed in at time t (LT-1)
A similar expression also holds true for an unconfined aquifer, providing the storage coefficient S can be replaced by the effective porosity, E., and variations of head are small with respect to the saturated thickness of the aquifer. 2.2 The second model named IASO deals with the dispersion of a pollutant in an aquifer (in a two-dimensional case) according to the following fundamental expression, the assumption being made that the pollutant does not react chemically with the solid matrix:
resources management
153
Starting from an assessment of the various phenomenological aspects /3/, computer programs are now being developed for the integration of equation (2) on the basis of a groundwater velocity map produced by the IDROSIM model. 2.3 The third model i.e. IDROPROG, embraces mathematical prograsiming procedures capable of determining that piezometric surface pattern which may be considered as "optimal" in relation to some pre-determined targets of groundwater withdrawal. As is known, formulating a mathematical programming model entails maximizing a function of the problem variables, i.e., the so-called "objective function", respecting suitable boundary conditions for the same variables. These conditions are expressed by other functions called "constraints". Of all the programming methods available, linear programming appears to be one of the most useful. In this method both the objective function and the constraints are expressed as linear functions of the variables. In Darcy conditions /4/, considering the aquifer as being discretized into a certain number of grid elements, a head value "h" can be attributed to the i-th grid element. Any sort of limitations on the use of groundwater can then be reduced to the constraints applied to the head, hi2
h. 1
with h. determined "exogenously" depending on the patticular conditions set for the i-th grid element. Sometimes the limitation may also be imposed on the gradient, in which case the constraint can be expressed as
where:c = pollutant concentration D = dispersion coefficient V x = seepage velocity components V i Y This model takes into account the way a pollutant disperses into an aquifer with reference to the determination of the dispersion coefficient and the groundwater velocity.
h
i
-hk"
Iik
for every pair of grid elements, i and k, with a given Iik. Under stationary flow conditions, hydraulic head can be expressed in the form of n hi = hi0 -
2
i diaij
1 where h. is the head corresponding to the initial con'8. itions while Q is the constant flow i
M. Benedini et al. /Methodology
154
study of groundwater
rate at well i. The "a..II term is an "influence coefficient, denoting ia t e effect caused on the i-th grid element by the water draw-off at the j-th grid element. Establishing
resources management
The corresponding costs of Q and QI, C and C can be assumed to be linear Punctions 0!?Q. I This means that the cost for District i can be expressed as:
that
pi
= cpi + CIi
b. = h. -l-Y 1 10 io with a predetermined h. , b. expresses the piezometric level drop ini?he :-th grid element due to aquifer tapping. Thence the final form of the linear constraints is given by
Qi aij 5 r
bi
i
Considering flow rate Q. as "production" and b. as being the upper limit, the above expressioni shows how withdrawal rates should be determined. The objective function may be expressed as: max
Consequently, considering all the districts in the area, the following total cost function will be:
The "manager's" problem will then be to determine the "best" combination between groundwater extracted inside the area and of the exogenous water so as to minimize overall cost function 0 taking into account the boundaries of Q and QI. P Once the minimization will lead to:
District 2
1 -yl=
has been completed, this
Q&Qpmaxl+ qlQI,
'i
subject to the abovementioned
constraints.
2.4 The fourth model, i.e. IDRINTEGRO, deals with problems of groundwater management in the framework of a multi-objective programing scheme. In order to clarify the matter let us suppose that an area, which is under a regional authority's jurisdiction for water resources management, may be divided up into n water districts. For each district monthly demand requirements may be defined that are inclusive of agricultural, urban, industrial and recreational uses. Separately, by means of simulation and programming models the "manager" can decide on the maximum amount of groundwater to be withdrawn in every district. He can also determine the maximum withdrawal rate (Q max) for the whole acquifer falling under his gurisdiction. The deficit, i.e. the difference between the amount of groundwater available monthly and the demand for the different uses, will have to be compensated month-by-month by bringing in water from external sources (QI), e.g. long distance supply, desalination, waste water treatment and so forth.
CX Q P2 pmax2+
District 2 -Cp,=
=
. . . District n -y
n
=
... o( Q n2 pmaxn +
72912
... %lQIn
It is thus possible to indicate how much water must be pumped from the aquifer month-by-month and district-by-district and how much must be imported from outside districts with adequate groundwater supplies. Furthermore, this system enables the "manager" to fix taxation levels for the volume of water required by the users of every district. Our attention will now be focussed on the mathematical model which is more relevant to the subjects of this Working Conference, i.e. the IDROSIM model.
3. THE IDROSIM MODEL This model is designed to simulate piezometric head variations over a certain period of time for a coastal aquifer subject to exploitation and seawater intrusion. IDROSIM is actually a cluster of models which
M. Benedini et al. /Methodology
study of groundwater resources management
perform different calculations. The entire model development may be broken up into two phases, (Fig. 1). Phase 1 - Preparation of a tape containing bulky input data and - Processing of the basic data fed in (data includes aquifer pattern, boundary conditions, transmissivity, storage coefficient and undisturbed piezometric levels). Phase 2 - Access to the tape by other programs controlled by card inputs with final printout. Phase 1 is further sub-divided units:
into several
a) input of aquifer pattern data via Program DATACFG,
r-‘--'--
Fig. 1 - IDROSIM-Program
System Structure.
155
b) input of Transmissivity (T) distribution via the ITMC Identification Program or via Direct Assignment Program ASSIGN T, c) input of Storage Coefficient (S) distribution via the IMMAG Identification Program or via Direct Assignment Program ASSIGN S, d) calculation of ,an Optimal Relaxation Factor (characteristic of every aquifer) allowing the number of iterations in calculating a piezometric surface for pre-set conditions to be minimized. This is achieved using the FINDWE Program, e) calculation of piezometric surface for undisturbed conditions via Program SOLVEHO. Phase 2, on the other hand, envisages the following units: a) elaboration of piezometric head patterns of the aquifer under well pumping rates either for a steady rate of seawater intrusion,via
-. -.-_
-.
(s= card input; c= print output; n= tape output and input).
156
M. Benedini et al. /Methodology
study of groundwater
Program PIEZO, or a varying rate of seawater intrusion, via Program PINDER, b) calculation of Aquifer Influence Coefficient via Program CD 1, of tidal effects on the aquifer C) calculation using either a unidimensional approximation method, via Program DH 1, or a two-dimensional approximation method, via Program DH 2. Of the first phase programs, ITMC is designed to determine transmissivity values over a rectangular area on the assumption that piezometric head data are available for some points within that area. The ITMC Program originates from Equation (1) and assumes that the right hand side is negligible, thus giving:
(3)
resources management
The IMMAG Program sets out to solve the problem of identifying Storage Coefficient S (or effective porosity) in a non-homogeneous aquifer starting from Equation (l), Transmissivity T (ITMCor ASSIGN T - derived) and the in-aquifer tide propagation being known. Identifying S therefore consists of two distinct phases,namely working out piezometric head oscillations on the basis of field measurements, establishing the minimum of the error function between calculated and measured patterns. Here, too, a part of the program extrapolates the results obtained for a part of the aquifer to the adjacent areas. ASSIGN S program uses the same extrapolation method as the ASSIGN T program wherever storage coefficient data are obtained from field data. Before describing the other programs, Equation (1) for steady flow conditions may conveniently be recalled:
From the above equation the function of T can be obtained by the "Method of Characteristics", whenever the function h(x,y) is known and is continuous with its first and second partial derivatives for all points of the aquifer. However, the piezometric head is usually known only in very few points so it becomes necessary to interpolate discrete data in order to obtain an analytical form of h. The "Cubic Natural Spline" was chosen as interpolating function as it has the above mentioned analytical properties (continuity of function and of first and second derivatives) and at the same time gives a good . approximation of both the function h and its derivatives. Consequently, identification of transmissivity values consists of two quite distinct aspects, namely working out the interpolating Spline func tion and determining transmissivity values along characteristic lines. In other words, the problem consists in (a) determining how Parameter T is distributed over a certain area, (b) picking out some physical plausibility constraints for a part of the aquifer /7/, (c) extrapolating these data for the remaining part.
(1’)
The solution approach, based on finite differences, transforms the problem of solving Equation (1') into a problem of solving a system of linear algebraic equations of the general form: A.x = b where: A= x= b=
coefficients matrix unknown vector of heads known vector related to boundary conditions
The i-th equation of System (5) can be expressed by writing a flow rate continuity equation for the i-th grid element of the discretized aquifer. Referring to Fig. 2, the following
__i_____’ i __
w
‘______~__ I N __
c
E
__
__
s These extrapolation methods are also used for the ASSIGN T Program wherever there is a distribution of direct transmissivity data for the several parts of the aquifer.
(5)
__i._.._ @
I _____~_
Fig. '2- Discretizing
Grid
M. Benedini et al. /Methodology
study of groundwater
resources management
1.57
The PINDER Program may be considered as an offholds:
T~~(h~-h~) + TCW(hC-hW)+ TCE(hC-hE)+TCS(hC-hS)= =
(6)
qC
where; T
CN
T CW T CE T
CS
= (TN + TC)/2 = (TW + TC)/Z = (TE + TC)/2
shoot of the PIEZO Program as it may deal with fresh and seawater flows. Most authors have tackled the problem from an analytical point of view, while PINDER and COOPER /lo/ propose a different approach which, although somewhat elaborate, seems sufficiently all-embracing. The problem can be reduced to a system of three equations, the first of which describes fresh water flow, the second concerns saline water flow and the third expresses the relationship between the interface and the fresh and saline water piezometric heads. After some manipulations, the following final equations can be obtained:
= (TS + T&/Z A - For saline water
Reordering, we have,
=2q
C
Moving along the aquifer grid row-by-row, n equations may Be obtained where n indicates the total number of grid elements. Whenever one of the grid elements falls on the boundary, the corresponding term will be expressed as a known quantity.
d%(X-Xi
;
y-Yj
\
)
(8)
/
B - For fresh water
To have a sufficiently detailed numerical solution, discretization in real cases calls for a large system of linear equations (5), often in the order of thousands. Classical direct methodology (elimination, triangular decomposition, etc.) was discarded due to features of coefficient matrix, the possibility of being limited to approximate solutions, and some operative considerations to implement the algorithm on a computer, Iterative calculation methods (Jacobi, Gauss-Seidel, S.O.R., etc.) were preferred. The IDROSEM model utilizes the S.O.R. iterative method and so, in conjunction with the FINDW Program, it can determine the value of the Optimal Relaxation Factor, w , for the aquifer being studied. The Program System can then go on to solve Equation (5) using the SOLVEHO Program where no pumping occurs and the PIEZO Program where it does. The latter program also enables seawater intrusion to be calculated if the sea water conditions can be assumed to be stationary.
(9) where
I/“’ = I
htd’cl
,T_
b’“’ = ),“I+
zi
b’d’
As usual, A C the water water,E the dary of the
is the sea le~el,B the ground level., table,D the piezometric of saline interface and F an impervious bounsaline water field fl,
M. Benedini et al. /Methodology
158
study of groundwater resources management
and 6'= R = b = h = K = S = p =
transmission in a simplified unidimensional approach /8/ or in a more realistic two-dimensio&approach /7/.
withdrawal flow rate per unit area recharging flow rate per unit area Dirac function piezometric head permeability storage coefficient density
Finally, Program CD 1 enables calculating the "influence coefficients" that establish to what extent one pumping well is influenced by another one operating inside the same aquifer. As mentioned in the preceding paragraph, this program acts as a link between the IDROSIM and IDROPROG models with the aim of defining the "influence coefficients" which establish the constraints of the linear programming model.
(9)
6
= porosity
Subscripts x, y and z represent components in their respective directions. Subscripts s and d indicate saline and fresh water, respectively (Fig. 3). The PINDER Program discretizes Equations (8) and (9), and, once integrated by means of the finite difference method, it provides a picture of the piezometric surfaces of both fresh and saline water as well as of the interface in a coastal aquifer at different utilization rates flf. If, on the other hand, an overall picture is required of how a coastal aquifer is affected at its borders, for example, to what extent groundwater and sea are connected, then Programs DH 1 and DH 2 are used to solve the problem of tidal
4. APPLYING NEW METHODS TO GROUNDWATER RESOURCES MANAGEMENT The purpose of this paper is restricted to the presentation of methodologies and is not concerned with any numerical results of their applications. The latter will be dealt with in separate papers in preparation. The described methods are applied o a coastal 5 aquifer extending over some 600 km on the Ionian side (Fig. 4 of Italy's Salento peninsula which was chosen as a case study. The aquifer in question floats on a body of continent-invading seawater. Aquifer thickness tends progressively to zero in proximity of the
I -:_: -
_.-__--____---_ __ ~-_---__-_._- - __ --
_-_---
:.,_
. . . .
Fig. 3 - Fresh water/saline
.
..
water interface in a coastal aquifer.
-
-
-_:
__
M. Benedini et al. /Methodology
GU
study of groundwater
L 0
159
resources management
F
F
TARANTO
i:
Fig. 4 - Location of the investigated aquifer.
coast line where the permeability of the outcropping rocks allows a free seaward outflow of the groundwater. Studies already carried out /9/ on the dynamics of the diffusion zone between the aquifer and underlying seawater provide the necessary data for establishing the aquifer pattern /ll/. The area chosen is bounded upstream by the water divide separating the flow of groundwater towards the Ionian or the Adriatic Seas. This divide is not, however, permanently fixed but is subject to lateral shifts of up to about one kilometer due to the difference in level between the two seas, differences that at some times last even for long periods, To the northwest, the boundary ?s tdentified by the appearance of another bordering hydrogeolo-
gical system (the Murgia system), i.e. by the sudden variation of certain significant features such as piezometric gradients and heads, water flow rate, etc. To the southeast however, a conventional boundary has been adopted taking into account the flow of groundwater entering the zone from the bordering area. In a joint effort with Bari University, Applied Geology and Geotechnics Institute, several datataking campaigns have been carried out on the aquifer to obtain the necessary inputs for the four models dealt with in this paper. The datataking consisted basically in obtaining: - piezometric heads in the different selected water points - both undisturbed and subject to pumping or tide effects,
160
M. Benedini et al. /Methodology
study of groundwater
- local concentration of inorganic chemicals and specific water quality indicators, by means of on site sampling and specific laboratory analysis. - concentration of bacteria, e,g, total coliforms, fecal coliforms and fecal streptococchi. The aquifer was divided up into a grid of 16,000 elements, each measuring 200m x 2OOm, on which pumping and discharge points were also located. ACKNOWLEDGEMENTS The authors particularly wish to thank the technical assistants B. Soria and A. De Fazio for their contribution in the computer work, as well as Mr. R, Jacovelli and the people at SVIM S.p. A. for their fruitful help in preparing the computer programs.
REFERENCES /l/ QUADERNO IRSA 30 - Un esperrmento metodologico di PSanificazione e Gestione delle risorse idriche in un grande bacino idrografico, Roma, 1976, 232 pages. (A Methodology Experiment in the Planning and Managing of the Water Resources of a Large Hydrographic Basic) /2/ CASTELLANO, L. - TROISI, S,: Alcune considerazioni sulla dispersione di un inquinante in una falda sotterranea. Rapport0 tecni‘co IRSA R/79, 1980. /3/ BENEDINI, M. & TROISI, S. - Procedimenti nuo vi per lo studio della gestione delle risorse idriche sotterranee. Geol. Appl. e Idrogeol. Vol. XII, Parte 1, pgs. 79-109, 1977. (New Procedures for the Study of Groundwater Resources Management). 141 SCHWARZ, J. -I Linear Models for Groundwater Management. Jour. of Hydro. 28 (1976), pgs. 377-397. /5/ YU, W & HAINES, Y.Y. - Multi-level Optimization for Conjunctive Use of Groundwater, Water Resources Research, Vol. 10, 11, 4 August, 1974, pgs. 625-635. /6/ CAROTENUTO, L.; DI PILLO, G.; RAICONI, C. & TROISI, S. - Modelli matematici e identificazione parametrica di falde acquifere. Quaderno IRSA 41 pgs.: 96, 1979. (Mathematical Models and the Identification of Aquifer Parameters). /7/ SIRANGELO, B. & TROISI, S. - Studio metodologico di nuove tecniche di gestione delle risorse idriche sotterranee - Applicazione ad
resources management
una falda reale, Quaderno IRSA 50, 1980. (Methodology Study for New Groundwater Management Techniques - A Case Study). /8/ MAGRI, G. & TROISI, S. - Sulla influenza delle fluttuazioni di specchi d'acqua sui livelli delle falde costiere - Applicazione allo studio della circolazione idrica sotterranea nella penisola salentina, Geol. APP~. e Idrogeol., vol. IV, pgs. 25-40,1966. (On Row Free Water Surface Body Fluctuations Influence Coastal Aquifer Levels - Study of Groundwater Circulation in the Salento Peninsula). /9/ COTECCHIA, V. - Studi e ricerche sulle acque sotterranee e sull'intrusione marina in Puglie (Penisola Salentina), 345 pages, Quademo IRSA 20, 1977. (Studies and Research into Groundwater and Seawater Intrusion in Apulia (the Salento Peninsula, Italy). /lO/PINDER, G.F. & COOPER, H.H. - A Numerical Technique for Calculating the Transient Position of the Saltwater Front, Water Resources Research, 6 (3), 1970. /~~/PINDER, G.F. - A Galerkin- Finite Element Simulation of Groundwater Contamination of Long Island, Ne York, Water Resources Research, 9 (6), 1973.