Coupled Sharp-interface and Density-dependent Model for Simultaneous Optimization of Production Well Locations and Pumping in Coastal Aquifer

The main management challenge in coastal aquifers is to prevent saltwater intrusion, ensuring ample freshwater supply. Saltwater intrusion happens due to unregulated pumping from production wells. Therefore, it is essential to have an effective management policy, which ensures the requisite amount of freshwater to be withdrawn from coastal aquifers without causing saltwater intrusion. A methodology for optimizing production well locations and maximizing pumping from production wells is presented to achieve these conflicting objectives. The location of production wells directly affects the amount of freshwater pumped out of the coastal aquifer. Simultaneous optimization of production well locations and pumping from the same is achieved by linking mathematical simulation models with the optimization algorithm. A new methodology using coupled sharp-interface and density-dependent simulation models is developed to find optimal well locations and optimize the amount of freshwater pumped from the coastal aquifer. The performance of the developed methodology is evaluated for saltwater intrusion in the coastal city of Puri, India. The performance evaluation results show the developed methodology's applicability for managing saltwater intrusion while maximizing freshwater pumping in coastal aquifers under constraints of well location.


Introduction
Coastal regions have been the seat of large, economically viable settlements with vibrant populations along its shorelines throughout history. Large amounts of freshwater are required for sustaining this population, which is met from pumping of groundwater from coastal aquifers. The reliance on groundwater grows due to increased uncertainties around reliable quantity of quality surface water. These coastal groundwater aquifers are regularly overexploited to meet increased demand. The water table falls due to overexploitation of groundwater resources, disrupting the delicate balance between freshwater and saltwater. Falling groundwater table results in the reversal of groundwater flow, allowing saline seawater to occupy the space of freshwater causing saltwater intrusion (Frind 1982).
The application of Linked Optimization Simulation (LOS) in pumping optimization with fixed well locations has been documented by various authors (Mantoglou and Papantoniou 2008;Park and Aral 2004). The problem is frequently treated as a constrained optimization problem to optimize freshwater pumping while keeping pumping rates within prescribed limits, and salt concentration around well sites below a certain level (Abd-Elaty et al. 2021;Sreekanth and Datta 2014), or toe of the interface (Mantoglou and Papantoniou 2008) does not cross the well locations. The fundamental advantage of a simulation-based optimization method is that the effect of any pumping strategy on a saltwater intrusion can be realized in advance. Various optimization algorithms ranging from linear programming (Culver and Shoemaker 1993) to heuristic optimization techniques such as genetic algorithm (Abd-Elhamid and Javadi 2011; Cheng et al. 2000), quadratic programming (Willis and Finney 1988), non-shorting genetic algorithm-II (Dhar and Datta 2009), evolutionary annealing simplex (Kourakos and Mantoglou 2009) have been used in coastal aquifers for pumping optimization from existing well locations.
Thus, location of the production wells has a significant impact on the amount of freshwater that can be pumped from the coastal aquifer. Thus, it is imperative that the locations of the production wells be simultaneously optimized along with the total pumping for sustainable coastal aquifer management. Real field implementation of production wells involves large number of constraints to be satisfied in designing new pumping locations (Guan and Aral 1999). For example, what should be the minimal distance between two production wells? What should be the limiting position of the saltwater-freshwater front from a production well? Where should the well be located that it does not capture any pollutant from a poor water quality zone, if any? What is the safe/maximum pumping rate from the designed wells?
The above objective can be achieved by simultaneous optimization of production well locations and pumping rates (Dey and Prakash 2019b;Todd 1959). This study presents a LOS based methodology for simultaneously optimizing the production well locations and the pumping rates from the production wells so as to maximize the freshwater pumping without causing saltwater intrusion. In the developed methodology the amount of pumping is adjusted for a corresponding set of well locations, adhering to the constraints mentioned above to determine the maximum freshwater pumping without causing saltwater intrusion. Both, production well locations and pumping are treated as explicit decision variables, which increases the complexity of the problem. Not much work has been reported wherein production well locations and pumping from these production wells have been simultaneously optimized for managing saltwater intrusion in a coastal aquifer problem. This study tries to address this specific gap.
A new method based on the sharp-interface and density-dependent models is developed (Christelis and Mantoglou 2016) to simultaneously optimize production well locations and pumping rates in a coastal aquifer. Pumping rates from the production well and the location of the production wells are considered as explicit decision variables to be optimally determined by the optimization model. Various factors used as constrained in designing optimal production well locations are distance between the production wells, location of the wells from the aquifer boundary, safe distance from poor water quality zone, and intrusion of saltwater interface (Ahmadi et al. 2021). The optimization algorithm essentially tries to navigate efficiently through the space of feasible solutions to find the most optimal solution. The developed methodology is first validated against two hypothetical scenarios and then applied as a proof of concept for the coastal aquifer of Puri City in India.

Methodology
The main objective of the developed methodology is to maximize the freshwater pumping from optimally located production wells in a coastal aquifer to meet the freshwater demand without causing saltwater intrusion. The optimal values of the pumping rates and the optimal location of the wells are simultaneously estimated by the developed LOS model.
In this method, first, sharp-interface model based on Strack (1976) potential is linked with Simulated Annealing (SA) (Goffe 1996) to find optimal pumping for varying well locations. Subsequently, density-dependent simulation model (SEAWAT) (Guo and Langevin 2002) is linked with Particle Swarm Optimization (PSO) (Shi and Eberhart 1998) to adjust the outcomes of the sharp-interface model (Dey and Prakash 2020). This iterative process is continued until the total pumping of freshwater is maximized without causing saltwater intrusion at proposed production well locations optimally chosen by the methodology. Coupling of density-dependent model to sharp-interface model is described in Sect. 2.2.

Formulation of Pumping Optimization Problem
The main objective is to maximize freshwater pumping for a given time period without causing saltwater intrusion at the designed pumping well locations. Thus, pumping well locations become an explicit decision variable in the formulation. The objective is mathematically expressed using the following equations: where, Q t is the total freshwater pumped from the coastal aquifer, Q i is the pumping from i th well, d i is the distance of the well from the coast, d i is the distance of the toe of sharpinterface from the coast, N is the total number of pumping wells, x i , y i is the location of (1) (2) Subjected to ∶ Q Min the i th well, x τi location of the interface at the x th row corresponding to the i th well, Q Min i and Q Max i is the minimum and maximum allowable pumping from i th well. D is the minimum safe distance between toe of interface and well location.
The freshwater pumping from the production well is maximized satisfying the constraints specified by Eqs. (2) and (3). Equation (2) limits the the maximum and minimum allowable pumping from any production well such that the salt water never reaches the production well (Eq. (3)), thus preventing saltwater intrusion. Locations of the pumping wells x i , y i is an explicit decision variable, which is simultaneously optimized to maximize freshwater pumping. While deciding the optimal well locations, the simulation model needs to satisfy the constraints on the location of the pumping wells given by Eqs. (4), (5) and (6), respectively. Constraint in Eq. (4) ensures that there is minimum separation between two production wells, and a safe distance is maintained between the production well and any poor quality water zone in the study area, or any restricted area where the well should not be installed. Equation (5) and (6) ensures that the pumping wells are located within the boundary of the study area. If the wells are located close to each other or near the boundary of the study area, wells may not be able to realize their full freshwater yield potential due to overlapping of the cone of influences, and overall reduce the freshwater production from the coastal aquifer. Maintaining a safe distance from poor quality water zone would ensure that the pumped water is free from contamination.

Simultaneous Pumping and Location Optimization Using Coupled
Sharp-interface and Density-dependent Model Pool and Carrera (2011) and Christelis and Mantoglou (2016) have shown that integrating the effects of dispersion and diffusion in the sharp-interface model enhances its prediction accuracy. This concept is utilized in this study by employing a density-dependent model to include the effects of dispersion and diffusion into a sharp-interface model, therefore improving precision and computation speed. The density-dependent model is used to modify the density ratio of shap interface model. Sharp-interface model is a function of head and density ratio. Since, the head is governed by the pumping rate and the pumping well locations, the density ratio is modified within an optimization framework (PSO) such that the outcomes of the sharp-interface model is same as that from a density-dependent model. This is achieved by minimizing the residual distance τ d (Fig. 1) i.e. toe location of the of sharp-interface model and the corresponding isosalinity contour from density-dependent model given by Eq. (7) subject to constrain given in Eq. (8).
where x i is the distance of the toe location of sharp-interface from the coast at i th well, x ci is the distance of the toe location of isosalinity contour from the coast at the i th well.
The developed approach is a three-step iterative procedure. The sharp-interface model is linked with SA optimization in the first step to maximize pumping and simultaneously optimize the location of production wells. In the second step, the results (optimized pumping values and the corresponding well locations) are used as input in the density dendent model to calculate the iososalinity contour. In the last step, an optimal density ratio is calculated in an optimization framework that minimizes the residual distance τ d (Eq. (7)). This optimal density ratio is again used in the sharp-interface model in step 1. This porcess continues until the total pumping obtained from the sharp-interface model stabalizes and there is no significant change from one iteration to the other. This iterative use of coupled sharp-interface and density-dependent models is described in schematic chart shown in Fig. 2.
In a classical linked simulation optimization, the simulation model (in this case, SEA-WAT) would have to be iteratively evaluated for each candidate solution to find the most optimal pumping rates and the pumping locations. As density-dependent models are computationally expensive, running several iterations of the same would make the solution computationally infeasible. The more precise SEWAT simulation model is substituted by a more efficient sharp-interface model in this developed technique. This efficiency comes at the expense of precision since sharp-interface models do not account for the effect of hydrodynamic dispersion. To compensate for this inaccuracy, the density ratio is modified until the toe location generated by the sharp-interface model aligns with the toe location of the isosalinity contour from SEWAT model. The specific advantage of this method is that it is no longer necessary to run the density-dependent model iteratively but only once to estimate the isosalinity contour. Therefore, this method, in principle can significantly decreases the computing cost while generating results that are as accurate as densitydependent model.

Performance Evaluation of the Developed Methodology
The performance of the developed methodology is evaluated against two hypothetical study area and a case study of coastal aquifer of Puri City (Fig. 3).

Description of Hypothetical Study Areas
The hypothetical study areas comprises of a heterogeneous, isotropic three dimensional unconfined costal aquifer having an area of 2.0 × 10 6 m 2 (1000 m × 2000 m × 15 m). The study areas are discretized into finite difference grids (10 m × 10 m × 3 m). The right side boundary is considered as the sea. The left hand boundary is land and considers as 2.5 m fixed head boundary (Fig. 3a, b). Therefore, at the start the head is considered linearly varying from land side at 2.5 m to the seaside at 0 m. The remaining two boundaries are considers as no-flow boundary. Initially, the salt concentration in the aquifer is assumed to be zero. Dirichlet boundary condition is assigned at seaside and landside with concentration 0 kg/m 3 and 35 kg/m 3 . It is assumed that the aquifer is pumped  Fig. 2 The developed iterative proses uniformly throughout the depth of the aquifer at the designed well locations. The aquifer is studied for a period of 10 years. The details of the hydrogeological parameters used in the study scenarios are presented in Table 1.
In the two hypothetical study scenarios considered, six number of pumping wells need to be installed for pumping freshwater. The locations of the pumping wells along with the pumping rates will be optimally obtained from the developed methodology. It is assumed that the installed wells will be uniformly pumped throughout the stress period. The upper and lower limit of the pumping is considered as 1000 m 3 /day to 50 m 3 /day respectively, for each well. To ensure that the installed well do not interfere with each other and contribute to its full potential a minimum distance of 350 m is kept between the two wells and the wells are 350 m from the boundary. Finite difference

Description of the Coastal Aquifer of Puri City
Puri City is situated in the east coast of India spread between 19 • 47 ′ to 19 • 50 ′ N latitudes and 85 • 48 ′ to 85 • 52 ′ E longitudes (Fig. 3c) and has hot humid climate. The topography of the area slopes gently towards the sea. The maximum height at the north eastern side is 18.4 m. The aquifer is sandy unconfind and extends below 40 m form Mean Sea Level (MSL). Puri City receives an average of 1,517 mm precipitation annually. Most of the preciptation happens in the four months of monsoon (July to September) (Vijay et al. 2011a). Groundwater is the main source of freshwater for Puri City because rivers near the city are intermitent in nature and not able to supply required amount of freshwater. The city is supplied with 20.5 × 10 3 m 3 /day through 36 pumping wells (Fig. 3c) from two pumping field in the city i.e. Chakratirtha and Balia Panda. There are 19 pumping wells in Chakratirtha and 13 wells in Balia Panda water field. There also four wells spread across the city. Apart from the pumping wells, there are several hand pumps. The details of the hydrogeological parameters used in the study scenarios are presented in Table 1. MODFLOW is used to simulate the flow model and SEAWAT is used to calculate salt concentration. The aquifer is divided with 30 m × 30 m finite difference grid with single layer. A fixed head value of 0 m and constant salt concentration of 35 kg/m 3 is taken along the coast. For the landside, the constant salt concentration of 0 kg/m 3 is taken. The head at the boundary is determined based on the contour from water table during post-monsoon of 2006 and premonsoon of 2007 provided by Vijay and Mohapatra (2016). Though, there are 19 production wells in Chakratirtha water field and 13 production wells in Balia Panda water field, only eleven production wells in Chakratirtha water field and three pumping well in Balia Panda water field are considered in the finite difference model. This is done in case when more than one well is located in one discretized grid or are located very close to each other. Therefore, final well count in the finite difference model for Puri City aquifer is 17 (11 in Chakratirtha water field, 3 in Balia Panda water field and 4 wells spread across the field) (Fig. 3c).  Vijay et al. (2011b) reported that total pumping from the aquifer is 20.5 × 10 3 m 3 /day from 36 pumping wells. Since the pumping from each well is not reported therefore, the pumping is adjusted to calibrate the head values. Pre and post monsoon water table information of 13 well locations are used to calibrate the model with head interval of 1 m and confidence level is 95%. The calibration results for the groundwater flow model are shown using box plots (Fig. 4). The absolute difference between the estimated and observed heads was well within the appropriate target interval, as indicated by the green box plots. Since no information is available for the initial position of the saltwater-freshwater front, therefore, modified Ghyben-Herzberg density ratio (Dey and Prakash 2019a) is used to extrapolate the initial saltwater-freshwater front. The calibrated model is run for 30 years with the current rate of pumping to realize the impact of saltwater intrusion in the Puri coastal aquifer (Fig. 5). It can be clearly seen that at current rate of pumping there would be saltwater intrusion by year 2037 and some of the existing wells will be rendered useless. Hence, it is necessary to optimize the pumping as well as find optimal locations for new wells to be installed to augment the existing supply of freshwater to meet the increasing demand. The locations of the new wells to be installed is further constrained by an existing poor water quality zone due to pollution from anthropogenic sources (Fig. 3c).

Evaluation Results and Discussion
The iterative process of pumping optimization followed by density ratio optimization is run 15 and 16 times for scenarios 1 and scenario 2 respectively, to find the optimal pumping well locations and pumping rates from the installed wells. Corresponding density ratio is also calculated that would reflect the effect of dispersion and diffusion in a sharp-interface simulation model, thereby giving a more realistic estimate of the actual scenario without causing saltwater intrusion. The performance evaluation results of the optimal pumping value and the corresponding density ratio is presented in Fig. 6. It is observed that after eight iterations for scenario 1 and seven iterations for scenario 2, there is no significant change in values of total daily pumping and density ratio. The maximum total pumping for scenario 1 at eighth iteration was 2060.22 m 3 /day, and for scenario 2 the maximum total pumping after seventh iteration was 1560.88 m 3 /day. The maximum total pumping in case of scenario 1 is more than scenario 2. This could be explained by the locations of the pumping wells. Since the location of the pumping wells are not the same in both the cases, therefore, the maximum total pumping also varies. The methodology finds different optimal locations for wells due to varying hydrogeological properties in the two scenarios. Optimal locations of each well, isosalinity contour corresponding to a salinity level of 5 kg/m 3 salt concentration, and toe of sharp-interface is plotted in the Fig. 7. It is evident from the figure that the optimal locations of the wells found by the methodology satisfies the required criteria of minimum distance for each other i.e. 350 m and a minimum distance of 350 m from the boundary. This ensures that there is no interference between the cone of influence of one well with the other and the wells can be pumped to its full potential without causing saltwater intrusion. The toe location from the sharp-interface is in agreement with the corresponding isosalinity contour (salinity level of 5 kg/m 3 salt concentration) obtained from SEA-WAT. This shows that optimizing the density ratio can improve the accuracy of prediction of the saltwater-freshwater front in a sharp-interface model without compromising the efficiency.
The performance evaluation results of pumping optimization and the corresponding density ratio for every iteration for Puri coastal aquifer is shown in Fig. 8. Maximum pumping value is achieved after four iterations with a total pumping of 57.3 × 10 3 m 3 /day and corresponding density ratio of 0.000992. The iterative process is run for ten iterations to ensure that the optimization does not converge to a local optima giving suboptimal results. The maximum total pumping and density ratio did not show much significant change in further iterations and was stable up to 10 iterations. The final well locations, location of toe of sharp-interface and isosalinity contour are plotted in Fig. 8b. It is clearly evident from the Fig. 8b that none of the wells located by this methodology is intruded by saltwater. Toe location from the sharp-interface is found to be in close agreement with the isosalinity contour from SEAWAT. Results obtained from this method shows that by optimally orienting the well locations the total freshwater extraction can be increased almost 3 times of the current extraction rate (20.5 × 10 3 m 3 /day). This methodology also ensures that the well locations found have a minimum distance of 750 m between any two wells and a minimum of 300 m from the boundary of the study area. The methodology also ensures that the extraction wells are such positioned that they maintain a safe distance of 90 m from the poor quality water zone. Wells no i.e. well-10, well-6, well -15, well-17 and well-4 as seen in Fig. 8b are either very close to the poor quality water zone or the saltwater freshwater interface. Therefore, the respective distance of these wells from interface or the boundary of the poor water quality zone is calculated (Table 2) to check if there is any violation in terms of placing the wells. The calculated distances show that all the wells are safely located. Correlation between the estimated freshwater-saltwater front and the isosalinity contour from density-dependent model are statistically compared (Fig. 8c). The R 2 value (R 2 = 0.873) clearly indicates that the freshwater-saltwater front predicted by sharp interface model closely match with the isosalinity contour from densitydependent model.

Conclusions
A coupled sharp-interface linked to SA optimization and density-dependent model linked to PSO is developed to maximize total pumping from a coastal aquifer without causing saltwater intrusion and simultaneously find the optimal well locations. The developed methodology performs satisfactorily in optimizing pumping well locations and maximizing total pumping from a coastal aquifer without causing saltwater intrusion. The coupled approach reduces the computational effort by eliminating the need to run SEAWAT model to evaluate each candidate solution. The developed methodology clearly establishes that location of the pumping wells have a significant role in deciding the maximum total pumping from a coastal aquifer. It is also seen that the locations of the production wells are sensitive to hydrogeological settings and impact the total pumping from the aquifer. It can be concluded that the developed methodology is a viable alternative in managing coastal aquifer with greater accuracy compared to linked simulation optimization with sharp-interface model for limited computational budget. Performance evaluation results show that the method can simultaneously optimize both the locations of the pumping wells as well as pumping rates as tool for managing saltwater intrusion. Constrains of location in design of any well field is quite common and often restrictive in terms of placement of production wells. This developed methodology is robust to handle such restrictive constraints and still find the optimal pumping rates to manage saltwater intrusion. While taking costal aquifer of Puri City as a proof of concept, it is clear that this methodology can potentially increase the total freshwater pumping by almost 3 times the existing pumping without causing saltwater intrusion by optimally choosing the well locations. The methodology in its current form can be directly applied for new well fields design for sustainable freshwater extraction from coastal aquifers. The method can also be applied to scenarios where the existing freshwater extraction is to be increased by augmenting new wells in coastal aquifer without causing salt water intrusion. The developed methodology can bolster the freshwater supply in meeting the future demands without causing saltwater intrusion.