A Stochastic Optimization Approach in the Design of an Aquifer ...

Water Resources Management 13: 335–351, 1999. © 2000 Kluwer Academic Publishers. Printed in the Netherlands.

335

A Stochastic Optimization Approach in the Design of an Aquifer Remediation under Hydrogeologic Uncertainty Y. A. MYLOPOULOS1, N. THEODOSIOU1 and N. A. MYLOPOULOS2 1 Division of Hydraulics and Environmental Engineering, Department of Civil Engineering,

Aristotle University of Thessaloniki, GR 540 06, Thessaloniki, Greece 2 Department of Civil Engineering, University of Thessaly, GR 383 34, Volos, Greece

(Received: 10 August 1998; in final form: 29 September 1999) Abstract. A stochastic optimization approach is presented for the remediation design of a contaminated aquifer with limited hydrogeologic information. Stochastic simulation using the Monte Carlo technique, produces a series of equally probable realisations of the spatially varying random hydraulic conductivity field. The stochastic flow and transport simulation model is coupled, using the response matrix approach, with a nonlinear optimization algorithm. The whole process is integrated into an algorithm which is effectively applied in the case study of the Kalamaria aquifer, Chalkidiki, Greece. The stochastic optimization procedure is followed by a reliability analysis, giving useful information to the decision makers concerning the effectiveness of the optimal results. Key words: aquifer management, groundwater decontamination, hydrogeologic uncertainty, reliability analysis, stochastic optimization.

1. Introduction Groundwater contamination has become a severe environmental risk not only in the developed, but also in the developing countries all over the world, imposing in our days a limiting factor to the use of one of the most important freshwater resources. This explains the increased public concern about groundwater quality, as well as the high interest of the governmental authorities and the scientific community in protecting and remediating groundwater resources. Aquifer remediation design is a typical engineering procedure in which the technical objectives, that is the plume stabilisation followed by the clean-up of the groundwater in the long run, must be met under a set of technical, environmental, legal, and/or political constraints. Therefore the design takes the form of selecting among a series of alternative solutions which are both technically sound and costeffective. This means that the solution to be selected should satisfy the constraints and the objective of the management problem either in an optimal or in a ‘best’ possible way. Management models, combining groundwater flow and transport simulation algorithms with optimization techniques, have been successfully ap-

336

Y. A. MYLOPOULOS ET AL.

plied for the evaluation of the alternative remediation scenarios and selection of the optimal schemes (Gorelick et al., 1984; Lefcoff and Gorelick, 1985; Ahfeld et al., 1988; Mylopoulos et al., 1991). In all the above mentioned works, it was assumed that the aquifers under study were fully characterised as far as the knowledge of the hydrogeologic and the hydraulic parameters is concerned, and therefore they all adopted the deterministic approach in the design of the alternative decontamination schemes. In recent years, however, it has been recognised that the uncertainties emerging from the heterogeneity of the hydrogeologic environment, in combination with the associated scarcity of information and data, play a very critical role in the engineering design of groundwater remediation schemes. Among various sources of uncertainty the one related with the spatial distribution of the hydraulic conductivity field is the most important. Therefore stochastic groundwater flow and transport simulation, accounting for aquifer parameter uncertainty, has been widely used as a more realistic approach to groundwater resources remediation design (Gelhar, 1986; Mylopoulos, N., 1994). Two different approaches of groundwater quality management under uncertainty have been presented in the literature, each having its advantages, as well advantages as drawbacks: the stochastic optimization and the decision analysis approach. The decision analysis approach can be characterised as better suited to riskbased engineering design than the stochastic optimization one, as its structure permits the inclusion of all costs, benefits, as well as risks of failure due to uncertainties, into a single objective function. This can be achieved by easily coupling a hydrogeologic uncertainty model with a decision one. The risk-cost-benefit objective function developing in this way allows for immediate comparisons and a final decision making, among a number of well specific, real alternative design schemes (Massman and Freeze, 1987; Covello, 1987; Freeze et al., 1990; Massman et al., 1991; Mylopoulos et al., 1994; Latinopoulos et al., 1997). The fact, however, that decision analysis involves the determination of the best from a small set of previously specified alternatives, is the main drawback of this methodology compared to the stochastic optimization. By incorporating the effects of the uncertainty in the set of the constraints and not in the objective function, the latter method involves determination of a globally optimal solution. This main drawback of the stochastic optimization approach, has been addressed by many researchers, who were concerned of the method’s computational complexities, especially in case of large-scale contamination problems. In the work of Gorelick (1990), a comprehensive analysis on the conceptual and practical issues of the overall problem of combined groundwater simulation and nonlinear optimization was presented, while the value of quantifying the uncertainty and the need to incorporate a concept of reliability were emphasised. In this work, a stochastic optimization approach will be presented and applied in the engineering remediation design of the contaminated aquifer of Kalamaria,

A STOCHASTIC OPTIMIZATION APPROACH

337

Figure 1. Chalkidiki peninsula and the Kalamaria plain.

located in the northern part of Greece, in Chalkidiki peninsula. Information concerning the spatially varying hydraulic conductivity parameter of the aquifer is limited, leading to a high degree of hydrogeologic uncertainty. Stochastic simulation using the Monte Carlo technique, produces a series of equally probable realisations of the random hydraulic conductivity field. Using the response matrix approach, the stochastic simulation model is coupled with a nonlinear optimization algorithm. The stochastic optimization model is applied multiply for many realisations of the random field (200 in this case study), as well as for one overall approach, subjected simultaneously to the constraints corresponding to the results of all the simulated realisations. The process is integrated into an algorithm which is developed and applied in the case study of the Kalamaria aquifer. The simplicity in the development and the efficiency in the application in a real case study, are the algorithm’s main advantages. The stochastic optimization procedure is followed by a reliability study, giving useful information to the decision maker concerning the effectiveness of the optimal results.

2. Description of the Problem 2.1.

THE AQUIFER UNDER STUDY

The plain of Kalamaria is situated at the southwestern part of the Chalkidiki peninsula in Northern Greece, as shown in Figure 1. It constitutes the main agricultural region of the Chalkidiki County and as a result of this, it is intensively cultivated and irrigated. During the summer a large number of tourists visit the area. Taking into consideration the high water demand for both domestic and agricultural use, the importance of the protection of the underlying aquifer is obvious since it is the sole water resource in the area. The particular aquifer has been studied by Latinopoulos et al. (1996).

338


Figure 2. The aquifer under study with the contaminant plume.

Figure 3. The hydraulic head distribution in the Kalamaria aquifer.


339

Figure 4. The thickness distribution in the Kalamaria aquifer.

The part of the aquifer that has been contaminated is shown in Figure 2. It covers an area equal to 90 km2 (10 × 9 km), while the area of the contaminant plume is less than 1 km2 . The hydraulic head distribution of the whole aquifer, as well as of its part under study, is shown in Figure 3. Direction of groundwater flow is from the top to the bottom of the figures. The aquifer thickness distribution, ranging from 20 to 60 m, is shown in Figure 4. The aquifer is bounded by no-flow boundaries along the left and right sides, by constant flux boundary along its top side (Neumann type boundary condition), while its bottom side is bounded by seawater, having thus a constant hydraulic head equal to 0 m. The effective porosity of the aquifer is taken as n = 0.30. Information about the hydrogeology of the aquifer, especially on spatial variability of hydraulic conductivity, is limited. Assuming that all other aquifer characteristics can be estimated with good approximation by the existing measurements and historic data available, the stochastic approach is limited to hydraulic conductivity only. The random hydraulic conductivity field is assumed to follow log-normal distribution, with expected mean value of hydraulic conductivity Km = 1 × 10−4 m s−1 , standard deviation σ k = 1 × 10−4 m s−1 and correlation length λ = 500 m.

340


These values were estimated using the limited field measurements of hydraulic conductivity on the one hand and our knowledge about the aquifer material on the other. 2.2.

FORMULATION OF THE REMEDIATION PROBLEM

Due to an accident in an underground industrial pipeline, a spill of waste was released in the region. After leakage the waste passed through the upper geological formations of the aquifer and managed to reach the groundwater, resulting to the contaminant plume shown in Figure 2. The polluted water migrates along with local groundwater flow downwards, threatening the groundwater resource with contamination for a long time-period. The concept of the aquifer decontamination is based on the operation of two pumping wells placed inside the plume area. The locations of the pumping wells are considered to be fixed. The primary purpose of installing the two wells is to stabilise the plume and assure that contamination will not be extended further on the aquifer. In the long run abstraction of the wells should decontaminate the groundwater. Through the proposed stochastic optimization procedure, the optimal pumping rates of the wells will be determined. The hydraulic gradient control is performed at a number of check points at the plume boundary. At these points the reversion of the flow field is required in order to assure that the groundwater hydraulic gradient is always towards the interior of the plume. 3. The Stochastic Simulation Model 3.1.

THE GOVERNING EQUATIONS

The governing equation for the two-dimensional steady state flow in a horizontal confined aquifer is ∂ ∂h Tij = W, i, j = 1, 2, (1) (∂xi ) (∂xj ) where h is the hydraulic head, Tij is the transmissivity tensor (Tij = bKij , where b is the aquifer thickness and Kij is the hydraulic conductivity tensor), W is the source-sink term and xi , xj are the spatial coordinates. The governing equation for reactive single species solute transport in a twodimensional horizontal aquifer is described by ∂(bc) ∂ ∂c c0 W ∂ R bDij − (bcvi ) − = , i, j = 1, 2, (2) ∂(t) ∂xi ∂xj ∂xi n where c is the solute concentration, R is the retardation factor, Dij the hydrodynamic dispersion tensor, vi is the groundwater velocity vector, c0 is the solute concentration in fluid source s−1 ink, n is the effective porosity and t is time.


341

In order to solve the above Equations (1) and (2), the USGS MOC (Method Of Characteristics) model of two-dimensional solute transport and dispersion in groundwater (Konikow et al., 1978), was chosen. This well-known model is a very efficient numerical code with a structure allowing for reliable simulation of realworld problems in both the deterministic and stochastic contexts.

3.2.

THE STOCHASTIC SIMULATION MODEL

In the present study it is assumed that the uncertainty of the hydraulic conductivity is the only source of uncertainty and, therefore, the stochastic approach is focused upon the spatial distribution of this parameter. The quantification of the uncertainty has traditionally been handled by employing a stochastic model of groundwater flow and solute transport. It is generally accepted that the Monte Carlo technique is very suitable for doing this, (Massmann et al., 1987; Freeze, 1990; Mylopoulos et al., 1994; Latinopoulos, 1997). Within the context of the stochastic approach the true spatial distribution of the hydraulic conductivity is considered to be a single realisation of the spatially random field of this heterogeneous variable. By invoking the assumptions of stationarity and ergodicity, the moments of the random field can be inferred using information and/or in situ measurements from the single realisation, i.e., the spatial distribution of the hydraulic conductivity of the aquifer under study (Gelhar, 1986). A typical assumption made for most aquifers is that stationarity is restricted to the first two moments, a fact which guarantees also a stationary covariance. Consequently the mean of the random filed is constant in space and the covariance between two locations, xi and xj , is a function of their separation vector, 1xij , only. As it is widely accepted, (Hoeksema et al., 1985) the distribution of hydraulic conductivity can be assumed to be lognormal. The random log-hydraulic conductivity field is thus statistically characterised by three parameters, namely µY , σY2 and λY , that is by the mean, the covariance and the correlation length. In the present application an isotropic exponential covariance is assumed: Cov[Yi , Yj ] = σY2 exp{−|1xij |/λY },

(3)

where Yi − log Ki is a natural logarithm of the hydraulic conductivity at any point xi , Cov[·] is the covariance of the random variable and |1xij | is the distance between two arbitrary points xi and xj in the field. The Monte Carlo procedure used in this study consists of the following steps: (a) Given the three statistical parameters of the stationary, correlated, 2-D random field, Yi , a number of realisations is generated by applying the method of turning bands (Mantoglou et al., 1982; Mylopoulos N., 1994; Latinopoulos et al., 1997).

342


(b) Each hydraulic conductivity realisation is used as input data to the MOC model together with the set of the rest of parameters, boundary condition, etc., which are considered deterministic in all simulations. The random hydraulic conductivity field is assumed to follow lognormal distribution with expected mean value of hydraulic conductivity as previously mentioned Km = 1 × 10−4 m s−1 , standard deviation σk = 1 × 10−4 m s−1 and correlation length λ = 500 m. To perform the numerical simulations with the MOC model the aquifer is descretized into 1440 square cells, with a discretization of 250 × 250 m2 . Using Monte Carlo method, a number of realisations of the random field of hydraulic conductivity has been produced (200 for this case study), resulting to an equal number of relevant probable aquifer simulations. 4. The Stochastic Optimization Model 4.1.

COUPLING SIMULATION WITH OPTIMIZATION MODELS

In the design of large–scale plume stabilisation problems, groundwater flow simulation models must be combined with optimization techniques. The two main methods used to combine the simulation and the optimization models are the ‘embedding method’ and the ‘response matrix approach’. In the embedding method, finite differences or finite elements approximations of the flow equations are treated as part of the constraints of the optimization model. Although such management models have been successfully developed, according to Gorelick (1983) there are no large-scale applications of this method, due to numerical difficulties that are likely to arise. In the response matrix method, an external flow simulation model is used to develop coefficients, each of which relates unit ‘stresses’ at one location of the aquifer, with their impacts at another. The assemblage of these coefficients, which is called a response matrix, is included in the management model as a substitute of the simulation model. In the current study, the response matrix method was selected to combine the simulation and the optimization models. According to the decontamination scenario, two pumping wells are to assure the stabilisation of the plume, through controlling the plume at a number of check points at the plume boundary. Thus, for the formulation of the response matrix, the correlation between the pumping rates of the two wells and the drop of the groundwater head at the check points, is needed. The procedure for the development of the response matrix for the aquifer of Kalamaria is as follows: for each of the 200 different realisations of the random hydraulic conductivity field, the simulation model is applied for three different scenarios: (a) unit pumping rate at well 1 only, (b) unit pumping rate at well 2 only, and (c) no pumping at all.


343

From all simulations, values of hydraulic heads at the check points and at the well sites are recorded. Linearity between pumping rates and head drawdowns, permits the use of superposition principle and applies the response matrix methodology. Head drawdowns for scenario (a) can be computed by subtracting the hydraulic heads that resulted from this scenario, from the respective heads, at the same points, of scenario (c). The same process can be repeated for the calculation of head drawdowns for scenario (b) , respectively. The coefficients derived from these subtractions express the correlation between groundwater head drawdown and pumping rates at the wells. The composition of these coefficients represents the response matrix of the aquifer. Through this procedure 200 response matrices were developed corresponding to the 200 different hydraulic conductivity distributions of the aquifer. It must be noted, that the value of the pumping rates that is used as ‘unit pumping rate’, is not exactly 1 m3 s−1 but it is actually 0.0566 m3 s−1 . The reason for selecting this value, is that it is closer to the expected value of the pumping rates and thus it implements a much smaller numerical error to the calculations.

4.2.

FORMULATION OF THE OPTIMIZATION MODEL

The basic characteristic of the proposed procedure is the assumption that the stabilisation of the plume can be achieved through the development of a series of equally probable optimization problems, aiming to the hydraulic gradient control along the boundary of the contaminant plume. This technique was introduced by Lefkoff and Gorelick (1986), who presented a quadratic programming formulation using groundwater velocity constraints assuming that pumping costs are a quadratic function of pumping rates. In their work Mylopoulos Y. et al. (1991), extended this approach to a sequential optimization procedure during which the hydraulic gradient control was performed over a time varying plume. The required simultaneous solution of the flow and contaminant transport equations under the framework of an optimization problem, constitutes a highly nonlinear mathematical programming problem, which is very difficult to solve. So the management problem was formulated, as described in the following paragraphs, in such a way as to maximise the simplicity and the efficiency of the algorithm. 4.2.1. The Decision Variables The pumping rates of the two wells, operating inside the contaminant plume for the sake of its stabilisation and removal, were considered as the decision variables of the management problem.

344


4.2.2. The Objective Function One of the most common and realistic measures for the effectiveness of the decontamination system, is the minimisation of the total operating cost. This is achieved by employing a quadratic objective function, expressed as min

M X

qm 1hm ,

(4.1)

m=1

where M is the total number of wells, 1hm is the drawdown at well m due to pumping, and qm is the respective pumping rate at the same well. It is noted that the above objective function is quadratic, in respect to the decision variables q, because the drawdown is in principle a linear function of the pumping rate. This quadratic function represents the energy cost of pumping below the water level (Lefkoff and Gorelick, 1986; Latinopoulos et al., 1994). 4.2.3. The Constraints The linear constraints used in the formulation of the optimization model are of two types: (a) Constraints which are designed to keep the flow confined in the whole aquifer. Thus, at each well location m, lower bounds hu are imposed upon the hydraulic head in order to control excessive drawdowns and assure the under pressure functioning of the confined aquifer: hm > hu ,

m = 1, 2.

(4.2)

(b) Constraints that force the hydraulic gradient at a number of control locations to turn in toward of the plume. This is accomplished by selecting pairs of points that are located close, and at each side, of the boundary of the plume, and imposing constraints of the form: h1k − h2k > h∗k ,

k = 1 to 10,

(4.3)

where points 1 and 2 are located outside and inside the plume respectively. The right hand side term is a user defined lower bound on the difference between the heads at points 1 and 2, ensuring a marginal value for the hydraulic gradient. This value is defined as the product of the minimum gradient required at the specific control location and the distance between the two points. These constraints were applied over ten pairs of check points located around the downstream boundary of the plume.

345


4.3.

INTEGRATING THE PROCEDURE THROUGH AN ALGORITHM

Due to the large number of iterations needed to solve all 200 management problems, a computer program that could combine and automate the procedure, was developed. The program executed the following steps: (1) Read the output files produced from the application of the turning bands method describing the hydraulic conductivity distribution over the aquifer. (2) Produced the input files necessary for the application of the flow simulation model. (3) Executed the flow simulation model for the three scenarios described in paragraph 4.1 (unit pumping rate at well1, unit pumping rate at well 2, no pumping). (4) From the output files produced from the application of the flow simulation model, the computer program recognised the elements that expressed the correlation of the pumping rates and the hydraulic head drop at the wells and at the check points. (5) Produced the response matrix of the aquifer. (6) Using the response matrix, it generated the input files for the optimization model based on the objective function (3.1) and the set of constraints (3.2) and (3.3). The optimization model used was the well-known program MINOS (Murtagh and Saunders, 1987). (7) Executed MINOS. (8) From the output files of MINOS it recognised the results of the management model (the pumping rates at the two wells) as well as the characterisation by MINOS of the status of the final solution (if it was an optimal solution, a nearly optimal solution, or an infeasible problem). (9) The computer program, organised these results in a final output file in order to have a complete picture of the procedure. (10) The procedure was repeated automatically for all 200 realisations of the permeability coefficients distribution fields. It must be noted that for the solution of the optimization problem except from the necessary input files, MINOS required a subroutine describing the nonlinear objective function. In order to accomplish this, during step 6, the algorithm produced the computer code for this subroutine, compiled it and linked it to the MINOS without any user interference. 5. Results and Discussion 5.1.

THE RESULTS FROM THE

200

INDEPENDENT OPTIMIZATION PROBLEMS

The results, i.e., the optimal solutions derived from the procedure described above for the 200 independent realisations of the simulation model, are presented in this

346


Table I. Characteristic values of the decision variables

Maximum value Minimum value Mean value

Pumping rate at well 1 (m3 s−1 )


Total pumping rate (m3 s−1 )

Objective function (m m3 s−1 )

0.313 0.000 0.101

0.257 0.000 0.092

0.414 0.063 0.193

18.781 1.790 6.534

Figure 5. Variations of the total pumping rate (y-axis) for the 200 optimization problems. The x-axis denotes the consecutive number of each problem.

paragraph. Table I illustrates some characteristic values after statistical analysis of the 200 optimal results, in the form of the minimum, maximum and mean optimal values of the decision variables of the problem. The following figures present a general view of the results. Figure 5 presents the variation of the total pumping rates (i.e., the sum of the pumping rates of the two wells), as derived from the solution of the 200 independent optimization problems. Figure 6 shows the corresponding variation of the objective function and of the operating cost that it represents. In these two figures, one can see the magnitude

Figure 6. Variation of the value of the objective function (y-axis) for the 200 optimization problems. The x-axis denotes the consecutive number of each problem.

347


Table II. Results from the overall approach

Overall approach



Total pumping rate (m3 s−1 )

0.276

0.209

0.485

of the influence that the spatial distribution of the hydraulic conductivity field has upon the optimal pumping rates of the clean-up wells. 5.2.

THE RESULTS FROM THE OVERALL APPROACH

At the second stage of the procedure, an overall optimization problem was solved, which was simultaneously subjected to all the constraints derived from the 200 realisations of the distribution of the permeability coefficient. This optimization problem represents the so called ‘worst case’, as it involves no probabilities of failure at all, by assuring satisfaction of the constraints in any possible realisation of the random field. In this way instead of solving independently 200 optimization problems, each subjected to 10 constraints of type (3.3) and 2 constraints of type (3.2), one big overall optimization problem was solved, subjected to 2000 (= 200·10) constraints of type (3.3) and 400 (= 200·2) constraints of type (3.2). In reality the single optimization problem that was solved through this procedure is facing the worst conditions among the conditions of the 200 independent problems. The results of this overall approach are shown in Table II. Comparing the above results with the characteristic results presented in Table I, one can observe that the pumping rates of wells 1 and 2 of the overall solution do not exceed the maximum values of the 200 problems even though the total pumping rate does. 5.3.

RELIABILITY ANALYSIS

In order to estimate the effectiveness of the results of the optimization procedure, a reliability analysis has been conducted. During this analysis, for each of all the 200 optimal pairs of pumping rates produced by the repeated application of the optimization model, their effectiveness in stabilising the plume also under other possible realisations of the random field is investigated. In other words, the effectiveness of each optimal solution is evaluated through its ability to satisfy the constraints of other optimization models, except of the one from which they have been derived. This means that for each pair of optimal pumping rates, all the 200 simulations corresponding to the 200 realisations of the random field have to be

348


Figure 7. Reliability curve for the optimal results. The x-axis denotes the consecutive number of each pair of optimal pumping rates, while on y-axis the number of successful simulations for each pair is shown. The order is ascending with regard to the number of successful approaches (or to the y-axis).

Figure 8. Correlation between the total pumping rate (x-axis) and the reliability of the optimal results. The y-axis denotes the number of successful simulations for each rate. It is noted that one value of the rates of the x-axis may correspond to more than one optimization problems.

executed, in order to check in which of them the plume stabilisation is achieved. The satisfaction of all the constraints of each problem is considered as a successful approach, while the violation of even one constraint is considered a failure. Figure 7 presents graphically the results of the above reliability investigation, in ascending order, with regard to the number of successful approaches. As it can be seen, there are solutions with a very low degree of reliability – even with only one successful approach, obviously corresponding to the problem from which they were derived. On the other hand, there are other solutions that have a very high degree of reliability. The maximum degree of reliability recorded was 197 successful approaches. No solution was found with reliability of 100%, except of course of the overall solution, which was designed to be able to satisfy all the 200 realisations. In Figure 8 the total pumping rate corresponding to the 200 optimal pairs of pumping rates is presented. The optimal solutions are presented in the same as-


349

cending order with regard to their reliability, as in Figure 7. As it can be derived, no correlation can be made between the degree of reliability on one hand, and the total pumping rate on the other. This means that a large value of the total pumping rate does not necessarily lead to a high degree of reliability and vice versa. This becomes true as there are cases, where although two optimal solutions have the same degree of reliability, the difference in their total pumping rate is up to 30%. This happens because, as a result of the hydraulic conductivity’s spatial variability, it is not the sum, but the distribution of the pumping rates in the two wells that counts in the remediation design of the aquifer.

6. Conclusions A stochastic nonlinear optimization algorithm was developed and applied for the remediation design of a contaminated aquifer, with limited information and data concerning the spatially varying hydraulic conductivity parameter. Stochastic simulation using the Monte Carlo technique, produced a series of equally probable realisations of the random hydraulic conductivity field. The stochastic flow and transport simulation model was coupled, using the response matrix approach, with a nonlinear optimization algorithm. The optimization model aimed to determine the optimal pumping rates of two wells, under a set of technical and environmental constraints. The latter assure the stabilisation of the contaminant plume by controlling the hydraulic gradient along its boundary, as well as the maintenance of the under pressure operation of the confined aquifer. The stochastic optimization model was applied multiply for all the realisations of the random field (200 in this case study). The model was applied once more, for the overall approach, which was simultaneously subjected to the sum of the constraints of the 200 independent problems. The whole process was integrated into an algorithm which was developed and applied in the case study of the Kalamaria aquifer, in Chalkidiki, Greece. The simplicity in the development and the efficiency in the application, are the main advantages of the presented algorithm, which proved to be able to produce accurate results even in large scale groundwater remediation problems. Finally, 200 optimal solutions were derived, corresponding to the equal in number realisations of the random hydraulic conductivity field, and a statistical analysis helped in their study and presentation. The reliability analysis that followed provides the decision makers with useful information concerning the effectiveness of the optimal results. The reliability of the optimal solutions was evaluated through their ability to satisfy the constraints of as more as possible optimization problems corresponding to possible realisations of the random field of hydraulic conductivity. The analysis showed that the reliability of the optimal solutions range from 0.5% up to 98.5%, except of course of the ‘conservative’ overall solution, which was designed to be 100% successful. Another important conclusion of the analysis was that the increase of the total pumping rate of the clean-up wells, does not necessarily lead to a higher degree

350


of reliability, and vice versa. This remark shows the magnitude of the influence that the spatial variation of the hydraulic conductivity parameter may have upon the optimal remediation schemes, as the important issue is not the total energy spent for the remediation only, but also the distribution of the pumping rates to the various wells. The stochastic optimization methodology presented in this work, can be part of an integrated process, where the optimal pumping rates of the clean-up wells can be used for the formulation of the alternative scenarios, that will be finally evaluated by the decision analysis approach (Mylopoulos N., 1996). This can be an effective way of combining the main advantages of the two methodologies, namely the determination of globally optimal solutions offered by the stochastic optimization, and the inclusion of all costs, benefits, as well as risks of failure due to uncertainties, into a single objective function, offered by the decision analysis.

References Ahfeld, D. P., Mulvey, J. M. and Pinder, G. F.: 1988, Contaminated groundwater remediation design using simulation, optimization and sensitivity theory 1. Model Development, Water Resour. Res. 24(3), 431–441. Freeze, A., Massmann, J., Smith, L., Sperling, T. and James, B.: 1990, Hydrogeological Decision Analysis 1. A Framework, Ground Water 28(5), 738–766. Gelhar, L. W.: 1986, Stochastic subsurface hydrology. From theory to applications, Water Resour. Res. 22(9), 135S–145S. Gorelick, S. M.: 1983, A review of distributed parameter groundwater management modeling methods, Water Resources Research, 19(2), 305–319. Gorelick, S. M., Voss, C. I. and Gill, P. E.: 1984, Aquifer reclamation design: The use of contaminant transport simulation coupled with nonlinear programming, Water Resour. Res. 20(4), 415–427. Gorelick, S. M.: 1990, Large-scale non linear deterministic and stochastic optimization: Formulations involving simulation of subsurface contamination, Math. Prog. 48, 19–39. Hoeksema, R. J. and Kitanidis, P. K.: 1985, Analysis of the spatial structure of properties of selected aquifers, Water Resour. Res. 21(4), 563–572. Konikow L. F. and Bredehoeft, J. D.: 1978, Computer model of two-dimensional solute transport and dispersion in ground water, Techn. Water Resour. Invest., U.S. Geol. Surv., Book 7, Ch. C2. Latinopoulos, P., Theodossiou, N., Mylopoulos, Y. and Mylopoulos, N.: 1994, A sensitivity analysis and parametric study for the evaluation of the optimal management of a contaminated aquifer, Water Resour. Manage. 8, 11–31. Latinopoulos, P., Xefteris, A., Anastasiadis, P. and Theodossiou, N.: 1996, Application of modeling techniques in a groundwater nitrate pollution problem, In: Proc. of 6th International Conference on Hydraulic Engineering Software, HYDROSOFT 96, Penang, Malaysia, pp. 101–110. Latinopoulos, P. D., Mylopoulos, N. A. and Mylopoulos, Y. A.: 1997, Risk-based decision analysis in the design of water supply projects, Water Resour. Manage. 11, 263–281. Lefkoff, L. J. and Gorelick, S. M.: 1986, Design and cost analysis of rapid aquifer restoration systems using flow simulation and quadratic programming, Ground Water 24(6), pp. 777–790. Mantoglou, A. and Wilson, J. L.: 1982, The turning bands method for simulation of random fields using line generation by a spectral method, Water Resour. Res. 18(5), 1379–1394. Massmann, J. and Freeze, A.: 1987, Groundwater contamination from waste management sites: The interaction between risk-based engineering and regulatory policy. 1. Methodology, Water Resour. Res. 23(2), 351–367.


351

Mylopoulos, N. A.: 1994, Stochastic simulation and decision modeling in groundwater resources management under risk conditions, PhD Thesis, Dept. of Civil Engineering, Aristotle University of Thessaloniki, Greece, p. 201 (in Greek). Mylopoulos, Y. A., Latinopoulos P. D. and Theodosiou, N.: 1991, A combined use of simulation and optimization techniques in the solution of aquifer restoration problems, In: L. C. Wrobel and C. A. Brebbia (eds), Water Pollution: Modelling, Measuring and Prediction, Computational Mechanics Publications, Southampton, pp. 59–72. Mylopoulos, Y. A., Mylopoulos, N. A. and Latinopoulos, P. D.: 1994, Risk-based decision modelling in the engineering design of groundwater pollution problems, In: Tsakiris and Santos (eds), Advances in Water Resources Technology and Management, Balkema, Rotterdam, pp. 269–274. Mylopoulos, Y. A.: 1996, Decision analysis in groundwater resources management, results of a research project, Ministry of Research and Technology, Athens, Greece, p. 291 (in Greek). Murtagh, B. A. and Saunders, M. A.: 1987, MINOS 5.1 User’s Guide, Tech. Rep. SOL 83-20R, Dept. Oper. Res., Stanford Univ., Stanford, CA. Wagner, B. J. and Gorelick, S. M.: 1989, Reliable aquifer remediation in the presence of spatially variable hydraulic conductivity: From data to design, Water Resour. Res. 25(10), 2211–2225.