Stochastic Optimization for an Analytical Model of Saltwater Intrusion in Coastal Aquifers

The present study implements a stochastic optimization technique to optimally manage freshwater pumping from coastal aquifers. Our simulations utilize the well-known sharp interface model for saltwater intrusion in coastal aquifers together with its known analytical solution. The objective is to maximize the total volume of freshwater pumped by the wells from the aquifer while, at the same time, protecting the aquifer from saltwater intrusion. In the direction of dealing with this problem in real time, the ALOPEX stochastic optimization method is used, to optimize the pumping rates of the wells, coupled with a penalty-based strategy that keeps the saltwater front at a safe distance from the wells. Several numerical optimization results, that simulate a known real aquifer case, are presented. The results explore the computational performance of the chosen stochastic optimization method as well as its abilities to manage freshwater pumping in real aquifer environments.


Introduction
The present study Saltwater intrusion in freshwater aquifers is a major concern in coastal areas around the globe. Apart from natural disastrous phenomena, such as earthquakes or floods, human activities and pumping policies nearby the coastline can seriously affect the chemical constitution of water aquifers.
Under predevelopment conditions, groundwater systems are in long-term equilibrium (cf. [1][2][3]). That is, averaged over some period of time, the amount of water entering or recharging the system is approximately equal to the amount of water leaving or discharging from the system. Water discharges and recharges are more or less equal and the total volume of the water within the aquifer remains almost constant. Groundwater levels fluctuate in time over a relatively modest, natural range. Once pumping begins, this equilibrium changes. The main factors for this are the consumptive water use and the irrigation for agricultural needs at nearby areas, that peak during summer months.
People's access to high quality water is a priority for the UN not only due to the climate changes effects, especially severe in coastal regions, but also due to the "unfair" allocation of water resources. The increased worldwide concern to protect groundwater resources in coastal regions is a motivation for the research community to better understand the behavior of complex physical systems and explore new areas of optimization that will improve the evaluation of integrated management solutions as they are absolutely essential to maintain freshwater quality and avoid all the economical implications from salinization of freshwater aquifers. Scientific and engineering community's interest in this problem remains unrelenting for over a century (cf. [4][5][6][7][8] for a review) and has increased over the recent years, as the demand for optimal solutions becomes eminent.
In the present study the sharp-interface model [2] for saltwater intrusion in coastal unconfined aquifers of finite size and its analytical solution (cf. [14]) are used as foundation tools for the optimal management of freshwater pumping by means of stochastic optimization. Our main objective is to produce a simple method able to rapidly maximize the total pumping rate of the aquifer wells, and, at the same time, provide full control over both the safety of the aquifer wells and the pumping policies, as far as salinization is concerned. For this, we adapt ALO-PEX (cf. [28]) stochastic unconstrained optimization procedure coupled by a penalty system, that implements the necessary constraint conditions to avoid saltwater intrusion, to guide and manage the pumping process. The behavior and the computational efficiency of the whole optimization procedure are being investigated through the examination of several pumping and weather-soil conditions on a rectangular shaped aquifer. Our test cases, presented here, approximate a known problem from the literature (cf. [14,15]) which simulates a real aquifer at Vathy area of the Greek island of Kalymnos.

ALOPEX stochastic optimization
ALOPEX (ALgorithm Of Pattern EXtraction [28,29]) is a stochastic iterative process seeking the optimum values of the control variables, in order to maximize an objective (profit, cost, energy) function. It is based on a balanced combination of a bias feedback term and a stochastic term to update the values of the control variables. The bias term uses local correlations between changes in individual control variables as well as changes in the global error measure and has the tendency to move each control variable to the direction which was successful in the past. The stochastic term is a random number, generated for each control variable in each iteration, and provides the opportunity to move each variable against the direction of recent success and avoid local extrema. Balance between the bias and the stochastic term is absolutely necessary as the method reduces to a simple gradient descent in case the stochastic term is omitted, and to a simple random walk in case the bias term is omitted.
ALOPEX was originally devised and introduced in [28] (see also [29]) for the purpose of experimentally determining receptive fields of individual neurons in the visual pathway. Later on, the method was studied in detail and new versions were introduced aiming in accelerating its convergence (cf. [30] and [31] for example). Application areas, where the method was successfully applied, include neurophysiology [28], pattern recognition (cf. [32]), image enhancement in adaptive optics (cf. [33]), cardiovascular or biomechanical applications (cf. [34]) and neural networks (cf. [32,[35][36][37]). Some preliminary and limited results pertaining to the utilization of the ALOPEX II form for coastal aquifers may be found in [38]. The method, through its applications has been compared with known methods, demonstrating success in variety of problems and showing specific advantages over other optimization techniques (cf. [39]).
ALOPEX is an easily implemented stochastic algorithm and presents an alternative to Simulated Annealing [40] (SA). Some of its main features include: • Requires no knowledge of the dynamics of the system, it is essentially gradient-free, and it does not explicitly depend on the functional form of the objective function but rather on estimates of its values in previous steps. This fact alone has a major contribution to the applicability of the algorithm in real time environments and its implementation by use of analog or digital devices (cf. [41]).
• Updates simultaneously all control variables, while SA does not. Each variable is computed using only local computations and independently of the other variables. The computation of the whole process can be carried out completely in parallel (cf. [42]) contributing to its computational efficiency.
• The magnitude of the stochastic term does not depend on the amount it raises or lowers the response, as in SA via the Boltzmann distribution. This makes it easier for the variables to traverse wide gaps between extrema. However, if its magnitude is much higher than the magnitude of the bias term, then the process is mainly driven by randomness, a situation similar to the melting process in SA.
Coastal aquifer management problems are formed as optimization problems with typically nonlinear objectives and constraints. Stochastic, heuristic and adaptive search algorithms, like ALOPEX, SA or Genetic Algorithms, have been developed and suggested for problems with nonlinear objectives and constraints, mainly to overcome the difficulties gradient-based methods are facing with local extrema. Adding to this the ease of implementation and low computational cost of the ALOPEX stochastic algorithm, we believe that it might be potentially useful to the engineering community.

Materials and Methods
As we have already mentioned, the model we consider for the saltwater intrusion problem in coastal unconfined aquifers of finite size, is based on the sharp-interface simplification, which assumes that there is no mixing (transition) zone between freshwater and saltwater inside the aquifer (see Fig 1), and on the Ghyben-Herzberg relation, which at steady state estimates the position of the interface assuming that horizontally floating freshwater floats over static saltwater. The sharp-interface approximation is reasonable in regional scale problems when the transition zone is narrow relative to the scale of the problem (cf. [9,15]). In this situation, the approximation is more appropriate depending on the hydro-geological parameters, the aquifer geometry (especially the layer thickness), the location of the pumping wells and the recharge/pumping scenarios, as it has been recently studied in [9]. For the model aquifer we consider in this work it will be shown, in later sections (in particular in section that investigates validity issues of the optimization results), that there are ranges of the hydro-geological parameters that narrow the width of the transition zone without compromising the stability of the whole computation. In this sense, one may regard the model aquifer as an appropriate one to be studied under the sharp interface and the Ghyben-Herzberg assumptions. In the section that follows, the governing equation of the sharp-interface approach and its analytical solution are described.

Model equations based on a single potential formulation
Consider the unconfined aquifer, depicted in Figs 1 and 2, and observe that the interface between saltwater and freshwater areas is assumed sharp. Freshwater movement from inland towards the coastline repels saltwater towards the sea. When the total discharge rate Q increases due to pumping activities, saltwater toe penetration x τ moves towards the interior of the aquifer, while, at the same time, the stagnation point x s , that separates water flowing into the sea from water flowing into the well, moves towards the coastline (see Fig 2). And whenever Q overpasses a critical value Q c , toe penetration x τ crosses over the stagnation point x s salinizing the well. Clearly, therefore, the detection of the saltwater wedge toes x τ would provide a measure for the distance of saltwater interface from the active pumping sources (wells) and alarm us for the danger of pumping salinized water.
Working towards this direction, Strack (1976, cf. [43]), considering the sharp-interface approximation and the Ghyben-Herzberg relation, as well as the Dupuit hydraulic assumption, derived a two dimensional single-potential model and used it, together with the method of images, to analytically describe the saltwater interface. Since then, many authors have been also used the single potential formulation (cf. [8, 14-16, 19, 21, 44, 45]) to study saltwater intrusion.
To briefly describe Strack's model and fix notation, let us first consider the Ghyben-Herzberg relation: where ρ s and and ρ f denote the density of salt and fresh water, respectively. Moreover, it is well known that the governing equation of steady flow in an unconfined aquifer can be expressed as (cf. [2,14]): where b satisfies Then, Strack's flow potential ϕ = ϕ(x, y), defined by (cf. [43]): is a continuous and smooth function across the boundary between zones 1 and 2, satisfying the differential equation with Dirichlet boundary conditions ϕ = 0 at the coast boundary (x = 0) and Neumann boundary conditions ϕ η = 0 at the no flow boundaries (η is the normal vector). Furthermore, at the location of the interface toes, where ξ = d hence h f = (1 + δ)d from relation Eq (1), there holds If, now, the values of K, N, Q and the boundary conditions are known, the differential Eq (5) can be solved for ϕ(x, y) using analytical or numerical methods. Once ϕ = ϕ(x, y) is determined, the position z = z(x, y) of the interface surface, as well as the piezometric surface h f , can be calculated as a function of ϕ, by ð7Þ and the locus of the toes of the interface area can be determined by solving, for x τ , the nonlinear equation described in relation Eq (6) above.

Analytical solution
Let us now consider the case of homogeneous (constant hydraulic conductivity K) rectangular aquifers of finite size with a fixed head boundary at the coastline, and three other impermeable boundaries, as it is visually demonstrated in Fig 3 that follows. For this case, Mantoglou [14], using the method of images and superposition to incorporate both ambient flow and surface accretion, extended Strack's work and presented a generalized analytical solution of the differential Eq (5). The analytical representation of the flow potential ϕ = ϕ(x, y) included in [14] and used in our simulations is given by: where and Q j denotes the pumping rate (m 3 /day) of the j-th aquifer well w j , j = 1, . . ., M, with coordinates (x j , y j ) and considered to be non-negative (Q j ! 0). We point out that, throughout this paper, the wells are assumed to be numbered in a way such that

Pumping Management
As we have already mentioned in previous sections, our main objective is to maximize the total volume of fresh water pumped by the wells while keeping at the same time the wells safe from saltwater intrusion. To numerically simulate said task, let us define the normalized objective (or profit) function P: and where Q j and Q j denote the minimum and maximum, respectively, pumping rate (m 3 /day) capabilities of the jth-well and may depend on several technical parameters (e.g. pumping equipment, pipes, hoses) as well as on local pumping policies imposed to regulate social, agricultural or industrial needs, while Q A denotes the maximum volume (m 3 /day) of fresh water available for pumping from the aquifer.
Observe that, if we let then we may also write the objective function as where, of course, Q ¼ ½Q 1 . . . Q M T denotes the vector of maximum pumping rates.
It can be, now, easily verified that, as hence @PðQÞ @Q i ¼ 0 for SðQÞ ¼ SðQÞ and sign @PðQÞ the basic properties of the objective function are described by: PðQÞ # for SðQÞ > SðQÞ ð18Þ Therefore, the objective function P(Q) attains its maximum To protect, now, the wells from saltwater intrusion, we couple the objective function P(Q) with a reinforced Toe constraint. For this, consider two circular areas around every well, shown schematically in Fig 4 (cf. [38]). The inner circular area represents the area of the cone of depression, where water-pumping results in an accelerated movement of the water in this area towards the well. The radius of the cone's of depression area is denoted by r c and, in the literature, varies from 200m to 1200m depending on the maximum well pumping rates and the aquifer natural recharge capabilities. The outer circular area represents a safety distance barrier around the cone of depression of radius r c + d s , with d s > 0, to prevent saltwater reaching the cone's of depression area. Then, the Toe constraint, which couples the objective function P, takes the form where (x τ,j , y j ) denotes the saltwater front's point across the j-th aquifer well w j with coordinates (x j , y j ). We remark that it might be necessary, especially for aquifers with more than one barriers at the sea, to consider a small number of adjacent points (e.g. on the circle of radius r c + d s ) to effectively protect the wells from saltwater intrusion. Finally, we would like to point out that the toe constraint condition may also be coupled (cf. [15]) by the potential constraints which force the free surface of the aquifer to be maintained above the sea level. These constraints mainly protect wells located sufficiently far from the coast (cf. [15]). In this work we will adopt them, together with the rest of the constraints, in an attempt to resolve certain stability issues of the optimization results.

ALOPEX for optimized pumping management
Following last section's discussion, pumping management can be expressed in the following nonlinear optimization form: We will also consider the case of reinforcing the above set of constraints by the potential constraint, described in Eq (22), namely and will report separately on its effect on the total pumping rate.

ALOPEX optimization procedure
For the solution of the above problem we consider the employment of the ALOPEX stochastic optimization method along with appropriate penalty conditions to implement the constraints described above. Taking into account the fact that the objective function is already normalized, among the several different variations of ALOPEX stochastic optimization a modified ALO-PEX II version (cf. [31]) is adopted and used in our simulations. In each iteration step of this ALOPEX stochastic process all control variables Q j , j = 1, . . ., M, are updated simultaneously by means of the following vector rule with provided, of course, initial values of the control variables Q (k) for k = 0, 1. Furthermore, c k is a real parameter controlling the amplitude of the feedback term, while g (k) is the noise vector, with values uniformly distributed in an appropriately chosen interval, to provide the necessary agitation needed to drive the process to global extrema avoiding local problems. As we further discuss in Remarks 2 and 3 at the end of this section, the values c k and g ðkÞ j , attained throughout the whole ALOPEX process, are given by and g ðkÞ where γ denotes a small percentage 1% − 2% of Q ðkÀ 1Þ j (here we use γ = 0.015) and X ðkÞ j is a random variable with values uniformly distributed in the interval (−0.5, 1). We remark that, with c k defined as in Eq (26) above, the feedback term of the ALOPEX rule in Eq (24) reduces to c k ΔP (k−1) ΔQ (k−1) = sign(ΔP (k−1) )ΔQ (k−1) , while the noise term corresponds to a few m 3/ day. The next paragraph is being devoted to the development of a simple and effective deterministic frame of rectifications applied to the main ALOPEX step, described in relations Eqs (24) and (25) above, in order to implement the constraints described in Eq (23).

Constraints implementation
In each iteration step k = 2, 3, . . ., the values of the control variables Q ðkÞ j , j = 1, . . ., M, obtained by the ALOPEX relationship in Eq (24), are examined to satisfy the constraints described in Eq (23). In case of violation, an appropriate penalty strategy is implemented by means of the following rectifications: Rectification where 0 μ 1 is a real parameter.
Rectification for global maximum pumping rate constraint violation To ensure no violation of the maximum volume of fresh water available from the aquifer constraint, described in Eq (12), will occur at any time during the optimization procedure, let us assume that k − 1 iteration steps (k ! 2) of the optimization process have been carried out successfully (that is without any constraint violation or appropriately rectified). Furthermore, during the current k iteration step, assuming that all values Q ðkÞ i , i = 1, . . ., p − 1 and 1 p M, have been successfully evaluated, let S ðkÞ pÀ 1 be defined by  (24), the global maximum pumping constraint is being violated wheneverS in which case the value of Q ðkÞ p is being rectified by where ν ! 1 is a real parameter. The ALOPEX stochastic iteration step, described in Eq (24), together with the above three local constraint rectifications, described in Eqs (28)-(32), constitute MODULE I of the proposed technique. Its data flow is depicted in Fig 5, while the corresponding algorithmic implementation is demonstrated in Fig 6 below through the procedural flow chart that follows. One may easily observe the pipelined mode of execution through MODULE I.

Rectification for the toe constraint violation
The strategy of rectifying the pumping rates of the wells in case of toe constraint violation unfolds in two stages. More specifically, let us assume that all values Q ðkÞ j , i = 1, . . ., M have been evaluated by means of the ALOPEX main step Eq (24) and have been further rectified by relations Eqs (28), (29) and (32) whenever a corresponding violation has occurred. Then, for each index p 2 {1, 2, . . ., M} such that a toe constraint violation has occurred pertaining to the corresponding p-well, namely the first stage of rectification is implemented by where 0 < z p < 1 is a real parameter. In case the above correction does not resolve the toe violation problem occurred, and relation Eq (33) still holds, then the second stage of rectification is activated and implemented by where 0 < z i < 1 are real parameters. The two stages of the tow constraint rectification strategy, described above, is implemented algorithmically through MODULE II. Its procedural flow chart and data flow are depicted in Fig 7 that follows.
Rectification for the potential constraint violation A similar simple strategy may also be applied for rectifying the pumping rates of the wells in case of the potential constraint violation. More specifically, let us assume that all values Q ðkÞ j , i = 1, . . ., M have been evaluated by means of the ALOPEX main step Eq (24) and have been Stochastic Optimization for Saltwater Intrusion in Coastal Aquifers further rectified by relations Eqs (28), (29), (32) and (34) where ξ < 1 is a real parameter. At this point we would like to conclude this section by means of the following remarks: Remark 1-Stopping Criterion Besides our primary purpose to demonstrate the asymptotic converge properties of the stochastic constrained optimization strategy, described in the previous paragraphs, we also introduce an effective stopping criterion of the whole iterative process. The stochastic nature of the optimization procedure led us to use a combination of the standard deviation σ Λ of the objective function's values in the last Λ iterations (with Λ relatively small) and the differential of the objective function's mean value μ Λ in a window of the last 2Λ iterations. To be more specific, let us define the windowed mean value μ Λ and standard deviation σ Λ as where k denotes the current iteration and k − 2Λ > 0. Then, for positive small tolerances and δ, whenever the value of the boolean expression is true, where Δμ Λ = jμ Λ − μ 2Λ j with m 2L ¼ 1 L P kÀ L i¼kÀ 2L P ðiÞ , all values of the objective function are in a narrow neighborhood of their mean value while the rate of change of the mean values of the objective function tends to zero. Thus, assuming also that all constrains are satisfied, the optimization process may be effectively terminated. Although outside of the scope of this work, we point out that the above strategy may also be coupled with slope related criteria.
Remark 2-Feedback amplitude The value of the parameter c k which, together with the term ΔP (k−1) in Eq (24), determines the amplitude of the ALOPEX feedback mechanism, owes to honor simplicity and low cost computations, main characteristics of the ALOPEX process. For this, recall relation Eq (24) and observe that, in the absence of noise, one obtains and by extension or a Taylor expansion of S(Q) in Eq (13), we also have allows the feedback term to effectively couple the noise term and guide the random walk and, therefore, accelerates the process considerably. Remark 3-Noise amplitude Certainly, several different strategies may be adopted for the effective determination of an appropriate noise interval. Traditionally, ALOPEX algorithms are combined with constant amplitude noise intervals, such as, for example, where g ðkÞ j denotes the noise term in Eq (24), 0 < γ < 1 represents a small percentage of Q j , and X ðkÞ j denotes a random variable with values uniformly distributed in an interval G, such as G ¼ ðÀ 0:5 ; 1Þ. The above rule has been successfully used in our simulations. In this work, we slightly deviate from the above rule and combine ALOPEX with varying amplitude noise interval. More precisely, we considered the noise interval which takes into consideration the current values of the pumping parameters and implements an "annealing" strategy for the compromised front (that is, closest to the coast) wells. This varying amplitude noise interval proved to be at least as effective as the constant amplitude one and adopted in our simulations.

Numerical simulations
Among the large number of simulations we performed to test the stochastic pumping management process, described in the previous sections, in this section we include the results of selected numerical simulations involving mainly low number of pumping wells, that better illustrate the qualitative characteristics of the sea front as well as the behavior of pumping management process. The test case aquifer, used in all experiments, approximates a known real case aquifer at Vathi area of Kalymnos island in Greece (cf. [14,15]) and is considered to be a rectangular aquifer with one barrier at the sea and three impermeable barriers on the other sides (see Fig  3). Its characteristics are reported in the Table 1 that follows.

Stochastic Optimization for Saltwater Intrusion in Coastal Aquifers
The values of the pumping management parameters are included in Table 2 and have been kept constant in all numerical simulations. For their determination we have implemented a uniform strategy of ±5% correction. Namely, the value of any control variable, obtained after the enforcement of the penalty described in relations Eqs (28) or (32), is increased by 5%, while, after the enforcement of the penalties described in relations Eqs (29) or (34)- (35) or (37), is decreased by 5%. Therefore, the values of the corresponding parameters in the penalty relations Eqs (28)-(37) are 1±0.05 as reported in Table 2 that follows. The value of 5% has been chosen to maintain balance between two different targets. Namely, on one hand it has to be large enough to resolve the associated constraint violations and make the penalty system effective and, on the other hand, it has to be adequately small to allow the ALOPEX stochastic algorithm to guide the optimization process and converge within small fluctuations. Furthermore, as it pertains to the stopping criterion parameters in Eq (39), we remark that the optimization process terminates whenever all values of the objective function, during the last Λ iterations fluctuate around their mean value by a distance of at most = 10 −2 , an indication that the system converges at least locally to a value, and at the same time two consecutive windowed mean values of the objective function have at least three identical decimal points (δ = 10 −4 ).
All experiments were conducted on an Intel quad-core i7 2.6Mhz PC with 8 Gb DDR3 RAM using MATLAB environment.

Simulation Cases and Optimization Profiles
The results from two groups of simulation cases are reported in this section. The first group includes four artificial test cases characterized by two or three pumping wells, located at positions (in meters) (x i , y i ) (see Table 3), which, besides from illustrating the qualitative characteristics of the sea front, are used as well to test the optimization algorithm's ability to sense small differences of the data profile and reflect them in the proposed pumping plan. The second group includes two test cases from the literature, characterized by five and eleven (cf. [14,15] respectively) pumping wells (see Table 4), used for comparison purposes.
In our simulations we consider two different optimization profiles depending on whether or not we make use of the potential constraints. Namely, • Profile 1: The ALOPEX is coupled by the constraints described in Eq (23) • Profile 2: The ALOPEX is coupled by the constraints described in Eqs (22) and (23). Profile 1 is the constraint profile adapted in all simulations used to demonstrate the performance of the ALOPEX optimization algorithm. Profile 1 is being reinforced, in the sequel, by the potential constraints, hence evolves to Profile 2, in order to report on its effect in resolving stability issues. To avoid the singularities of ϕ(x, y) at the well coordinates (x i , y i ) in implementing the potential constraint Eq (22), the x-coordinate x i is being approximated by x i %x i ¼ x i À 1, as we have already mentioned in Eq (36).

ALOPEX Pumping Management Performance
The main results pertaining to the performance of the ALOPEX pumping management process for all simulations cases, when constraint Profile 1 is implemented, are summarized in Tables  5-7 and Figs 8-10. To be more precise, aiming in testing ALOPEX asymptotic behavior, Table 5 includes the optimal values of the main parameters involved in the pumping management process obtained after 500 iterations, that is the values of the parameters that maximize the objective function max 1 k 500 PðQ ðkÞ Þ and do not violate any of the constraints described in Eq (23).
To test the effectiveness of the stopping criterion, introduced in the previous section, Table 6 includes the optimal values of the main parameters obtained within k SC iterations, that is the values of the parameters that maximize the objective function max 1 k k SC PðQ ðkÞ Þ and do not violate any of the constraints described in Eq (23), where k SC denotes the smallest iteration number for which the stopping criterion relation Eq (39) holds.
Observe that the optimal values of P(Q (k) ) and S(Q (k) ) included in Table 6 are practically the same to the corresponding values included in Table 5. Furthermore, the computational time included in Table 6 corresponds only to a small fraction of the computational time included in Table 5. These properties indicate the effectiveness of the stopping criterion used to terminate the ALOPEX optimization process.
The fact that the differences, between the asymptotic and stopping criterion optimal values, are small also yields small differences in the associated contours of the sharp interface between saltwater and freshwater at the bottom of the aquifer as depicted in Fig 8. Table 7, that comes along with Fig 8, contains the x-coordinates x τ,i of the toe critical points (x τ,i ,y i ) that correspond to the asymptotic optimal values (included in Table 5) of all simulation cases. The coordinates of the critical points that correspond to the stopping criterion optimal values are very close to the ones reported in Table 7, as may easily be verified by inspection of Fig 8, and therefore have been omitted.
An observation that merits highlighting is that the optimal position of the saltwater-freshwater interface, produced by the ALOPEX pumping management process and reported in Fig  8, remains very close to the perimeter of the guarding circular area of the frontal (closest to the sea) compromised well (w 1 ) for all simulation cases. The exact difference of x 1 À x t;1 , where Table 7 and verifies the observation. Figs 9 and 10, that follow, further explore the behavior of the ALOPEX stochastic optimization process. Fig 9 depicts the evolution of the objective function P(Q (k) ), with respect of k, for all simulation cases. Observe that, in all cases, ALOPEX drives expeditiously the objective function P(Q) to its, restricted by the constrains, maximum value and remains close to it within relatively small amplitude fluctuations. Fig 10 depicts the evolution of the well pumping rates Q ðkÞ i ðm 3 =dayÞ, with respect of k, for all simulation cases. Observe that, the balance between algorithm's deterministic terms (feedback and constrains) on one hand, and the stochastic term on the other, in order to achieve convergence within small fluctuations of the objective function around its maximum, is being achieved by several different motifs as it pertains to the behavior of the control variables Q i . For instance, the value of a control variable may converge within fluctuations at specific leveled value, as it happens for both variables of Case Ib in Fig  10. In cases, though, such as Cases Ia and IIa, in which local characteristics of the wells (x-coordinates and distance from the sea, minimum and maximum pumping rates) force their pumping rates to be close enough, it is possible to witness an "interweave" behavior during the optimization process. In Case Ia this happens only momentarily, as the algorithm is sensitive enough to realize even small differences (in the x-coordinates in this case) that affect the optimal value and "correct" its way to determine it.
Finally, as it pertains to the stability of the optimization results presented above, Table 8, that follows, contains the x-coordinates x s,i (see Fig 2) of the stagnation points (x s,i , y i ), i = 1, . . ., M, that correspond to the asymptotic optimal values (included in Table 5) of all simulation cases. As it can be verified by direct calculation, the difference x s,1 − x τ,1 , reported in the last  Table 8, satisfies Dx st x s;1 À x t;1 ¼ min hence all stagnation points remain between their associated critical toe interface points and the corresponding wells and, therefore, the optimization results satisfy the necessary stability requirement.
Observe, however, that as the magnitude of Δx sτ is relatively small, the optimal results presented in Table 5, although stable, are sensitive to small augmentative perturbations. For instance, increasing the optimal values of all control variables Q i reported for Cases IIIa and IIIb in Table 5 by 1.10% and 3.85%, respectively, would lead to unstable results. Sensitivity issues of the optimization results are further discussed in the section that follows.

Validity/Stability Issues of the Optimization Results
As we have already mentioned, results have shown that the sharp-interface approximation is reasonable in regional scale problems when the transition zone is narrow relative to the scale of the problem (cf. [9,15]). The validity study of the sharp-interface approximation, performed in [9], has concluded that the width of the transition zone decreases for: • higher values of the recharge parameter N, • lower values of the horizontal hydraulic conductivity parameter K, • lower distances x w of the wells from the coast, • higher values of the aquifer's height d.
In this section, in an attempt to create scenarios that enhance the validity of the optimization results presented in the previous section, we consider the model problems described by the Test Cases IIIa and IIIb and investigate validity, sensitivity and stability issues of the optimization results as it pertains to the recharge and hydraulic conductivity parameters N and K.
The optimization results obtained are summarized through Figs 11 and 12 and refer to the usage of optimization Profiles 1 and 2 respectively. Both figures explore the dependence of the optimal total pumping rate S(Q), and its associated stability related difference Δx sτ = x s,1 − x τ,1 , on the parameters N and K, while keeping the rest of the parameters fixed at their values included in Tables 1 and 2.
Inspecting Fig 11, that refers to the usage of the optimization Profile 1, one may observe that: • The optimal pumping rate S(Q), in both test cases IIIa and IIIb (Fig 11A and 11C respectively), increases linearly with respect of N and decreases with respect of K. After a relatively large value of N, S(Q) reaches the aquifer's maximum pumping capability Q A ¼ 10000m 3 =day activating the corresponding constraint.
• The stability related difference Δx sτ , as it pertains to the test case IIIb (Fig 11D), remains positive and therefore the necessary stability criterion is being satisfied for all values of N and K.
Observe, however, that as N increases and K decreases its magnitude is getting relatively close to zero. Therefore, the corresponding optimization results are sensitive to small augmentative perturbations and close to instability regions. As it pertains to the test case IIIa (Fig 11B) observe that, for large intervals of N and for K = 10, 25, the difference Δx sτ becomes negative hence the corresponding optimization results are unstable.
In summary, one may remark that, using optimization Profile 1, higher values of N and lower values of K decrease the width of the transition zone and increase the total pumping rate but at the risk of instability.
Utilization of optimization Profile 2, namely augmenting the set of constraints described in Eq (23) -and used in Profile 1-by the potential constraint Eq (22), may help to overcome the above instability problems. Indeed, inspecting Fig 12, that refers to the usage of the optimization Profile 2, one may observe that: • The optimal pumping rate S(Q), in both test cases IIIa and IIIb (Fig 12A and 12C respectively), increases with respect of N, while reaches a limiting value as K decreases. After a relatively large value of N, S(Q) reaches the aquifer's maximum pumping capability Q A ¼ 10000m 3 =day activating the corresponding constraint.
• The stability related difference Δx sτ , as it pertains to both test cases IIIa and IIIb (Fig 12B and  12D respectively), remains positive and therefore the necessary stability criterion is being satisfied for all values of N and K. Moreover, in all practical cases there exists a non-negative number N K ! 0 such that for N > N K the value of Δx sτ increases, as N increases and K decreases, improving stability. In summary, one may remark that, using optimization Profile 2, higher values of N and lower values of K decrease the width of the transition zone, increase the total pumping rate up to a limiting upper bound and improves stability. Therefore, utilization of optimization Profile 2 should be preferred over optimization Profile 1 as it may improve both validity and stability of the optimization results.
To conclude this section we point out that a complete sensitivity analysis is necessary, since the validity and stability of the optimization results depend on a number of other parameters too. This, however, exceeds the purpose of this work and will be presented elsewhere.

Conclusions
This work presents a new version of the ALOPEX stochastic optimization algorithm and studies its performance when applied on analytical models of saltwater intrusion of coastal aquifers. The objective of the optimization problem was to maximize the volume of the fresh groundwater pumped from the aquifer subject to a set of constraints that prevent the intrusion of the saltwater.
The mathematical formulation of the saltwater intrusion problem was based on the work presented by Mantoglou [14] using an analytical model that is based on the sharp interface approach and the use of the Ghyben-Herzberg relation.
During the course of our study, we: • Determined the amplitude of the algorithm's feedback and noise terms that accelerate the process considerably, • Presented a penalty system to implement the constraints, • Introduced a stopping criterion for termination of the iterative stochastic process.
The model case aquifer simulates the coastal aquifer at Vathi on the island of Kalymnos in Greece. The model aquifer's characteristics (Table 1) are identical to the ones used in [14]. Several test cases, consisting of two, three, five and eleven wells, were examined. Utilizing the proposed ALOPEX stochastic algorithm and for all test case problems we: • Derived the optimal pumping rates ( Table 5) of all active wells within a fixed large number of iterations (asymptotic case). The optimal results strongly compete the corresponding results from the literature while the CPU time for their computation varied from 1.97 sec to 17.36 sec.
• Tested the effectiveness of the stopping criteria. The corresponding optimal results obtained ( Table 6) found in complete agreement with the corresponding asymptotic ones and computed at a fraction of the CPU time (from 0.33 sec to 1.74 sec).
• Explored the convergence profile of the objective function to verify that in all cases the algorithm drives expeditiously the objective function to its maximum and remains close to it within small amplitude fluctuations (Fig 9).
• Computed the coordinates of the critical points (Table 7) as well as the stagnation points (Table 8) to verify the stability of the optimal optimization results.
Finally, we investigated the dependence of stability of the optimization results on the recharge and hydraulic conductivity parameters N and K, and showed that the usage of the potential constraint profile (Profile 2) improves the stability of the optimization results (Fig 12).