Multi-objective optimization for free-phase LNAPL recovery using evolutionary computation algorithms

Abstract

A nonlinear, multi-objective optimization methodology is presented that seeks to maximize free product recovery of light non-aqueous phase liquids (LNAPLs) while minimizing operation cost, by introducing the novel concept of optimal alternating pumping and resting periods. This process allows more oil to flow towards the extraction wells, ensuring maximum free product removal at the end of the remediation period with minimum groundwater extraction. The methodology presented here combines FEHM (Finite Element Heat and Mass transfer code), a multiphase groundwater model that simulates LNAPL transport, with three evolutionary algorithms: the genetic algorithm (GA), the differential evolution (DE) algorithm and the particle swarm optimization (PSO) algorithm. The proposed optimal free-phase recovery strategy was tested using data from a field site, located near Athens, Greece. The PSO and DE solutions were very similar, while that provided by the GA was inferior, although the computation time was roughly the same for all algorithms. One of the most efficient algorithms (PSO) was chosen to approximate the optimal Pareto front, a method that provides multiple options to decision makers. When the optimal strategy is implemented, although a significant amount of LNAPL free product is captured, a spreading of the LNAPL plume occurs.

Editor Z.W. Kundzewicz; Associate editor L. See

Citation Dokou, Z. and Karatzas, G.P., 2013. Multi-objective optimization for free-phase LNAPL recovery using evolutionary computation algorithms. Hydrological Sciences Journal, 58 (3), 671–685.

Résumé

Nous présentons une méthode d'optimisation non-linéaire et multi-objectif visant à maximiser la récupération de liquides légers en phase non aqueuse (LLPNA) tout en minimisant les coûts d'exploitation, en introduisant une nouvelle stratégie de périodes optimales de pompage et de pause alternées. Ce procédé permet l’écoulement de davantage de pétrole vers les puits d'extraction, et assure une élimination maximum de produit à la fin de la période de dépollution tout en minimisant l'extraction d'eaux souterraines. La méthodologie présentée ici combine FEHM, un modèle d’écoulement souterrain multiphasique simulant le transport des LLPNA, avec trois algorithmes d’évolution: l'algorithme génétique (AG), l'algorithme à évolution différentielle (ED) et l'optimisation par essaims particulaires (OEP). Cette stratégie de récupération a été testée à l'aide des données d'un site situé près d'Athènes, en Grèce. Les solutions de l'OEP et l'ED étaient très similaires tandis que celle de l'AG était moins intéressante, alors que le temps d'exécution était à peu près le même pour tous les algorithmes. L'algorithme le plus performant (OEP) a été choisi pour se rapprocher du front de Pareto optimal, méthode qui offre plusieurs options aux décideurs. Lorsque la stratégie optimale est mise en œuvre, même si une quantité importante de LLPNA est récupérée, on observe cependant la propagation d'un panache de LLPNA.

Key words:

• free product LNAPL recovery
• multi-objective optimization
• Pareto front
• genetic algorithm
• differential evolution
• particle swarm optimization

Mots clefs:

• récupération de LLPNA libres
• optimisation multi-objectif
• front de Pareto
• algorithme génétique
• algorithme à évolution différentielle
• optimisation par essaims particulaires

INTRODUCTION

Light non-aqueous phase liquids (LNAPLs) such as petroleum hydrocarbons (oil, gasoline, diesel etc.) and organic solvents are responsible for the contamination of groundwater reserves worldwide. Problems associated with these contaminants arise from disposal dumps, petroleum refineries, leaking storage tanks and pipelines and accidental spills (Chrysikopoulos et al. 2003). There exist numerous sites worldwide (Hinchee et al. 1991, Lu et al. 1999, Bass et al. 2000, Banks et al. 2003, Khaitan et al. 2006, Gidarakos and Aivalioti 2007, Qin et al. 2008) that have been contaminated by LNAPLs and the number is growing rapidly.

LNAPL migration is a complex process that is affected by various parameters, the most important of them being gravity forces, which cause vertical migration, and capillary forces, which cause horizontal spreading. During the vertical migration of LNAPL some will remain behind, trapped in the vadose zone (residual saturation). When LNAPL reaches the capillary zone it accumulates, creating an oil table that rests directly above the water table, often causing a depression on the water table due to its weight. Finally, some of the LNAPL partitions into the groundwater through dissolution causing long-term groundwater contamination (Qin et al. 2009).

Successful remediation of contaminated sites is essential for the protection of human health, but it can prove to be a costly and time consuming task. To this end, researchers have directed their efforts towards developing tools that can improve the time-efficiency and cost-effectiveness of groundwater remediation strategies. Such tools are algorithms that couple groundwater contaminant transport simulation models with optimization techniques (Dokou and Karatzas 2010).

Since the LNAPL remediation process involves different tasks, such as clean-up of the residual LNAPL phase, free-product recovery and removal of the dissolved phase, various techniques have been developed to accomplish the optimization of each task. In recent years, most researchers have focused on optimizing the removal of LNAPLs dissolved in groundwater using the pump-and-treat method, which has been the most commonly-used remediation strategy (Gorelick et al. 1984, Ahlfeld et al. 1988, Chang et al. 1992, Culver and Shoemaker 1992, Karatzas and Pinder 1993, 1996, Johnson and Rogers 1995, Rizzo and Dougherty 1996, Papadopoulou et al. 2003, 2007, Hilton and Culver 2005, Babbar and Minisker 2006).

Several studies have also been devoted to optimizing bioremediation designs using a variety of methods, such as: analytical derivatives (Minsker and Shoemaker 1996, 1998); evolutionary algorithms (binary-coded genetic algorithm, real-coded genetic algorithm and derandomized evolution strategy); direct search methods and derivative-based optimization methods (Yoon and Shoemaker 1999); multiscale derivatives (Liu and Minisker 2002, 2004); a genetic algorithm integrated with constrained differential dynamic programming (Hsiao and Chang 2002); and a combination of stepwise cluster analysis, nonlinear optimization and artificial neural networks (Huang et al. 2006).

Other researchers have focused on optimizing different remediation techniques, such as bioslurping—using a set of regression submodels describing the system response to water and gas pumping and oil skimming in hydrocarbon recovery (Yen and Chang 2003), bioventing (Diele et al. 2002) and soil vapour extraction (SVE)—using a mixed-integer programming model to determine the optimum number of wells, their locations and pumping rates (Sawyer and Kamakoti 1998).

There are a few studies that have attempted to optimize free product recovery of LNAPLs. Specifically, Cooper et al. (1998) presented a methodology for optimizing free product recovery from a single well that combines simulation, nonlinear regression and optimization, neglecting the economic aspects of the problem. They used MINOS, an optimization algorithm that solves linear or nonlinear problems. In order to solve nonlinear problems it utilizes projected Lagrangian and augmented reduced-gradient algorithms. Qin et al. (2007) coupled numerical modelling, a multivariate regression method and a genetic algorithm to optimize a vacuum-enhanced free product recovery (VFPR) process. Their method included environmental and economic effects, and provided a means to analyse the trade-offs between them. A related work by Qin et al. (2008) optimized a similar process, namely the dual-phase vacuum extraction (DPVE), using a multiphase flow simulator and cluster analysis in conjunction with a genetic algorithm to solve a multi-objective optimization problem.

Optimizing free-phase product recovery is a complex problem that usually involves multiple conflicting objectives, such as the requirement of minimizing the remediation cost while taking into account the environmental aspects (maximization of free product recovery). Thus, multi-objective formulations are of great use when trying to solve such problems. Multi-objective optimization models have been employed by researchers in order to solve similar problems. Ritzel et al. (1994) used a genetic algorithm to solve a multi-objective groundwater pollution containment problem (vector-evaluated GA and Pareto GA). Erickson et al. (2002) used a niched Pareto genetic algorithm for the optimization of a pump-and-treat system that simultaneously minimizes the remedial cost and the contaminant mass which remains at the end of the remediation period. McPhee and Yeh (2006) developed an experimental design-based methodology for groundwater management using multi-objective programming for parameter estimation that was solved by combination of a genetic algorithm and gradient-based optimization techniques.

Recently, there has been a growing interest in using evolutionary algorithms for solving groundwater remediation problems. These models are global, population-based methods that provide an alternative to traditional linear or nonlinear optimization models. The disadvantage of using linear optimization methods is that the linearity assumption is often unrealistic in field applications. While conventional nonlinear methods (classical, dynamic, mixed-integer) have a strong mathematical representation, they are not very attractive for real-world applications due to the large requirements in computational effort related to complex nonlinear manipulations (Qin et al. 2009). Specific to the problems of pump-and-treat, or free product recovery, the functions of various groundwater system components (e.g. system cost) may be discontinuous or highly complicated; thus, the calculation of the related derivatives with respect to the decision variables can be problematic (McKinney and Lin 1994). This renders evolutionary algorithms more attractive for solving this type of problems. The interested reader is referred to Mayer et al. (2002) and Qin et al. (2009) for an extensive literature review of the recent developments associated with optimization techniques applied to site remediation.

As a continuation of the work performed previously, the focus of this paper is to develop a simulation-optimization model that couples a multiphase flow simulation model using three evolutionary algorithms: the genetic algorithm (GA), the differential evolution (DE) algorithm and the particle swarm optimization (PSO) method, taking into account both the environmental and economic aspects of the remediation problem. More specifically, the goal is to achieve maximum LNAPL free product removal at least cost by optimizing both the pumping rates on each well and the remediation time. A concept that has not been previously modelled or optimized is introduced in this work. This concept involves optimal alternating pumping and resting periods, to allow more oil to flow towards the extraction wells, ensuring maximum free product removal at the end of the remediation period with minimum groundwater extraction. The performance of the proposed optimal free-phase recovery algorithms (GA, DE and PSO) is compared using data from a field site contaminated with LNAPLs, located near Athens, Greece. Then, one of the most efficient algorithms (PSO) was chosen in order to approximate the optimal Pareto front as a means of providing insight into the trade-offs between the two objectives, information critical for the decision-making process.

Methodology

Multiphase flow simulator

For the purpose of simulating the LNAPL transport in the subsurface, FEHM (Finite Element Heat and Mass transfer code), a groundwater model developed by the Los Alamos National Laboratory, USA, was used. This model was initially designed to assist in the understanding of flow fields and mass transport in the saturated and unsaturated zones below the area of the Yucca Mountain, which is both hydrologically and geologically complex. Its purpose is to simulate mass transfer for multiphase flow within porous and permeable media, and noncondensible gas flow within porous and permeable media (Zyvoloski et al. 1999).

In this work, isothermal NAPL-water transport was assumed. This assumption is generally valid for shallow subsurface transport, where pressure is practically constant and physicochemical properties are not affected significantly by temperature fluctuations (Chen et al. 2006). The model uses a pressure formulation and solves two conservation equations; one for liquid (free-phase) NAPL and one for liquid water:

(1)

(2)

(3)

The above equations are subject to the following initial conditions:

(4)

where l denotes liquids (l = w for water and l = n for NAPL); S_l are the water or NAPL saturations, with S_w + S_n  = 1; q _l ( x ,t) are the water or NAPL fluxes; F _l ( x ,t) is a source or sink term; _l ( x ,t) is liquid mobility; P _l ( x ,t) is the fluid pressure; ρ_l is fluid density; k ( x ) is the intrinsic permeability of the porous media; k_rl is the water or NAPL relative permeability; μ _l is the liquid dynamic viscosity; P _l0( x ) is the initial pressure in the domain; P _lt ( x ,t) is the prescribed pressure on a Dirichlet boundary segment Γ _A ; Q _l ( x ,t) is the prescribed fluid flux across Neumann boundary segments Γ _B ; g is the gravity vector; n ( x ) is the outward unit vector normal to the boundary Γ _B and ϵ is the porosity. The input to the model consists of an initial description of the fluid pressure as well as media properties. The output consists of the final fluid pressure and the volume fraction of water-NAPL (Zyvoloski et al. 1999).

For the capillary pressure and saturation relationship, the Brooks-Corey model was used based on the work of Reitsma and Kueper (1994), who conducted an experiment on fractured rock in order to measure the relationship between capillary pressure and saturation, and concluded that the Brooks-Corey model best fits the laboratory data. The Brooks-Corey model is described by the following equations:

(5)

(6)

where P _c is the capillary pressure; P _d is the entry pressure; S _e is the normalized water saturation; S _w is the water saturation; λ is a pore size distribution index; and S_r is a curve-fitting parameter representing the irreducible water saturation. The above model was chosen for this particular application because the main geological formation found in the area of interest is fractured limestone.

Optimization problem

Problem formulation

The objective of the proposed optimization algorithm is twofold: maximizing the remediation efficiency (LNAPL free product removal) while minimizing operation cost. The extraction well locations and number were assumed fixed; thus, the capital cost associated with their construction was not included in the objective function. It has been observed in the field that, in order to ensure maximum product removal, the pumping should be followed by a resting period that will enable the LNAPL product to flow towards the extraction wells before pumping is started again. Thus, in the optimization formulation the remediation period comprised two parts: a pumping period and a resting period. Consequently, the decision variables of the problem were the NAPL pumping rates for each well, and the pumping and resting times. The number of pumping–resting alternations was assumed fixed. LNAPL pumping rates instead of total liquid (oil plus water) rates were used in the optimization problem to model the operation of the oil pumps that were used in the field. These pumps were automatically turned off when the entire free product volume present in the well was pumped, and started pumping again when more oil flowed in the well at the end of the resting period. This way, a minimal amount of water is pumped from the aquifer.

The first objective of the optimization problem is associated with the economic aspect of the free product recovery, in this case the operation cost that is directly proportional to the pumping rate and duration, and is described by the following equation:

(7)

where c ₁ is the unit cost of operation; t_i is the duration of each pumping period; n is the number of pumping sub-periods; Q_j the pumping rate of each well; and m is the number of pumping wells.

The second objective involves the environmental considerations of the problem that are represented in this work by the maximization of free product removal from all extracting wells during the entire remediation period, which is described by the following equation:

(8)

where V_ij is the LNAPL free-phase volume recovered during each pumping period i from each extraction well j.

In single-objective optimization problems, the main focus is on the decision variable space, while in a multi-objective concept the interest switches to the objective space. Due to the contradiction of the objective functions, it is not possible to find a single optimal solution for all of them simultaneously (Miettinen et al. 1998). A simple and common approach to overcome this problem is to combine the objectives into a single objective function, assigning weights to each of them according to their relative importance (conventional weighted aggregation):

(9)

where w ₁ and w ₂ are weights that define the relative importance of the two terms of the objective function.

An alternative to optimizing the objective function with fixed weights as presented above, is to determine a set of “Pareto optimal” designs. When the Pareto optimal set is considered, a better understanding of the trade-offs between the two objectives is obtained. Usually, the weights are defined according to the decision maker's preferences. The main characteristic of the points on a Pareto optimal front is that these are dominant, meaning that there exist no other points that have smaller values for both the objectives.

In this work, two tests were done: in the first test, three optimization algorithms were compared (GA, DE and PSO) using fixed weights. In the second test, one of the most efficient algorithms (PSO) was used in order to approximate the Pareto optimal front, performing multiple runs in order to span a range of sets of weights, under the assumption that the summation of the weights should equal 1 (w ₁ + w ₂ = 1). A comparable work was performed by Schaerlaekens et al. (2006) for the remediation of DNAPL contamination using the shuffled complex evolution algorithm (SCE).

Solution algorithms Classical optimization algorithms are not convenient when dealing with multi-objective problems. The motivation behind the choice of the three evolutionary algorithms—the genetic algorithm (GA), the differential evolution (DE) algorithm and the particle swarm optimization (PSO)—for this problem was based on the fact that, unlike traditional gradient-based methods, population-based algorithms are less likely to get trapped in local minima; in addition, they are able to deal with non-differentiable and non-convex functions (Abido 2002).

Evolutionary algorithms are stochastic search methods inspired by natural biological evolution. A fixed-size population (NP) of potential solutions is required to initialize the process. The initial population should cover the entire parameter space and in most cases it is created randomly. The next step involves the evolution of the population using genetic operators such as crossover and mutation in order to find superior solutions. The newly created population is then used in the next iteration (generation) of the algorithm until a stopping criterion terminates the process. The stopping criterion can be either in the form of a maximum number of generations or a satisfactory fitness level. Even if the algorithm has terminated due to a maximum number of generations, a satisfactory solution cannot be always guaranteed (Qin et al. 2009).

Genetic algorithm

During each iteration of a genetic algorithm, individuals from the previous population are selected using some selection scheme, and are combined (crossover) in order to form a new population. Some of those individuals will undergo mutation, which is a random change in one or more of their chromosomes. The number of individuals that will be combined and/or mutated is determined by the crossover and mutation probabilities respectively. The rest of the individuals will enter the new population unchanged (Goldberg 1989).

Traditional genetic algorithms are binary coded. Nevertheless, when applied to real-world problems, real-coded evolutionary algorithms have proved more computationally efficient and easier to implement (Yoon and Shoemaker 2001). For this reason, in this work a real-coded genetic algorithm was used. In addition, special attention was paid, in adopting an “elitist” scheme by replacing the individual with the lowest fitness at each generation with the best individual found up to that generation, in order to ensure the preservation of the best solution throughout the procedure.

Differential evolution algorithm

In DE algorithms (created by Storn and Price 1997), new individuals are created through the process of mutation by adding the weighted difference between two individuals ( x _{r 2,G} , x _{r 3,G} ) and a third ( x _{r 1,G} ), in order to create what is called the mutated vector ( u _i,G+1) according to the following equation:

(10)

where r ₁, r ₂ and r ₃ (i = 1, 2, 3, … , NP) are mutually different random indexes and F is a scaling parameter between 0 and 2.

The mutated vector is combined with a randomly chosen individual of the population, called the target vector x _i,G , in order to produce a trial vector u _ji,G+1, through the process of crossover:

(11)

where randb(j) is the jth evaluation of a uniform random number generator with outcome [0,1] and rnbr(i) is a randomly chosen index that ensures that u _ji,G+1 gets at least one parameter from u ′ _ji,G+1 and D is the dimensional parameter space.

Finally, the fitness of the trial vector is compared to that of the target vector and if it is greater, the trial vector replaces the target vector in the next generation (selection process). During each iteration, every individual has to serve once as the target vector (Storn and Price 1997).

Particle swarm optimization algorithm

The PSO algorithm mimics the behaviour of swarms of animals. In this method, each individual is considered as a particle in multidimensional space with a specific position and velocity that keeps track of the best position it has achieved so far (Eberhart and Kennedy 1995). The collection of particles is called a swarm, a term analogous to the population term in GAs and DEs. At each iteration, a particle moves to a new position in space by adding a velocity to its current position. The velocity term is a random combination of three components: one causing the particle to continue moving in the direction it was moving in the previous iteration (inertia component), one causing the particle to move towards the best position it has ever been in (cognitive component) and a third steering the particle towards the best position of any particle of the entire swarm or in its neighbourhood (social component).

This can be summarized by:

(12)

(13)

where v _i is the velocity of the particle; φ₁ and φ₂ are positive numbers (weights); rand₁ and rand₂ are uniformly-distributed random numbers in the range [0,1]; and p _i and p _g are the best previously recorded positions of the ith particle and of the entire swarm, respectively. During the implementation of the PSO, there are certain parameters that need to be taken into account in order to avoid the “explosion” of the swarm, and to speed convergence. These include the selection of the maximum velocity, the acceleration constant and the constriction factor or inertia constant (del Valle et al. 2008). In the literature, different methods are given to define these parameters. In this work, the maximum velocity is limited using a multiplier that is usually between 0 and 1 as follows:

(14)

where k is the multiplier and x _max is the upper bound of the variable.

The acceleration constant controls the movement of the particle driven by its own experience (φ₁) and the experiences of the other particles of the swarm (φ₂). It has been found empirically that, even when the maximum velocity and the acceleration constant are defined correctly, the particles may still diverge. There are widely-used methods that can overcome this problem: the constriction factor method (Clerc and Kennedy 2002) and the inertia weight method (Shi and Eberhart 1998). The method selected herein is the constriction factor method which updates equation (12)(11) as follows:

(15)

(16)

Although PSO is considered an evolutionary algorithm there exist some differences between PSO and genetic and DE algorithms, which are considered more traditional evolutionary algorithms, but the procedure is analogous. In general, PSO utilizes information exchange only among the particle's own previous experience and the experience of the best particle in the swarm (or swarm neighbourhood), instead of being carried from “parents” to offspring as in evolutionary algorithms. In evolutionary techniques, there are three main operators: selection, crossover and mutation. In PSO, the crossover operator does not exist, but the concept of accelerating a particle towards its previous best position, as well as towards the best overall position or the best position of the particles in the neighbourhood is analogous to the crossover. The PSO approach utilizes a directional position updating operation with a kind of built-in memory that resembles mutation. PSO does not utilize selection, nor the survival of the fittest concept, given that all particles are members of the swarm for the entire run (Eberhart and Shi 1998).

While the genetic and DE algorithms have been successfully applied in the field of groundwater hydrology and contaminant transport by many researchers (Rogers et al. 1995, Yoon and Shoemaker 1999, Erickson et al. 2001, Zhang et al. 2005, Babbar and Minisker 2006, Karterakis et al. 2007, Mayer and Endres 2007, Trichakis et al. 2009, Cardiff et al. 2010, Papadopoulou et al. 2010), there are fewer and more recent applications of the PSO algorithm in these fields (Matott et al. 2006, Gaur et al. 2011, Tian et al. 2011).

Field application and results

LNAPL model

The applicability of the proposed strategy and the relative effectiveness of the three evolutionary algorithms were demonstrated through field application. The study area is located near Athens, Greece. The environmental assessment performed at the site revealed significant hydrocarbon contamination in all three phases (free-phase product, soil vapour and groundwater solutes). Soil, soil-gas, groundwater and oil product samples were collected, and piezometric level measurements were performed in the spring of 2001. Shallow boreholes of a large diameter were drilled in the area in order to collect groundwater and free product samples. The initial investigation detected a free-phase product plume with an area of 150 000 m², an estimated mean product thickness of 0.2 m and a maximum apparent thickness of 1.11 m. According to the measurements, the general groundwater flow is towards the sea. The groundwater flow and transport model constructed for the site focuses on LNAPL free-phase transport and removal.

The geological background of the study area was determined by borehole data, which indicated that the main geological formation encountered in the area is a grey to greyish-white limestone consisting of a fractured upper part that extends 1–6 m below the ground surface. This part constitutes an unconfined aquifer which was vertically discretized in two numerical layers. The horizontal discretization of the study area was implemented using a quadrilateral finite-element mesh consisting of 902 nodes and 832 elements.

The study aquifer is shallow; thus, the physicochemical properties are not affected significantly by temperature fluctuations. This makes the isothermal NAPL-water transport assumption valid for this work. In addition, since only the free product recovery is optimized, it can be assumed that the model results will not be significantly affected by this assumption.

A single porosity model was used that replaces the hydraulic conductivity and porosity values of the fractured system with their locally-averaged, smoothed values. Practically, the fractured formation is represented as an equivalent porous medium by replacing the primary and secondary porosity and hydraulic conductivity distribution with a continuous porous medium with equivalent hydraulic properties (Dokou and Pinder 2011). This method has been widely used in similar applications because of its simplicity compared to the dual porosity (Zyvolovski et al. 2008) or discrete fracture (Dokou and Karatzas 2012) methods, which require a large amount of data and more complex modelling. The influence of fractures on the LNAPL movement can be significant; however, since no such data were available for this particular site, the modelling approach had to be simplified. This field study serves as a demonstration of the applicability of the optimization methodology and the comparison of the three evolutionary algorithms in general; if the proposed optimization methodology is chosen to be applied in the field, a more detailed investigation of the nature of the fracture network at the site will be needed in order to achieve more accurate modelling of the fractured aquifer.

The initial hydraulic head distribution for each layer was obtained by interpolation of field measurements, corrected for the effect of the floating oil product. The correction was performed using the following equation:

(17)

where h_c is the corrected hydraulic head value, h_m the measured hydraulic head, H ₀ the LNAPL thickness measured at the well, ρ₀ the LNAPL density and ρ_w the water density.

The calibration process of the flow field, performed by fine tuning of the flow boundary conditions and using data from nine locations (H1–H9 in Fig. 1), achieved a good fit between modelled hydraulic head values and measured head data (R² = 0.85, Fig. 2). The calibration was performed using measurements taken in May 2001 as initial data and measurements collected in May 2008 as target data.

Fig. 2 Comparison between measured and calibrated hydraulic head values.

Fig. 1 Field site domain, flow boundary conditions and calibrated hydraulic head field.

An initial hydrocarbon free-phase distribution was created by interpolating existing LNAPL thickness field measurements from May 2008 (Fig. 3(a)). The estimated volume of free product in the aquifer is 472.3 m³. The location of the 20 pumping wells that are already installed on site (P1–P10 and M1–M10) and are used for the recovery of the LNAPL free product is also shown in Fig. 3(a) and (b).

Fig. 3 (a) Initial LNAPL distribution and location of extraction wells; (b) Final LNAPL thickness distribution after the PSO optimal pumping strategy has been implemented.

The hydraulic conductivity field of the model area is considered heterogeneous. The main model area has a hydraulic conductivity of 2 × 10⁻³ m s⁻¹ and, additionally, there are three lenses (of the same depth as the rest of the domain) with hydraulic conductivities of 3.1 × 10⁻⁴, 1.9 × 10⁻⁴ and 4.3 × 10⁻⁴ m s⁻¹ around wells M9, M10 and P9, respectively, and one with a value of 1.1 × 10⁻⁴ m s⁻¹ north of well M10. These values were determined by pumping tests performed on site.

Due to the lack of actual capillary pressure–saturation measurements at the study site, laboratory measurements in fractured limestone taken by Reitsma and Kueper (1994) were used herein. The Brooks-Corey type model was found to fit their data more closely as compared to the van Genuchten model. The actual fitting parameters (P_d is represented as equivalent height of water in cm) are presented in Table 1 along with other model parameters, such as the porosity, hydrocarbon density and viscosity that were obtained either from site measurements or the literature.

Table 1 Groundwater flow and LNAPL transport model parameters

Parameter	Value
Equivalent porosity	0.03
Hydrocarbon density	878.7 kg m^−1
Hydrocarbon viscosity	6.04 × 10^−4 Pa s
Number of pumping/resting periods	10
Brooks-Corey model
λ	2.4
Pd	2.9 cm
Sr	0.05

Optimization problem

Evolutionary algorithms evaluation and comparison In the first part of the field application, three optimization algorithms were compared (GA, DE and PSO) using fixed objective function weights (w ₁ = 0.1, w ₂ = 0.9). The calibrated FEHM model was used in conjunction with the optimization algorithms in order to select an optimal remediation strategy for the LNAPL contamination in each case. For all algorithms, a population of 25 individuals (or particles for the PSO) was used in each generation, and a maximum number of 500 iterations was defined as the stopping criterion in all three cases, corresponding to 12 500 calls to the simulation model. The initial population was created randomly in all cases. Previous research (Arifovic 1998, Abido 2002) has shown that the choice of initial values for evolutionary algorithms can affect their performance in the sense that they influence their convergence speed, but not the final optimal solution obtained by the algorithm. The number of iterations in this work was chosen to be large enough to minimize this effect.

Table 2 Test for the selection of the optimal algorithm parameter values

	PSO
k	0.3	0.2	0.4	0.6
Objective value	−154.172	−150.503	−153.901	−152.175
φ1 = φ2	2.05	2.5	3	3.5
Objective value	−154.172	−144.504	−144.646	−143.359
	GA
Pc	0.8	0.6	0.7	0.9
Objective value	−142.872	−142.498	−141.838	−141.752
Pm	0.1	0.05	0.2	0.3
Objective value	−142.872	−142.103	−141.976	−141.513
	DE
Cr	0.8	0.6	0.7	0.9
Objective value	−152.796	−152.549	−152.114	−146.989
F	0.5	0.1	0.8	1.2
Objective value	−152.796	−146.783	−146.566	−143.626

In order to select the parameter values that would produce the best results for the three algorithms, different sets of parameters were tested. For each algorithm two parameters were adjusted: the crossover (P_c ) and mutation probabilities (P_m ) for the GA, the crossover (C_r ) and scaling parameter (F) for the DE and the multiplier for the maximum velocity (k) and acceleration constants (φ₁ = φ₂) for the PSO. The parameters selected for this test were based on initial values taken from the literature (Goldberg 1989, Price et al. 2005, del Valle et al. 2008) and were kept within each parameter bounds. The number of individuals (or particles) was kept at 25, but fewer iterations (100) were performed to reduce the computational effort. The results of this test are summarized in Table 2 with the first column containing the optimal parameter values determined by this procedure: P_c  = 0.8 and P_m  = 0.1 for the GA, C_r  = 0.8 and F = 0.5 for the DE, and k = 0.3 and φ₁ = φ₂ = 2.05 (corresponding to a constriction factor of χ = 0.729) for the PSO. The result for the acceleration constants is consistent with findings of other researchers (del Valle et al. 2008). It should be noted here that, although there is no absolute guarantee that the same parameters will be the optimal when 500 iterations are performed, this test is able to provide good estimates in a reasonable timeframe. All other optimization model parameters are summarized in Table 3. Note that the cost coefficient was selected in order to make the two objectives comparable and does not correspond to the real cost of remediation; it is used only as a “relative” cost.

Table 3 Optimization parameters

Parameter	Value
Population number	25
Maximum number of iterations	500
Pumping rate bounds: Q min, Q max	0 kg s^−1, 0.2 kg s^−1
Time period bounds	0.2–10 d
Cost coefficient: c 1	1 € m^−3
Weights: w 1, w 2	0.1–0.9

Table 4 Optimal pumping strategies for the GA, DE and PSO regarding pumping rates and pumping and resting times

Pumpingwell	Q (kg s^−1)			Time	t (d)
	GA	DE	PSO		GA	DE	PSO
P1	0.06	0.38	0.37	Pump1	5.63	3.71	1.34
P2	0.02	0.16	0.20	Rest1	9.72	6.02	4.38
P3	0.03	0.24	0.29	Pump2	0.61	1.32	0.20
P4	0.01	0.43	0.01	Rest2	9.95	9.80	5.98
P5	0.14	0.11	0.33	Pump3	0.79	0.12	0.25
P6	0.03	0.08	0.38	Rest3	9.78	4.38	2.86
P7	0.10	0.23	0.50	Pump4	1.91	0.10	0.20
P8	0.04	0.49	0.07	Rest4	9.19	8.75	10.00
P9	0.04	0.14	0.07	Pump5	0.57	0.11	1.27
P10	0.01	0.15	0.19	Rest5	9.88	9.20	8.11
M1	0.03	0.01	0.17	Pump6	0.64	0.66	3.24
M2	0.06	0.15	0.35	Rest6	9.89	6.74	5.29
M3	0.08	0.19	0.14	Pump7	0.09	0.53	0.90
M4	0.01	0.37	0.01	Rest7	9.66	8.59	7.37
M5	0.08	0.03	0.17	Pump8	0.90	1.72	0.25
M6	0.01	0.23	0.01	Rest8	8.51	9.07	7.88
M7	0.02	0.04	0.38	Pump9	0.67	0.41	0.20
M8	0.02	0.16	0.01	Rest9	8.68	3.68	9.96
M9	0.01	0.09	0.22	Pump10	1.32	0.33	0.20
M10	0.03	0.16	0.03	Rest10	2.64	3.71	1.34
	Objective function value
	−144.43	−181.55	−180.59	Total	101.3	78.95	71.22

The optimal solutions obtained from the three optimization approaches are presented in Table 4. Although the three algorithms started from a similar initial objective value, the GA's performance was inferior to that of the other two algorithms, converging with a slow rate towards a significantly larger objective value than the other two algorithms. Specifically, the genetic algorithm converged to an objective value of –144.43, while the other two algorithms achieved very similar objective values, that of the DE being slightly smaller (–181.55) than that of the PSO (–180.59). This finding is consistent with the results of other researchers (Abido 2002) who have noted that the GA's inability to control the balance between local and global exploration of the search space (unlike PSO) can lead to premature convergence.

The optimal objective function values for the DE and PSO algorithms are very similar. The optimal value for the DE is slightly better than that of the PSO for the 500 iterations that were chosen as the convergence criterion. However, the PSO converges much faster to a nearly optimal solution and remains always better than the DE for the iterations 180–485, as can be observed in Fig. 4. This means that, if a smaller convergence criterion was used, PSO would have found a significantly better solution than the DE. The free product removal achieved by the GA, DE and PSO algorithms was 171.6, 205.6 and 204.2 m³, corresponding to removal of 36.3%, 43.5% and 42.3%, respectively. The total remediation time according to the optimal strategies chosen by each algorithm was 101.3 d for the GA, 78.95 d for the DE and 71.22 d for the PSO.

Fig. 4 Convergence rate of GA, DE and PSO algorithms.

In Fig. 3(b), the final LNAPL thickness distribution after implementing the optimal pumping strategy of the PSO is presented. Similar figures were produced for the other algorithms, but are not presented for brevity. While the LNAPL thickness is reduced to zero in all 20 pumping wells, as expected, the maximum LNAPL head observed in the contaminated area reaches 0.32 m. The relatively high LNAPL head is observed in the area down gradient from wells M3 and P5. This suggests that when using the existing pumping wells in the area, the contamination cannot be fully contained and some product might escape, moving towards the sea. This is also evident when comparing the initial and final LNAPL distributions (Fig. 3(a) and (b)); while the maximum LNAPL thickness has been reduced from 0.71 to 0.32, a relatively small spreading of the LNAPL plume has occurred, thus the plume's extent has been slightly increased.

Pareto front approximation In the second test, the PSO algorithm, which was found to be the most efficient, together with the DE, for the fixed-weight multi-objective optimization problem studied above, was used in order to approximate the Pareto optimal front of the same optimization problem by means of performing multiple runs. The goal of this test was to analyse the trade-offs between cost and LNAPL removal. Optimization was performed for the following sets of weights: 0–1, 0.1–0.9, 0.3–0.7, 0.5–0.5, 0.7–0.3, 0.9–0.1 and 1–0. The seven Pareto points obtained from this procedure are shown in Fig. 5. Although the original optimization problem involved the maximization of the recovered LNAPL as a second objective, in Fig. 5 the LNAPL that remains in the system after the end of the remediation period is plotted in order for the results to be comparable with previous works (Shaerlaekens et al. 2006). The first and last Pareto points are calculated by single objective optimization since one of the weights is zero in each case. The first point does not take into account the cost and achieves the maximum LNAPL removal, while the last point minimizes cost without taking into account the LNAPL removal; thus the objective values in this case are both zero, corresponding to the “no-action” scenario. The values of the two objectives (cost and remaining LNAPL) ranged significantly from 0 to 95.5 for the cost and from 257.0 to 472.3 m³ for the LNAPL volume remaining in the subsurface. From the graph in Fig. 5, it is evident that a large increase in cost is required in order to remove the last fraction of LNAPL free product. The best alternatives for the decision maker would be to choose one of the next two points on the Pareto front (either the 0.1–0.9 or the 0.3–0.7) that achieve an acceptable level of remediation with a significantly smaller cost, although from a mathematical point of view all Pareto points are equally good and the selection of a strategy should be performed taking into account other criteria, such as the project budget and the environmental regulations.

Fig. 5 Calculated Pareto points.

Conclusions

A methodology for modelling and optimizing LNAPL recovery was presented that included the calibration of a multiphase model using data from a field site near Athens, Greece, the construction of a multi-objective optimization problem that optimized both the pumping rates and pumping and resting times of the LNAPL free product taking into account the cost of remediation and environmental considerations, the comparison of three optimization algorithms in terms of efficiency and the calculation of the Pareto optimal front using one of the two algorithms that were found to be the most efficient.

The solutions obtained using the PSO and DE optimization algorithms were very similar concerning the pumping rates and times and the optimal objective function values. The GA did not perform as well, although the computation time needed to perform the same number of algorithm iterations for all algorithms was roughly the same. Specifically, the genetic algorithm converged to a larger objective value (–144.43) than the other two algorithms, which achieved very similar objective values, with that of the DE being slightly smaller (–181.55) than that of the PSO (–180.59). However, the PSO converges to a nearly optimal solution much faster and remains always better than the DE for a large fraction of the algorithm iterations (iterations 180–485). Thus, it can be said that, if a convergence stopping criterion was used, or if a smaller number of iterations was chosen as a convergence criterion, then the computation time needed by the PSO in order to find a satisfactory solution would be much smaller than by the DE.

The Pareto optimal front, which is a method of analysing the trade-offs between the two objectives, was also calculated using the PSO algorithm. All points on the Pareto front are considered mathematically optimal, but the choice lies in the hands of the decision makers and depends on different parameters that are not taken into account directly by the optimization problem, such as the project budget and the local environmental regulations regarding the maximum allowed contamination levels.

Although a significant amount of LNAPL free product is captured by the optimal strategies (about 43% of the product found in the system for the two most efficient algorithms), the graphical plot of the remaining LNAPL thickness in the study area indicates that, using the existing network of extraction wells, the maximum LNAPL thickness has been reduced from 0.71 to 0.32. However, a small spreading of the LNAPL plume has occurred; thus, the plume's extent has been slightly increased. Consequently, future work should focus on identifying optimal locations for drilling new wells that can satisfactory contain and remove the LNAPL free product in the area.

Acknowledgements

The Los Alamos National Laboratory is gratefully acknowledged, for allowing the use of the FEHM simulator, and especially the code development team leader, George Zyvoloski, for his valuable help with the program.