Optimal Well Locations using Genetic Algorithm for Tushki Project, Western Desert, Egypt

The well location problem is challenging due to the non-linear, discrete and often multi-modal objective function. The optimal distribution of productive well locations mitigates the groundwater productivity problem that threats the national projects in arid countries like Egypt. In this paper, a trial to mitigate this problem in Tushki National Project, south western desert, Egypt is supposed via delineating the optimal well locations and optimal pumping rates. The methodology of combination between simulation and optimization techniques is applied. The simulation outputs of groundwater flow system by Visual MODFLOW model is linked by the constructed Fortran Code for Optimal well Location using Genetic Algorithm (OLGA Code) for obtaining the optimum management of groundwater resources in this project. Two management cases are considered by running the model domain with adopted steady and transit calibrated parameters. The first management case deals with the present well locations and predicts the optimal value of the objective function (maximum pumping rate). In the second case, the optimal new well locations resulted from the OLGA Code is predicted from flexible well location with the moving well option. Also, the prediction of the future changes in both head and flow are studied in steady and transient states. Keywords— Groundwater, Optimization, Genetic Algorithm, Tushki project.


INTRODUCTION
Selecting optimal well sites is a valuable problem to solve; maximizing groundwater recovery increases oil reserves and minimizing costs improves profitability. One of the best methods of selecting optimal sites is the combined use of simulation and optimization (S/O) models. While simulation models basically provide solutions of the governing equations of groundwater flow, optimization models identify an optimal management and planning strategy from a set of feasible alternative strategies. The genetic algorithm (GA) is widely used to modify the parameters of groundwater flow models (Yan et al., 2003;Yao et al., 2003) and to solve the management models of groundwater resources (Mckinney and Lin, 1994;Liu et al., 2002;Zhu et al., 2003). Rana et al, (2008) conducted a study on a spatio-temporal optimization of agricultural drainage, using groundwater model and GA. Guan et al. (2008) proposed an improved genetic algorithm (IGA) to solve optimization problems with equality and inequality constraints. Hamid et al. (2009) presented a paper focusing on the S/O for conjunctive use of surface water and groundwater on a basin-wide scale, the Najafabad plain in west-central Iran. Gad et al., (2011), and Moharram et al., (2012) used an optimization model based on the combination of MODFLOW simulation with GA to maximize the total pumping rate from the Nubian sandstone aquifer in El-Dakhla and El-Farafra depressions, Egypt. Saafan et al., (2011) applied multi-objective genetic algorithm (MOGA) Code in El-Farafra depression, Egypt to study the maximum pumping rate and minimum operation cost. Gad and Khalaf, (2013) used MOGA model to develop the maximum pumping rate and minimum operation cost in the Miocene aquifer of Wadi El-Farigh, West Delta, Egypt.
In the other side, there are limited studies for the determination of optimal operating strategy, including unknown well locations and pumping rates, for groundwater systems to the best of our knowledge (Wang and Ahlfeld, 1994, Karahan and Ayvaz, 2005, and Ayvaz and Karahan, 2008). Saffi and Cheddadi (2007) developed an algebraic expression which gave the matrix of transient influence coefficients for one-dimensional semi-confined aquifer model and solved the governing equation by using a mixed compartment model. Tung and Chou (2004) integrated pattern classification and tabu search to optimize the average zonal groundwater pumping for an aquifer. As a general, the research conducted in this field dealt with the identification of locations and released histories of unknown groundwater pollution sources Datta, 2000, 2001;Aral et al., 2001;Ruperti, 2002;Singh et al., 2004;Sun et al., 2006), pumping well optimization for optimum remediation design (Wang and Ahlfeld, 1994;Huang and Mayer, 1997; Guan and Aral, 1999;Zheng and Wang, 1999;Mantoglou and Kourakos, 2007;Chang et al., 2007), and optimum well locations and pumping rates in the coastal sides (Park and Aral, 2004; Ferreira da Silva and Haie, 2007). Meira, et al. (2016) proposed technique for decision-making process in optimal well location including a step of risk analysis associated with the uncertainties present in the variables of the problem. The proposed technique was applied to two benchmark cases and the results, evaluated by experts in the field, indicated that the obtained solutions were richer than those identified by previously adopted manual approaches.
Due to the few of scientific researches related to optimal well locations especially in our national governmental and invest mental agricultural projects, a trial was carried out in this paper to through a light on the optimal location for new fourteen production wells in Tushki project required for irrigation of 25 thousand acres which will drill under the strategy of reclamation of 1.5 million acres project.

Site description
Tushki area, is a part of the Western Desert of Egypt; between latitudes 22° 14' 24.75" and 22° 50' 6.10" N and longitudes 31° 0' 44.06"and 31° 50' 4.05E, with total area of 2268 km 2 . Tushki project, one of the national governmental and invest mental agricultural projects, aims to reclaim 540,000 acres using both surface water from El-Sheikh Zaid artificial canal (1.5 km 3 /year) and groundwater extracted from about 316 productive wells for cultivating 36000 acres. Moreover, the agricultural development in the north and south of Abu-Simbel area, which depends only on the groundwater resource, will extract 22.5 million m 3 /year from the present 75 productive wells ( Figure 1). This operation policy is expected to have its side effect on groundwater regime of Nubia Sandstone Aquifer in Abu-Simbel Area (NSAASA) due to the absence of optimal well locations which will be the fulcrum of this research.

Geomorphological and geological settings
Referring to the geomorphological studies prepared by

Hydrogeological aspects
The NSAASA is mainly composed of three successive water-bearing units, namely; the Gilf, Abu Simbel and Lake Nasser. The middle unit is an aquiclude that extends under the lake and the lower unit overlying the granitic basement is a confined aquifer of fluvial sandstones ( Figure 3). The groundwater of NSAASA exists under unconfined to semi-confined conditions. The NSAASA maximum thickness reaches 400 m and decreases towards north and NW. The transmissivity of the NSAASA ranges between 389 and 1322 m 2 /day ( Table 2), and increases from SW to NE. The constructed equipotential contour maps ( Figure 4) show that the local recharge of the NSAASA depends upon the leakage from the Lake Nasser due west direction. The estimated groundwater flow rate based on these contour maps ( Figure 4) reaches 0.054 m/day near Nasser Lake and decreases to 0.044 m/day towards northwestern parts.

Field trip and data collection
The materials used in this paper are collected through carrying out four field trips in Abu Simbel area during the period 2010-12. The basic hydraulic parameters of NSAASA are estimated based on carrying out one step test and 13 long duration pumping tests (Tables 1 and 2). The different hydrologic data are obtained during the field trips. These materials also include collection of archival data (well drilling reports, WRRI 2002), registration of discharge, distribution of wells, proposed operating systems for both groundwater supply and reclaimed area beside recording depth to water for groundwater level changes.
The methodological approach used in this paper is based on the mathematical modeling techniques which combines between simulation and optimization model. A computer programming with FORTRAN language has been originally established to apply the principles of the genetic algorithm for studying the optimal location of the groundwater productive wells as a principal tool for groundwater resources management. The Optimal well Location applying Genetic Algorithm (OLGA Code) links MODFLOW with genetic algorithm technique to establish a simulation optimization groundwater model.

Simulation technique
Groundwater flow simulation is carried out applying MODFLOW software. It describes groundwater flow of constant density under non-equilibrium conditions in three-dimensional heterogeneous and anisotropic medium according to the following equation (Bear, 1979): (1) in which; K xx , K yy and K zz are values of hydraulic conductivity (L T -l ); along the x, y, and z coordinate axes; h is the potentiometric head (L); W is the volumetric flux per unit volume and represents sources and/or sinks of water (T -l ); S s is the specific storage of the porous material (L -1 ); and t is time (T).

Conceptual Model
The conceptual model of the NSAASA is based on the geology of the study area which is comprised of Nubian sandstone sediments. The hydrogeologic system is concerned two hydraulically connected through fault planes layers. The top hydrostratigraphic unit (Gilf Formation) is considered to be a sandy unconfined aquifer zone up to 95 meters maximum thick. The second confined aquifer zone of fluvial sandstones is separated from the sandy unconfined aquifer zone by aquiclude, and had maximum thickness of 79 meters. The scattered new communities around the asphaltic road in the N-S direction of the Abu Simbel area are primarily planned for agriculture.

Model domain
The model domain of NSAASA is selected to cover 3315 km 2 (51x65 Km- Figure 3). It is discretized using 130 rows × 102 columns rectangular cells. This discrimination produces 13260 cells in the model layer. The width of the cells along rows (in x-direction) is 500m and along columns (in y-direction) is 500 m. The grid geometry is shown in (Figure 5).

Boundary conditions
The hydraulic conditions at the boundaries of the model domain are specified ( Figure 5). Based on the groundwater flow pattern of the NSAASA ( Figure 4) and the field measurements of the groundwater head in the present productive wells (Table 2), the initial and boundary conditions of the local model is defined. The lower boundary (water body of Lake Nasser) is considered a time-constant head at level of 175 m asl (DRC 2012). The equipotential contour line of 130m asl characterizing to the upper boundary ( Figure 5) is far enough from the well fields to be genral head boundary condition. The left and right boundaries are assumed to be no flow boundaries due to the regional flow direction from SE to the NW.
The impermeable Basement rock is assumed to be the bottom boundary of the model domain with no flow and maximum depth 236.45m (Table 3 and Figure 6). The hydraulic heads and hydraulic conductance are assigned to boundary cells. These heads vary during the simulation process according to different stresses applied on the modeled area.

Aquifer properties
The aquifer material properties of the NSAASA include the transmissivity, the storage coefficient and the hydraulic conductivity. The transmissivity ranges from 389 and 1322 m 2 /day ( Table 2) and the storage coefficient from 1.2x10 -4 to 5x10 -4 (El-Sabri et. Al., 2010). Vertical heterogeneity is assumed adequate to allow treatment as multi layered aquifer. The vertical hydraulic conductivity, estimated from the aquifer thickness (

Model calibration
Calibration involves modification of the parameters, and forcing functions in the model in such a way as to assure that the model reproduced observed-head values (Hill, 1998(Hill, & 2000. The calibration process typically involves calibrating to steadystate and unsteady state conditions. The performed steady state calibration permits the adjustment of the hydraulic conductivity, where NSAASA storage changes are not significant. The relation between the calculated and observed heads is checked from the calculated-observed head curve and the variance appears as 45.966 % (Figure 8). This large value indicates great difference between the heads calculated by the model and actual measured heads. After many times of changing the k values, the variance between the observed and calculated heads is minimized to 3.61 % and the calculated head values is very closely related to the field-measured head values (Figure 7). Based on the calibrated model, the estimated hydraulic conductivity values are ranging from 4.33 to 7.60 m/day. These values are generally more slightly than the values calculated from the field data. The difference may relate to the benchmarks, uncertainty and grid refinement. The range of the resulted specific yield after the final calibration of the transient state is found to be varying from 0.14 to 0.34 (Fig. 8). It can be seen that, in general, there is a good agreement between the observed and simulated head values. The natural evolutional processes of reproduction, selection, crossover, and mutation are applied using probability rules to evolve the new and better generations. The probabilistic rules allow some less fit individuals to survive. The objective of this study is to maximize the benefit function Z with respect to pumping rate to obtain optimal location of wells, Q j as design variable, where: Max (2) in which:, N w is the total number of pumping wells and P j is penalty.
Moreover, the management objectives must be achieved within a set of constraints. The constraints may be decision variables. The objective function is subjected to the following constrains:

Pumping constraint
The pumping rates at potential pumping wells in the water demand are constrained to values between some minimum (Q j min ) and maximum (Q j max ) permissible pumping rates as the following: In the GA simulation, this constraint can be easily satisfied by restricting the population space of the design variables to be within the proposed limits. Hence, no special treatment is needed for this constraint.

Drawdown constraint
This constraint normally means to protect the ecosystem by avoiding excessive drawdown. In this work, the drawdown constraints were formulated to avoid mining and formulated as follows:- Accordingly, it postulates the relation; (4) In which: r w1 is drawdown at point 1 due to pumping well no 1, r w2 is drawdown at point 1 due to pumping well no 2, r w3 is drawdown at point 1 due to pumping well no 3, N c no of control points, rj is the drawdown at control point i caused by a pumping rate from pumping well j (value of rj in Table 4 to 6), di is the permissible drawdown at control point i equals to 30 m.

Water demand constraint
The NSAASA is considered as the sole source of water. This, therefore, means that the designed optimal pumping strategy must supply at least the minimum water demand. It is formulated as follows: (5) In which Q D is water demand.

Distance between wells constraint
This constraint obeys the relation: Where D i is actual distance between wells and D all is minimum distance between wells.

Location of wells
The locations of wells constraint is to be decided by the model itself within a user defined region of the model grid until the optimal location is reached (Figure 9).

Optimization Procedure of the Simulation-Optimization model
In groundwater management problems, there are two sets of variables, decision variables and state variables, where the decision variable is the pumping and injection rates of wells. Also other decision variables include well locations and the on/off status of a well. By optimization techniques the decision variables can be managed to identify the best combination of them. Moreover, hydraulic head is the state variable. The simulation model updates the state variables, and the optimization model determines the optimal values for all decision variables. In this work, MODFLOW FORTRAN code is used as the simulation of groundwater flow which links with GA optimization (OLGA Code). The flowchart for simulation-optimization model is given in (Fig. 10).

FIGURE 10: FLOWCHART FOR THE SIMULATION/OPTIMIZATION MODEL
It must refer that MODFLOW simulation without wells is used to generate the initial head, but initial population is generated from GA via random number for each well according to constraints and number of chromosome.

Checking the termination criterion
The optimization algorithm continues to be executed by iterating until the given termination criteria are satisfied (Fig. 10). Note that different termination criteria can be used to stop the computation. These may be: Stopping the computation after a given number of iterations; reaching a specific objective function value; no improvement in the objective function value for a specified number of passed iterations; or after limited time. In this study, the relative error definition is used to check the convergence of iterations. Thus, the convergence criterion is defined as, Where ε is a tolerance for the convergence of iterations, if two consecutive objective values (f1, f o ) satisfy the criterion.

Testing scenarios
The testing scenarios include three proposed pumping rates and well locations policies. The first scenario estimates the pumping rates from 75 productive wells with their present locations. The second scenario checks the present locations for these 75 productive wells and the optimality of their pumping rates. The third one proposes water exploitation policy aimed at increasing the present productive wells by 14 wells, as a result of the increase in reclamation activities (new 25 thousand acres), and delineates the optimal locations of these new wells and keeping the present locations of the other 75 present wells besides estimating the optimal pumping rates for all.

IV. RESULTS AND DISCUSSIONS
The expected interference between cones of depression of the present productive wells due to their configurations around the highway roads of the Abu Simbel area will affect the groundwater productivity of the NSAASA. Accordingly, the need to determine optimal well location strategy for the new drilled wells required for future expansion of the project (1.5 million acres national project) is essential.
In addition, applying the constructed OLGA Code for solving the optimization problems related to the optimal pumping rates and the optimal well locations, the initial pumping rates are assigned to each well as more than 700 m 3 /day and less than 1000 m 3 /day and then the optimization procedure explained before is used. The GA parameters used in the optimization problem are given in Table 4. The optimal well location and pumping rate are determined, based on OLGA Code, from one feasible solution to the other according to the calculated residual errors values where they improve quickly in the early iterations and change slowly in subsequent iterations till satisfying the constraints. On the other side, the results of the three testing optimal policies applying OLGA Code show the future predictions in the optimal pumping rates after 5, 10, 15, 20, 25, 30, 40 and 50 years. The first scenario investigates the proposed pumping rates from the 75 productive wells penetrating the NSAASA under their current locations. The OLGA simulation results show that the optimal pumping rates, under the present conditions of water demand Q D of 60000 m 3 /day, range from 57585 to 50144 m 3 /day with corresponding total drawdown ranges from 3.42 to 13.73 m during the simulation period of 50 years (Table 5).

is number of operation wells, r is maximum drawdown, Q min /well is min optimal pumping rate for well, Q max /well is max optimal pumping rate for well and Q opt total optimal pumping rate for all wells.
Moreover, by running the developed OLGA Code for time step of 5 years, the predicted head distribution maps of the NSAASA under the current pumping rates at years 2015, 2025, 2035 and 2060 are shown in (Figure 11). It is noticed from (Figure 11) that the groundwater levels change from 164.5 masl to 161.3 masl in the center of the cone of depression. Few productive wells in the southwestern part of the model domain show relatively slight drawdown amounts to 13 m at East El-Emaratiah Company well field. This is mainly attributed to the effect of increasing the thickness of water bearing formation towards the northwestern and western parts. In addition, one large cone of depression will develop in the cultivated areas.
In the same side, the results of the OLGA Code based on this scenario show that the predicted optimal location of the productive wells differ more or less from the present location ( Figure 12). Although the proposed OLGA approach successfully identified all the well locations and pumping rates of the present numerical study, identification of all the well locations and pumping rates may not be feasible for real-world applications since the number of pumping wells is usually very high. Instead, the proposed procedure may be applied to the suspected illegal pumping areas as noticed in the field trips.
As an extension of this study, the performance of the proposed OLGA Code should be tested on a real field system in a future study.

FIGURE 12: OPTIMAL LOCATION OF PRESENT WELLS (2ND SCENARIO)
Otherwise, the predicted head distribution maps of the NSAASA under the 2 nd scenario ( Figure 13) show groundwater levels change from 165.5 masl to 162.5 masl in the center of the cone of depression. It is noticed from (Figure 13) that the normal decline in the groundwater levels (3 m) may attribute to the small thickness and low hydraulic conductivity of the NSAASA in these localities. Moreover, the cone of depression will appear in the cultivated areas during the simulation time (with diameter of 2.5 to 4 km approximately). While the southern and western cultivated parts of the model domain does not be affected with this scenario. This may attribute to the effect of the great aquifer thickness and the presence of thin clay layers in the southern and western parts rather than the geologic structure impact. Accordingly, this reflects the low potentiality of the northern cultivated parts for positive response to this scenario.

FIGURE 13: PREDICTED HEAD DISTRIBUTION MAPS OF THE NSAASA FOR OPTIMAL WELL LOCATIONS APPLYING 2ND SCENARIO A) AT 2015, B) AT 2025, C) AT 2035, AND D) AT 2060
The third scenario proposes water exploitation policy aimed at increasing the present productive wells by 14 wells, as a result of the increase in reclamation activities, under the constraints of 500 m as minimum distance between wells and 900000 m 3 /day as water demand Q D , and delineates the optimal locations of these new productive wells and keeping the present locations of the other 75 present productive wells, besides estimating the optimal pumping rates for all. The GA parameters used in this optimization case are given in (Table 7). The OLGA Code runs for time step of 5 years. The predicted optimal location of the new productive wells is checked ( Figure 14) and the predicted head distribution maps of the NSAASA for the optimal pumping rates at years 2015, 2025, 2035 and 2060 are shown in (Figure 15).

FIG. 14: OPTIMAL LOCATION OF THE NEW PRODUCTIVE WELLS (BLUE COLORED) (3RD SCENARIO)
It is noticed from ( Table 7) that the maximum drawdown r ranges from 4.27 to 14.56 m while Q min /well begins with 923.2 and ends by 800 m 3 /day in case of number of total operating wells N of 91 productive wells. Moreover, Q min /well begins with 923.2 and ends by 800 m 3 /day during the simulation period. The total optimal pumping rate for all wells (Q opt ) ranges between 87187 and 79747 m 3 /day.
However, it is noticed from the predicted head distribution maps ( Figure 15) that the groundwater level decreases from 163m asl in the beginning of the simulation period (the year 2015) to 158m asl in the end of the simulation period (the year 2060). The depression cones related to the first and the second scenarios are mitigated in the northern part of the model domain.
Moreover, limits are placed on the total optimal pumping rates of the productive groundwater wells beside its optimal locations in the reclaimed areas. Accordingly, it is preferred to extend the future reclamation activities parallel to the western side of the Khour Tushki area. Although the results of this scenario provide more or less good degree of confidence in the optimal location of the new productive wells and optimal pumping rates, the expected sharp decline in groundwater levels in the concerned area assumes more applied studies. The Q opt in the 3 rd scenario is actually lower than that in the 2 nd scenario because the second scenario studies the optimal locations for these 75 productive wells and their optimal pumping rates. The third one proposes water exploitation policy aimed at increasing the present productive wells by 14 wells. The predicted value of (r) based on the 3 rd scenario is more or less similar to the results of the other two scenarios although the number of the operating wells is increased by 19%. This reflects the great importance to apply the optimal well location concept in any new reclamation projects. As a result of the above discussion of the OLGA Code optimal well location and pumping rates results, it can be concluded that the groundwater withdrawal from the NSAASA under this optimization study could be safely conducted.

V. CONCLUSION AND RECOMMENDATION
This paper is devoted to propose an assisted methodology to identify representative models for optimal well locations. To do so, first a mathematical equations were developed to model the representativeness of a subset of constrains with respect to the full set that characterizes the well location problem. Then, an optimization tool was implemented to identify the representative constrains of any problem, considering not only the main output variables, but also the probability distribution of the attribute-levels of the problem. Accordingly, OLGA Code was linked MODFLOW with genetic algorithm technique to establish a simulation optimization groundwater model. This model was applied for NSAASA to develop the optimal location of the groundwater productive wells with maximum pumping rate. Three scenarios were tested to choose the proper water exploitation policy. For the first scenario, the predicted groundwater drawdown after simulation period of 50 years ranged from 3.42 to 13.73 m while the corresponding optimal pumping rates ranged from 57585 to 50144 m 3 /day. Under the second proposed scenario, the predicted drawdown (r) and the corresponding total optimal pumping rate for all wells (Q opt ) ranged from 2.22 to 12.13 m and between 102038 and 100773 m 3 /day, respectively. Moreover, the predicted value of (r) based on the 3 rd scenario was more or less similar to the results of the other two scenarios although the number of the operating wells was increased by 19%. Computational restrictions were common and likely the cause of the systemic lack of algorithm benchmarking for well location optimization. Still, some valuable insight can be gained from these results that could guide further investigation.
Based on the results of the OLGA model, it is highly recommended to choose the new productive well locations in staggered system parallel to the western side of the Khour Tushki area. More applied studies are needed for verification the results of this optimal well location model. Also, this study should be extended to contain more test cases and more comprehensive selection of algorithms including interpolation methods, ad joint methods and hybrid methods.