Pumping optimization of coastal aquifers using probabilistic search – case study: Quaternary aquifer of El-Arish Rafah, Egypt

Groundwater contamination through saltwater intrusion threatens the health of many people living in coastal areas. Coastal aquifers in arid and semi-arid zones which border the sea are an important source for water. Coastal areas are usually heavily urbanized zones, a fact that makes the need for freshwater even more acute. To meet growing demands, increased freshwater is extracted from these aquifers. Consequently, wells used for freshwater may be contaminated by saltwater intrusion.

In recent decades, stochastic optimization combined with flow simulation has been widely used to control saltwater intrusion and the general management of coastal groundwater. In most previous coastal groundwater management studies, simulated flow models are based on numerical methods such as the Finite Difference Method (FDM) (Bakker 2003) and the Finite Element Method (FEM) (El-Alfy et al. 2015) or analytical methods such as the theory of images (Abdel-Gawad & Shamaa 2004). The use of numerical or analytical methods for coastal aquifer simulation along with optimization algorithms is constrained by the huge computational burden involved. These approaches differ in the used methods to simulate the coastal groundwater flow system, the type of groundwater management problems and the algorithms used to solve these problems. Different objective functions and constraints have been employed to control this problem such as total pumping rate maximization, drawdown minimization, pumped water salinity minimization and pumping cost (Javadi et al. 2015; Christelis & Mantoglou 2018; Yang et al. 2018). In simulation–optimization modeling, the objective function is evaluated using an optimization algorithm, while the simulation model is utilized to satisfy the constraints.

Studies of the application of stochastic optimization algorithms for the management of coastal aquifers are reported in the literature. Cheng et al. (2000) reported the first application of the structured messy genetic algorithm (GA) combined with analytical solutions. GA and its variants are the most popular of stochastic optimization methods for such problems (Abd-Elhamid & Javadi 2011; Javadi et al. 2015). Examples of other applications include particle swarm optimization (Karatzas & Dokou 2015), ant colony optimization (Ataie-Ashtiani & Ketabchi 2011), differential evolution (Elçi & Ayvaz 2014), simulated annealing (Rao et al. 2004), and simplex simulated annealing (Kourakos & Mantoglou 2009). The literature review shows that some optimization algorithms were rarely used in groundwater resources management and, to the authors' knowledge, there is no report of the application of the Shuffled Frog Leaping Algorithm (SFLA) for coastal groundwater management.

Decision variables are related to pumping rate schemes (Cheng et al. 2000; Ataie-Ashtiani & Ketabchi 2011), water table level or drawdown (Katsifarakis & Tselepidou 2009), well locations (Elçi & Ayvaz 2014), the salinity of pumped water from wells (Ataie-Ashtiani et al. 2014), saltwater volume into the aquifer (Emch & Yeh 1998), operating costs (Abd-Elhamid & Javadi 2011), net benefits (Ferreira da Silva & Haie 2007), recharge rate schemes (Kourakos & Mantoglou 2011), trade-off between environmental and social issues, and interactions between surface and subsurface resources. Constraints include pumping limits (Cheng et al. 2000; Rao et al. 2004; Katsifarakis & Petala 2006), saltwater toe location (Cheng et al. 2000; Ferreira da Silva & Haie 2007; Ataie-Ashtiani & Ketabchi 2011), salt concentration (Kourakos & Mantoglou 2011; Ataie-Ashtiani et al. 2014), water table level (Katsifarakis & Tselepidou 2009), strategic nature of the aquifer, existing infrastructures, and historical rights rates (Abarca et al. 2006).

Over the years, a large number of models have been developed and used to study saltwater intrusion in coastal aquifers. FDM and FEM are widely used to simulate the groundwater flow. These models are SUTRA (Ojeda et al. 2004), SEAWAT (Kopsiaftis et al. 2017, 2019), CODESA-3D (Paniconi et al. 2001), SWIFT (Ma Tain et al. 1997), MOCDENS3D (Oude Essink 2001), SHARP (Barreto et al. 2001), Quasi-three-dimensional finite element (Zhou et al. 2003), and SHEMAT (Scholze et al. 2002).

Probabilistic Global Search Lausanne (PGSL) is a global search algorithm proposed by Raphael & Smith (2003a). It depends on selective sampling of the search space using a probability density function (PDF) that is varied dynamically during the search process. PGSL starts with a population of random solutions and ranks them based on their objective function values. Regions that contain good solutions, having higher probabilities, are searched with greater intensity. The application of PGSL on many benchmark problems with multi-variable non-linear objective functions demonstrated that PGSL performs better than GA, and simulated annealing (Raphael & Smith 2003b). There are a few applications of PGSL for water resources management such as minimum cost design of irrigation canals (Adarsh 2012; Adarsh & Sahana 2013), however it is not used in the literature for coastal groundwater management.

The objective of the current study is to calculate the optimal pumping rates in an unconfined coastal aquifer under steady-state conditions and assuming a sharp saltwater interface. Simulation-optimization models are developed using both the PGSL and SFLA optimization algorithms combined with FEM on Strack's (1976) formulation to simulate the groundwater flow. The results of the two optimization algorithms are compared. Both methods are then verified through their application on the aquifer system under the City of Miami Beach in northern Spain. Finally, they are applied to the Quaternary aquifer of El-Arish Rafah, Egypt.

In this study, for the unconfined aquifer, the following assumptions are considered (Figure 1): (1) sharp interface between saltwater and freshwater, (2) the three-dimensional geometry flow equation is reduced to two-dimensional by Dupuit's assumption, (3) steady-state flow condition, (4) the Ghyben–Herzberg assumption is utilized to interpret the interface location, (5) wells fully penetrate the aquifer thickness, (6) horizontal impervious aquifer bed, and (7) wells are only within the fresh water zone.

The seepage faces around different wells are ignored within the present study.

PGSL is a search algorithm based on sampling of the search space, without using special operators, using a probability distribution function (PDF) defined over the entire search space (Raphael & Smith 2003b). Potential solutions are generated randomly and evaluated using a defined objective. The basic assumption is that good solutions are more likely to be found in the regions with higher probabilities. Hence, the search is concentrated in regions containing good solutions. The evolution of solutions is carried out through updating the PDF in regions where good solutions are found. Unlike other techniques, there is no point-to-point movement towards a global optimum solution. Instead PGSL relies on a sampling procedure similar to the Monte-Carlo technique. In the beginning, each variable range (maximum value – minimum value) is divided into a number of intervals and a uniform probability is assumed within each interval. As the search progresses, intervals containing good solutions are assigned higher probability values and the search space along the variable is narrowed. As such, the search space is gradually narrowed down until convergence is achieved. The algorithm includes four nested cycles named Sampling, Probability updating, Focusing, and Sub-domain. The PGSL involves the following steps:

• 1.

  Initialization: The PDF of each decision variable is assumed having a uniform distribution over the entire domain.

• 2.

  Sampling cycle: Solutions are generated randomly using the current PDF. Then, each solution is evaluated and the best solution is selected.

• 3.

  Probability updating cycle: In this cycle, each decision variable PDF is modified using the probability-updating algorithm. This ensures that the probabilities of regions containing good solutions are increased.

• 4.

  Focusing cycle: After each iteration, the search is then focused on the interval containing the best point. For each decision variable, the interval containing the best point is then subdivided into smaller intervals.

• 5.

  Sub-domain cycle: At the beginning the entire space is searched, while in subsequent iterations a sub-domain of the problem is selected for search. The sub-domain size gradually decreases by changing the minimum and maximum values of each variable until the solution converges to a point.

The main parameters of the PGSL algorithm are: the number of evaluations of the objective function (numeval), number of focusing cycles (NFC) equals 20 × number of variables (Raphael & Smith 2003b), number of probability updating cycles (NPUC), Number of sampling cycles (NSC), and number of sub-domain cycles (NSDC) equals integer of numeval/(NSC × NPUC × NFC) (Raphael & Smith 2003b). The flowchart describing the PGSL algorithm is shown in Figure 2.

SFLA is used to seek a global optimal solution of combinatorial optimization problems. It combines the benefits of the idea of mixing information from parallel local searches and the social behavior-based particle swarm optimization algorithms to move towards a global solution. In this algorithm, the generation of a number of memeplexes is carried out, each one containing a number of random frogs (solutions). The frogs with the best fitness (X_B) and the worst fitness (X_W) in each memplex are determined. Also, in each cycle of evolution, the frog having the global best fitness (X_G) is determined and the worst fitness frog is updated. The SFLA involves the following steps:

1. Random initialization of initial population of possible solutions.

2. The fitness of each solution is calculated based on the given objective function and constraints violation penalties.

3. Determination of the global best frog position (X_G).

4. Population is partitioned into subsets named memeplexes representing different cultures of frogs. Then, the following steps are carried out for each memeplex:

   • Determine both the best and worst frog positions (i.e. X_B and X_W).

   • Update the worst frog position using Equations (13) and (14) (Elbeltagi et al. 2007):

     

     (13)

     

     (14)

     in which D(l) is a change in the position of the worst frog, X_W(l), located in memeplex l, X_B(l) is the position of the best frog located in memeplex l, X_{W new}(l) is the worst frog new position located in memeplex l, X_{W current}(l) is the worst frog current position located in memeplex l; c is a search acceleration factor proposed by Elbeltagi et al. (2007), r is a random number ranging from 0 to 1, n_m is the number of memeplexes and D_max is the maximum allowed change in the worst frog position.

   • If the obtained new position is better than the worst, change the worst frog with the new one. Otherwise, re-apply both Equations (13) and (14) replacing XB (l) by XG.

   • If the obtained new position is better than the worst, change the worst frog with the new one. Otherwise, generate a new frog randomly.

5. Perform step (4) for a specified number of trials until the position of the worst frog is improved.

6. Shuffling process is carried out for all frog positions in all memeplexes.

7. Steps 2–6 are repeated until no improvement occurs in the global best frog position (best solution).

The main parameters of the SFLA are: number of frogs, number of memeplexes, trials before shuffling and number of iterations after shuffling, maximum allowed change in frog position D_max, and the search acceleration factor c value. The flowchart describing the steps of SFLA is presented in Figure 3.

The PGSL and SFLA optimization algorithms are verified through applying them on a well-known benchmark problem. At the same time, the results of both models are compared. The two optimization models used the FEM to simulate the flow in coastal aquifers. The two algorithms are applied on the unconfined coastal aquifer at Miami Beach in the north of Spain (Benhachmi et al. 2003) to maximize both the net benefit on one scenario and maximize the total pumping on another scenario. In the two scenarios, pumping rates from a pre-existing system of wells are considered as the decision variables. The studied unconfined aquifer is considered among the most important sources of water supply for the city of Miami Beach. The Miami coastal aquifer consists of single-layered unconsolidated sediment of Quaternary age with pore permeability corresponding to coastal piedmonts and alluvial fans. The sediments consist of clay and gravel, and overlay blue clay of Pliocene age, constituting the effective lower impervious hydrologic boundary. The hydraulic parameters of the studied aquifer are its length along the sea coast (Y_D = 5,000 m) and perpendicular to the sea cost equals 4,500 m, hydraulic conductivity (K = 14 m/d), uniform discharge per unit width of the aquifer (q_u = 1.2 m³/d/m), depth of impervious bed from M.S.L. (d = 30 m), saltwater density (ρ_s = 1,025 kg/m³), and freshwater density (ρ_f = 1,000 kg/m³) (Benhachmi et al. 2003). Table 1 shows the coordinates of the aquifer pre-existing wells, their allowable maximum and minimum pumping rates, and their corresponding ground level.

The management strategy is settled to maximize the economic benefits from the pumped water and to minimize the cost corresponding to lifting the water, Equation (1). The values B and C shown in Equation (1) equal 0.01 $/m³ and 0.0002 $/m³/m respectively (Benhachmi et al. 2003). In Miami Beach there is a region near the sea coast with confined conditions and a groundwater level higher than the ground level (Abdel-Gawad & Shamaa 2004). Consequently, in cases where the water level is higher than the ground level, the value of the constant C is set equal to zero. Several trial runs are carried out to determine both the suitable parameters' values of the SFLA and the PGSL (Table 2). The largest values obtained for the net benefit (Equation (1)) are 24.93 and 23.07 $/d for both the SFLA and PGSL with corresponding total pumping rates equal to 3,849.4 and 3,999.8 m³/d, respectively. The obtained values of pumping rates represent 64.2 and 66.7% for the SFLA and PGSL of the total water entering the aquifer, respectively. The corresponding total numbers of the objective function evaluations, Equation (1), are 19,946 and 2,000 for the SFLA and the PGSL, respectively. Table 3 compares the discharge of each well corresponding to the largest obtained net benefits of the two models in the current study and other previous studies. It is noticed that most of the water is withdrawn from wells having low ground level to decrease the lifting cost of pumped water. As reported in Table 3, the SFLA gives a higher value for the benefit in comparison with the two previous studies while the PGSL gives a close value for the higher reported value. Also, it is found that the economic benefits obtained from the SFLA are larger than those obtained from the PGSL which has a higher total pumping rate.

In this scenario, the unconfined aquifer under Miami Beach is restudied to maximize the total pumping of freshwater, Equation (2). After performing several trial runs, the suitable values of SFLA and PGSL parameters are listed in Table 2. The largest total pumping rates obtained are 5,427.8 and 5,310.2 m³/d, respectively, for both the SFLA and PGSL that represent about 90.5 and 88.5% of the total water entering the aquifer. The corresponding total numbers of evaluations of the objective function, Equation (2), are 64,079 and 2,000 for the SFLA and PGSL, respectively. Table 4 presents a comparison between the discharge of each well corresponding to the largest total pumping of the two models in the current study and previous studies. As presented in the table, Abdel-Gawad & Shamaa (2004) obtained the largest total pumping rates. However, when checking the well pumping rates given by Abdel-Gawad & Shamaa (2004) against both constraints shown in Equations (8) and (9) using the developed FEM simulation method, it was found that seawater had intruded into several wells (Table 4). This may be due to the different method they adopted (the theory of images) to simulate the hydraulic response in the aquifer or due to the incorrect number of mirrors used to represent the theory of images. However, the two developed models also presented a good solution for the total pumping with the SFLA outperforming the PGSL.

As such, both the PGSL and SFLA optimization algorithms are able to reach greater solutions without violating the constraints when compared to previous studies.

The methodology adopted in the comparison of the two optimization algorithms' performance begins by performing 20 consecutive trial runs corresponding to each one. Then, a comparison is carried out between both based on a set of criteria, namely: (1) the success rate, S_r [=(N_s/N_t) × 100] (where N_s is the number of trials each model succeeds, to catch a value for the objective function greater than the maximum corresponding value reported in the literature and N_t is the total number of trials); (2) the minimum, maximum, mean, and standard deviation of the obtained solutions; and (3) mean numeval. The maximum values of the objective functions reported in the literature are 23.81 $/d and 5,116.1 m³/d respectively for scenarios 1 and 2 without constraints violations (feasible solution). It is worth noting that for fair comparison, the run time for each trial for both algorithms is kept constant (i.e. 45 minutes); this time was selected based on the recorded average run time during the previous experiments. The evaluation criteria values associated with both scenario 1 (i.e. maximization of economic benefit) and scenario 2 (i.e. maximization of total pumping rates) for the two developed models are summarized and presented in Table 5. From this table it is seen that the SFLA has the highest success rate, for the two scenarios of a value equal to 90% for scenario 1. This means that the SFLA is able to catch an optimal solution for the economic benefit greater than the corresponding maximum value reported in the literature for 18 runs out of the 20 carried runs. Also, the SFLA has higher minimum, maximum, and mean values for the two scenarios. The smaller standard deviation of the objective function is corresponding to the PGSL in both scenarios, with a value of 0.52 for scenario 1 and 116.3 for scenario 2. Consequently, the performance of the SFLA is better than the PGSL in terms of the success rate, minimum, maximum, and mean values for the two scenarios. The PGSL performed better in terms of standard deviation of the objective function for both scenarios. The mean numeval for the PGSL was significantly smaller than the SFLA in the two scenarios. In scenario 2, for example, the mean numeval for the PGSL was 3,753 while it was 23,654 for the SFLA. Thus, the PGSL would be very beneficial, especially when the number of variables increases, as the processing time increases exponentially.

Although the SFLA performs better than the PGSL, the latter outperforms other models reported in the literature.

The Quaternary aquifer in north Sinai consists of three areas, namely: the area from Rafah to east of El-Arish, the Delta El-Arish valley and the area from Bir El-Abd to west Rommana (Figure 4). The three areas have different hydrogeological conditions and depositional environments. The coastal aquifer at El-Arish Rafah (Figure 4) is considered to apply the proposed two optimization models with the purpose of maximizing the total pumping rates from a pre-existing system of wells similar to scenario 2 of the model verification.

The lithology of El-Arish Rafah coastal aquifer is grouped into three main units, namely, sand/gravel, clay and Kurkar (i.e. calcareous sandstone). Figure 5 shows a schematic sectional view for the conceptual model of Quaternary aquifer in El-Arish Rafah area. Based on the site characterization, the upper part of the aquifer is sand/gravel (sand dunes, old beach deposits, and alluvial), and the lower part consists of Pleistocene Kurkar formations. In many boreholes the clay layer between the two formations is absent or exists in lenses in which it works as a semi-confining layer for the Kurkar formations (Abdallah 2006). In this study, the coastal aquifer at El-Arish Rafah is considered as sallow and phreatic. Vertical homogeneity is assumed so the three layers have been merged to constitute a single equivalent layer with an equivalent thickness and hydraulic conductivity. Rainfall is the main source of groundwater recharge in the studied area. All hydraulic parameters of the studied aquifer are obtained from the central laboratory of the Desert Research Center, Egypt, and any missing data are reasonably assumed. The length of the aquifer along the shore line (Y_D) is 43,200 m and perpendicular to the shore line equals 8,150 m, the equivalent hydraulic conductivity (K) is 33 m/d, the depth of the impervious bed from the M.S.L. (d) is 40 m, saltwater density (ρ_s) is equal to 1,025 kg/m³, and freshwater density (ρ_f) is 1,000 kg/m³. Both the available point measurements of the flow parameters and hydrological measurements collected by the central laboratory of the Desert Research Center, Egypt, are adopted to calibrate the aquifer recharging flow (El-Sayed 2014). Recharge zones are selected based on the analysis of the colors of a satellite image (Figure 6). Such color composition shows vegetated areas where expected returned irrigated water can recharge the aquifer. Calibrated recharge values reflect recharge due to rainfall, returned water from different water uses and possible unknown recharge from other underlying formations are not considered. As shown in Figure 6, the recharge has a maximum value of 3 × 10^–3 m/d in the north-east of the study area, while it equals 1.2 × 10^–7 m/d in the southern boundary (Abdallah 2006). Table 6 shows the coordinates of the different pre-existing wells in the aquifer, their allowable maximum and minimum pumping rates, and their corresponding ground level.

In this application, the parameters setting for both the SFLA and the PGSL are summarized in Table 2. After performing several trial runs for each algorithm, the largest total pumping rates obtained were 42,861.5 and 46,387.72 m³/d, respectively, for both the SFLA and PGSL which represents about 22.1 and 23.9% from the total water entering the aquifer from rainfall. This is due to the non-uniform distribution of wells over the studied domain (the existing wells are not distributed over the whole studied area) and also the number of working wells located in low rainfall recharge zones (Figure 6). The corresponding total numbers of evaluations of the objective function, Equation (1), are 113,133 and 2,000 for the SFLA and PGSL, respectively. From the results obtained, the PGSL outperforms the SFLA and, consequently, the PGSL is more efficient for large scale problems, where the best results are obtained in a smaller runtime of about 35 minutes per one run. Table 6 lists a comparison between well discharges, in each well system, corresponding to the largest total pumping of the two algorithms. Using the SFLA optimization algorithm, the total pumping is obtained from only nine wells out of the 31 wells while 19 worked wells were determined by using the PGSL optimization algorithm.

Figure 7 shows the maps of the water levels contours (measured from the impervious bed) corresponding to each well system listed in Table 6. The locations of the wells have been plotted on the maps. The symbol (●) represents a working well and the symbol (+) represents a non-working well. The shaded areas represent the extent inland of the saline wedge and the contour line corresponding to at 41 m which represents the sharp interface toe. From the studied aquifer, both the SFLA and PGSL optimization algorithms can be used successfully for management of coastal aquifers with saltwater intrusion.

In this study, a probabilistic global search optimization algorithm, namely PGSL, is used along with a hydraulic simulation module to optimize the pumping rates from an unconfined coastal aquifer. The results of the PGSL are compared with the well-known Shuffled Frog Leaping Optimization Algorithm (SFLA). The FEM is applied on the linear formulation of Strack (1976) to simulate the hydraulic response of the steady state homogenous aquifer. A sharp interface between saltwater and freshwater is assumed. The two optimization algorithms are verified through applying them on the unconfined aquifer beneath Miami Beach to obtain the maximum pumping rate in one scenario and the maximum economic benefit in another scenario and comparisons are carried out with previous studies. Then, statistical analysis is carried out to compare the performance of each optimization algorithm. Both algorithms are then applied on the Quaternary aquifer of El-Arish Rafah, Egypt. The results obtained indicate that the performance of the SFLA is better than the PGSL in terms of the best results obtained. However, the PGSL, which is a relatively new optimization technique, outperformed other traditional stochastic optimization techniques (e.g. the GA). Also, the run time for the PGSL was significantly smaller than the SFLA. The two models can effectively and efficiently be used to solve real groundwater management problems. From the studied aquifers, it is noted that the PGSL does not give the optimum solution but the obtained results are close to the optimum. As such, the PGSL would be a suitable tool in solving large scale problems (i.e. real-life aquifers) to give an initial orientation, not a final decision.