Optimal Pollution Control and Pump-and-Fertilize Strategies in a Nitro-Polluted Aquifer, Using Genetic Algorithms and Modﬂow

: Nitro-pollution in a conﬁned aquifer may originate from its recharge area, e.g., agricultural sites, animal feedlots, septic tanks, and other waste disposal sites or from treated wastewater recharge wells. The latter case is studied. Existing water supply pumping wells should be protected for a given period. Instead of typically investigating optimal pump-and-treat or hydrodynamic pollution control management/remediation strategies, a novel combined pollution control and pump-and-fertilize (PAF) approach is proposed: protect existing wells with additional wells, considering pumped nitro-polluted groundwater as proﬁtable reusable fertilizer rather than a pollutant to be remediated; convey pumped polluted water to an irrigation reservoir, considering nitrogen (N) uptake by irrigated crops in nearby farmlands and proportional decrease in fertilizer application, meaning proﬁt. Optimization entails the operation of optimally located additional pumping wells with optimal ﬂow rates, minimizing the sum of (i) annual pumping cost, (ii) pipe network (connecting additional wells and reservoir) amortization cost, and (iii) proﬁt from retrieved N reuse. Modﬂow simulates a 3D ﬂow ﬁeld and advection-dispersion mass transport, while Genetic Algorithms (GAs) handle optimization. Various scenarios are simulated concerning crops’ retrieved N root uptake percentage, fertilizer, and energy market prices. The paper provides a data-ready optimization/decision support tool that creates a pool of alternative (sub)optimal management solutions/strategies.


Introduction
Groundwater nitrate pollution is a serious environmental issue. It is primarily related to human activities such as agricultural practices and the intensive use of nitrogen-based fertilizers or manure, sewage, industrial wastewater discharge, and leachates from landfill sites [1]. Other sources can be the overflow of septic tanks and mine tailings where ammonia-based explosives are used for detonation [2]. Nitrates can cause significant harm to human health and the environment, and their presence in groundwater can have longterm effects. Nitrate pollution has become a major concern in many parts of the world. Various strategies and methods are being developed to prevent or mitigate its impact. The nitrate in groundwater is usually remediated by either separation (electrocoagulation, nanofiltration, electrodialysis, electrodeionization, capacitive deionization, adsorption) or reduction-based approaches (chemical reduction, biological reduction) [3]. Especially in agricultural areas, the extensive use of fertilizers is the main source of nitrate pollution [4][5][6]. A promising agricultural practice that can significantly reduce the fertilizer quantities used by farmers is by accounting for the nitrogen (N) applied through their irrigation water. In this context, the "pump-and-fertilize" (PAF) approach is used to pump water from nitrateaffected wells and use it as a fertilizer in the irrigation system. However, many farmers use high water levels and N to meet plant requirements. Thus, the soil becomes vulnerable to leaching and transferring N to groundwater resources [7]. Because nitrate concentrations in pumped groundwater vary considerably in space and time, frequent monitoring and adaptive nutrient management are required [8]. Monitoring the N concentration in pumped groundwater and adjusting fertilizer applications accordingly to account for the quantity of N applied with the irrigation water could achieve a significant reduction of aquifer nitrate pollution risk [9]. The use of PAF and optimization models could reduce the consumption of N-based fertilizer, even by 88%, in a study in an agricultural area in Iran [10]. Irrigation is also critical in applying a PAF method, e.g., a reduced irrigation rate (490 mm from 750 mm) facilitated the use of approximately 42 kg N ha/yr of nitrates from groundwater without yield loss. It also gradually improved the groundwater quality in a corn-cultivated area [11]. However, the efficiency of utilizing nitrate-rich groundwater is related to optimal irrigation water and fertilization management. Therefore, technical guidelines should be provided to farmers by the regional agricultural policymakers.
In this work, the PAF approach is combined with nitrate pollution control for the optimal management of a nitrate-contaminated aquifer using a GA. The study case is a theoretically confined aquifer where the pollution could be the result of (a) long-term use of fertilizers due to intensive agricultural exploitation (diffuse source nitro-pollution) or treated wastewater discharge (point-source nitro-pollution) in the remote recharge area of the aquifer (the impact in the studied area would be similar), or (b) nitrogen-rich pointsource treated wastewater discharge directly into the studied confined aquifer (Figure 1). To increase the difficulty level and provide the most challenging pollution scenario for the optimization method, scenario (c) is presented in this paper. Two Gaussian pollution plumes are introduced. They approximate the result of long-term point-source leakages but are also specifically engineered to increase the optimization difficulty level. There are high (plume centers) and low concentration areas. The practical goal is to use additional wells (Aws) to protect the operation of an existing water supply pumping well scheme for a given period. Simultaneously, the pumped nitro-polluted water can be exploited for irrigation purposes through PAF strategies instead of pump-and-treat (PAT) or hydraulic control (HC) ones. The mathematical goal is the minimization of the aquifer's management cost, using Gas, protecting the existing water supply wells, considering: (a) pumping cost, (b) pipe network amortization cost (conveying pumped pollutes water from the AWs to a reservoir), and (c) retrieved N reuse profit. The last cost (profit) item uses pumped nitropolluted water (by the AWs) to irrigate nearby regional farmlands. The N from the polluted water that can be uptaken and assimilated by the plants can proportionally decrease the respective fertilizer application. This can save money and lead to a profit. A first approach was presented by Kontos (2020) [12], where the flow field was calculated analytically, and only advective mass transport was simulated with particle tracking. This research was itself based on the foundations of earlier research ("OptiManage" v4 software suite for optimal groundwater resources management, e.g., Kontos and Katsifarakis, 2012 [12]; 2017 [13]). There, the pollutant in the water, if pumped by AWs, should be treated, and the respective cost item was loss and not profit. The current paper drives the concept further: it uses finite-differences free software (Modflow 6; ref. [14]) while including mechanical dispersion in the mass transport mechanisms. Vertical section (sketch) of the theoretical studied confined aquifer for scenario "a" (diffuse agriculture-originated nitro-pollution in the remote recharge area) and scenario "b" (point-source pollution by wastewater discharge wells directly in the studied confined aquifer). Scenario b is presented in this paper as it is more computationally challenging for the optimization method.

Theoretical Problem and Study Field
The 2km x 2km top view of the theoretical flow field is presented in Figure 2, together with the initial nitrates' concentration distribution, the two existing water supply pumping wells (Wex,1 and Wex,2), and an irrigation reservoir (RES).
Nitrate concentration is assumed to be uniformly distributed with depth, that is, in all three layers of the confined aquifer. The initial nitrate concentration information is presented in Supplementary Material Files SM1 and SM2. The aquifer is assumed to have been polluted by two-point sources of wastewater. Specifically, two treated wastewater discharge wells are considered to have injected nitrate-rich wastewater directly into the studied confined aquifer (Figure 1, pollution scenario b). A series of tests were initiated concerning the form and dimensions of the plumes. A set of two approximate Gaussian plumes, one larger in size and concentration, mildly stretched along the y-axis direction, and a smaller one in size and concentration, strongly stretched along the x-axis direction, were used as initial conditions of the pollutant concentration. The plumes' form in the initial conditions ( Figure 2) does not accurately represent the expected form, given the assumed flow regime (i.e., initial and boundary conditions concerning the hydraulic head and layout of the pumping scheme). Nonetheless, this initial pollution scheme was selected in this work as it provides the most challenging case to test the performance of the optimization algorithm on the PAF strategies. Vertical section (sketch) of the theoretical studied confined aquifer for scenario "a" (diffuse agriculture-originated nitro-pollution in the remote recharge area) and scenario "b" (point-source pollution by wastewater discharge wells directly in the studied confined aquifer). Scenario b is presented in this paper as it is more computationally challenging for the optimization method.

Theoretical Problem and Study Field
The 2 km × 2 km top view of the theoretical flow field is presented in Figure 2, together with the initial nitrates' concentration distribution, the two existing water supply pumping wells (W ex,1 and W ex,2 ), and an irrigation reservoir (RES).
Nitrate concentration is assumed to be uniformly distributed with depth, that is, in all three layers of the confined aquifer. The initial nitrate concentration information is presented in Supplementary Material Files SM1 and SM2. The aquifer is assumed to have been polluted by two-point sources of wastewater. Specifically, two treated wastewater discharge wells are considered to have injected nitrate-rich wastewater directly into the studied confined aquifer (Figure 1, pollution scenario b). A series of tests were initiated concerning the form and dimensions of the plumes. A set of two approximate Gaussian plumes, one larger in size and concentration, mildly stretched along the y-axis direction, and a smaller one in size and concentration, strongly stretched along the x-axis direction, were used as initial conditions of the pollutant concentration. The plumes' form in the initial conditions ( Figure 2) does not accurately represent the expected form, given the assumed flow regime (i.e., initial and boundary conditions concerning the hydraulic head and layout of the pumping scheme). Nonetheless, this initial pollution scheme was selected in this work as it provides the most challenging case to test the performance of the optimization algorithm on the PAF strategies. The two existing pumping wells (EWs) must uninterruptedly supply fresh drinking water to an adjacent settlement at a constant flow rate of ΣQex = 200 L/s for 300 days. The computationally greedy nature of GAs requires a certain balance between computational accuracy and efficiency [15]. Thus, this pilot approach studies a simplified 3D flow field in a homogenous, isotropic, confined aquifer with a plane, single-phase, steady flow. The hydraulic conductivity of the aquifer is Kx = Ky = Kz = 8.64 m/d, aquifer thickness is b = 30 m, and efficient porosity ne = 0.2, throughout its extent.
It is assumed that currently, the total flow rate that is required to be pumped by the two EWs is evenly distributed (Qex,1 = Qex,2 = 100 L/s). If no counter-measures are taken, the pollution plumes will spread and pollute (be pumped by) the existing wells within 300 days, when the existing pumping scheme should be protected. This is proved by a simple simulation, presented in Figure 3a (colored map of nitrate concentration during the last 300th day of simulation) and Supplementary Materials SM3 (video of the full pollution spread progression). The nitrate pollution mass transport is simulated using Modflow (see Section 2.2 for details).  The two existing pumping wells (EWs) must uninterruptedly supply fresh drinking water to an adjacent settlement at a constant flow rate of ΣQ ex = 200 L/s for 300 days. The computationally greedy nature of GAs requires a certain balance between computational accuracy and efficiency [15]. Thus, this pilot approach studies a simplified 3D flow field in a homogenous, isotropic, confined aquifer with a plane, single-phase, steady flow. The hydraulic conductivity of the aquifer is K x = K y = K z = 8.64 m/d, aquifer thickness is b = 30 m, and efficient porosity n e = 0.2, throughout its extent.
It is assumed that currently, the total flow rate that is required to be pumped by the two EWs is evenly distributed (Q ex,1 = Q ex,2 = 100 L/s). If no counter-measures are taken, the pollution plumes will spread and pollute (be pumped by) the existing wells within 300 days, when the existing pumping scheme should be protected. This is proved by a simple simulation, presented in Figure 3a (colored map of nitrate concentration during the last 300th day of simulation) and Supplementary Materials SM3 (video of the full pollution spread progression). The nitrate pollution mass transport is simulated using Modflow (see Section 2.2 for details). Figure 3b specifically presents the (pumped) nitrate concentration (mg/L) vs. time (d) diagram in each EW.
According to the studied problem, up to two AWs (W ad,1 and W ad,2 ) can protect the EWs from nitro-pollution for 300 days. The notion of N reuse for irrigation is implemented in this paper through the following simple concept: the volume of pumped water by the AWs is collected in an existing irrigation reservoir ("RES" in Figure 2). This re-distributes water through an existing irrigation network. In this way, the seasonal irrigation needs of the regional farmlands (here 8 × 10 6 m 3 ) can be fully satisfied. According to the studied problem, up to two AWs (Wad,1 and Wad,2) can protect the EWs from nitro-pollution for 300 days. The notion of N reuse for irrigation is implemented in this paper through the following simple concept: the volume of pumped water by the AWs is collected in an existing irrigation reservoir ("RES" in Figure 2). This re-distributes water through an existing irrigation network. In this way, the seasonal irrigation needs of the regional farmlands (here 8 × 10 6 m 3 ) can be fully satisfied.
The conceptual model of N mass transport is based upon the following assumptions: • All N-based compounds from fertilizing the irrigated parts of the studied confined aquifer's recharge area have been converted to NO3− (similar to [16]).

•
The mass transport mechanisms only include advection and mechanical dispersion.

•
It is assumed that aerobic conditions apply in the aquifer, with dissolved oxygen concentrations DOC > 0.5 mg/L (oxic groundwater; similar to [17]), which do not favor the denitrification process (as described in [18,19]).

•
Nitrates in the aquifer are assumed non-reactive, and retardation is not considered; no additional chemical reaction (e.g., adsorption, desorption) is simulated (similar to [20]). Groundwater residence time or age is not considered.
To realistically investigate the studied environmental problem, various values of financial-management-related factors are considered: The conceptual model of N mass transport is based upon the following assumptions: • All N-based compounds from fertilizing the irrigated parts of the studied confined aquifer's recharge area have been converted to NO 3− (similar to [16]).

•
The mass transport mechanisms only include advection and mechanical dispersion. • It is assumed that aerobic conditions apply in the aquifer, with dissolved oxygen concentrations DOC > 0.5 mg/L (oxic groundwater; similar to [17]), which do not favor the denitrification process (as described in [18,19]).

•
Nitrates in the aquifer are assumed non-reactive, and retardation is not considered; no additional chemical reaction (e.g., adsorption, desorption) is simulated (similar to [20]). Groundwater residence time or age is not considered.
To realistically investigate the studied environmental problem, various values of financial-management-related factors are considered:

•
Low and high energy cost values (€/kWh) affect the pumping cost, including possible subsidization or discounts. • Low and high N root uptake values, namely the percentage of N mass in the nitratepolluted water pumped by the AWs and irrigated to crops assimilated by the crops, affect the profit by pumping nitro-polluted water with AWs.
• Low and high N mass cost values (N-based commercial fertilizers' market price) or final values after possible subsidization. Three scenarios of combinations of the above factors are investigated, presented in Section 3.

Hydraulic Simulation with Modflow 6
Simulating the flow field and mass transport are assigned to the Modflow 6 software of USGS [14]. It is an established finite differences tool. However, the need to automate the simulation and link the surrogate model with an optimization tool (GAs) of an iterative nature leads to the use of the "Flopy: Python Package for creating, running and postprocessing MODFLOW-Based Models" [25][26][27][28][29]. The 2 km × 2 km flow field is discretized in a square grid of 80 × 80 cells (2D top view), each with a 25 m side. The 30 m thick confined aquifer comprises three layers, 10 m thick each. The three layers provide a 3D model that adds to the accuracy of flow field calculation and mass transport simulation. It is also the only technical way in Modflow 6 to realistically simulate the fact that pumping wells do not usually pump along their full penetration length in the confined aquifer but through a well screen. Thus, the mid-layer is where all pumping wells sustain their well screens ( Figure 1). The initial piezometric surface is assumed to be 101.5 m high. The northern and southern boundaries are impermeable, while the western and eastern ones are constant head boundaries in the form of General Head Boundaries (GHB; [29]). Their hydraulic heads equal 103 and 100 m, respectively. This way, a West-East natural flow is described due to the hydraulic gradient of (103 m − 100 m)/2000 m = 1.5‰. The GHBs are assumed to be at a distance of 1 km from the physical field boundaries, with an intermediate conductivity C equal to [30]: where: C is the hydraulic conductivity of the intermediate layer between GHB and respective Western or Eastern physical boundary (C = 2.5 × 10 −5 m 2 /s); K is the hydraulic conductivity of confined aquifer (K x = K y = K z = 8.64 m/d); W is the vertical cell dimension (10 m); L is the horizontal cell dimension (25 m); M is the thickness of layer between GHB and respective Western or Eastern physical boundary (1000 m).
No chemical reactions are considered; molecular diffusion is assumed to be negligible. The longitudinal, transverse, and vertical dispersity values are assumed to be α L = 10 m, α T = 1 m, α V = 0.1 m, respectively, following previous research [24,31].

Objective Function
The studied problem is a classic constrained, nonlinear, stochastic, multi-objective optimization problem [32]. Specifically, it is a problem of minimization of Fitness Value (FV). FV is the total polluted aquifer management cost of each solution. The objective function to be minimized (including penalty imposition) is: Find: Q i , Q j , X j , Y j , i = 1,2, j = 3,4 so that FV = VB1 + VB2 + VB3 + Penalty = Min, (2) Subject to the constraints: C2: 0 L/s ≤ Q j ≤ 230 L/s, j = 3,4, C3: 0 m ≤ X j ,Y j ≤ 2000 m, j = 3,4, where: Q i are the flow rates of all four pumping wells (L/s); X j and Y j are the coordinates of the two additional pumping wells (in m; wells 1 and 2 are the EWs, while wells 3 and 4 are the AWs); FV is the Fitness Value (total cost including possible Penalty values; in €/yr); VB1 is the total annual pumping cost for all pumping wells (€/yr); VB2 is the amortization cost of pipe network that conveys pumped water from AWs to RES (€/yr); VB3 is the profit by the decreased use of commercially purchased N-based fertilizer, which is equal to the current market price of N mass retrieved from the nitro-polluted confined aquifer and reused as fertilizer (€/yr); Penalty is the total penalty imposed on a chromosome (proposed solution) in case of violation of one or more constraints (see Section 2.3.4); C i , i = 1,2 are the nitrates' concentrations in the respective EWs 1 and 2 (mg/L). Actually, due to the numerical method used to simulate the flow field and mass transport, the coordinates of the wells' locations are not accounted for in m. They are accounted as nominal numbers, representing consecutive columns and rows. Thus, constraint C3 is modified:

Genetic Representation, Chromosome Structure and GA Parameters
GAs are applied here with a binary genetic representation; chromosomes/proposed solutions are practically binary strings. Each chromosome represents AWs' coordinates (or rather several columns and rows), flow rates, and one of the EWs (since the total required flow rate by the existing pumping scheme is fixed). Each flow rate and coordinate can obtain certain max values in the decimal numeral system. Hence, their binary form defines the max number of digits (string length) each chromosome part requires. This way, the total number of digits (total chromosome/string length) can be calculated. According to C3 , as far as AWs wells are concerned, it is 1 ≤ X j ,Y j ≤ 80 (nominal; the number of columns and rows) and 0 ≤ Q j ≤ 230 (L/s) for j = 3,4 (AWs). As far as EWs are concerned, it is: Q 1 + Q 2 = 200 L/s, so it can marginally be accepted that: 0 ≤ Q 1 ,Q 2 ≤ 200 (L/s). Table 1 presents the maximum values of the decision variables in decimal and binary numeral form, the number of digits (string length) of the chromosome part that represents each decision variable, and the decimal value of the max binary number that the specific string part can obtain. For example, the length of each string part representing the consecutive number of a column/row of an AW is SLX = 7. Thus, the length of the part of all coordinates of AWs is SL1 = 4 × SLX = 28. The length of the part that represents the flow rate of W ex,1 is SL2 = SLQ ex = 8. The length of each part that represents a flow rate of an AW is SLQ = 8. Thus, the length of the part of all flow rates of AWs is SL2 = 2 × SLQ = 16. The total length of a typical binary chromosome is SL = SL1 + SL2 + SL3 = 28 + 8 + 16 = 52. Its structure is presented in Figure 4.  The genetic operators used are (a) selection, featuring the tournament procedure and elitism for the fittest chromosome, (b) crossover, and (c) mutation. The chromosome population is PS = 60, the number of generations is NG = 500, the selection constant is SC = 3, and the crossover probability is CRP = 0.4. Regarding the mutation probability (MP), the proposed methodology does not follow the old empirical rule proposing that 1/SL (here 1/52 ≈ 0.019) is a good first value for MP. Instead, it follows own previous research (Kontos and Katsifarakis, 2012 [12]) suggesting that MP values ranging from 2/SL to 2.5/SL (here 2/52 ≈ 0.038 to 2.5/52 ≈ 0.048) are more efficient in converging to the optimal solution. Finally, the value MP = 0.04 is used.

Pumping Cost VB1
Pumping cost VB1 is the market price of the energy needed to extract the desired groundwater volumes in the desired period (ΣQex = 200 L/s for 300 days). It is the main cost item in groundwater resources management problems, given by: where: P is the required power for pumping (kW); Tp is the time duration of pumping (h); CkWh is the cost of one kWh (€/kWh). Energy prices are volatile at the time of current research, with an increasing energy inflation rate. A current realistic value during the winter of 2021 in Greece for a business/professional pack by Public Power Corporation (PPC; "C21" pack [33]), considering the Adjustment Clause and Energy Transition Fund [34][35][36] is CkWh = 0.28244 €/kWh. The impact of the variation of this parameter on the optimal strategies produced is investigated. For this reason, the value CkWh = 0.06641 €/kWh is also used, representing the real cost before the economic crisis in 2008. Alternately, it represents the reduced cost if the energy price was subsidized/discounted for agriculture/irrigation. The genetic operators used are (a) selection, featuring the tournament procedure and elitism for the fittest chromosome, (b) crossover, and (c) mutation. The chromosome population is PS = 60, the number of generations is NG = 500, the selection constant is SC = 3, and the crossover probability is CRP = 0.4. Regarding the mutation probability (MP), the proposed methodology does not follow the old empirical rule proposing that 1/SL (here 1/52 ≈ 0.019) is a good first value for MP. Instead, it follows own previous research (Kontos and Katsifarakis, 2012 [12]) suggesting that MP values ranging from 2/SL to 2.5/SL (here 2/52 ≈ 0.038 to 2.5/52 ≈ 0.048) are more efficient in converging to the optimal solution. Finally, the value MP = 0.04 is used.

Pumping Cost VB1
Pumping cost VB1 is the market price of the energy needed to extract the desired groundwater volumes in the desired period (ΣQ ex = 200 L/s for 300 days). It is the main cost item in groundwater resources management problems, given by: (8) where: P is the required power for pumping (kW); T p is the time duration of pumping (h); C kWh is the cost of one kWh (€/kWh). Energy prices are volatile at the time of current research, with an increasing energy inflation rate. A current realistic value during the winter of 2021 in Greece for a business/professional pack by Public Power Corporation (PPC; "C21" pack [33]), considering the Adjustment Clause and Energy Transition Fund [34][35][36] is C kWh = 0.28244 €/kWh. The impact of the variation of this parameter on the optimal strategies produced is investigated. For this reason, the value C kWh = 0.06641 €/kWh is also used, representing the real cost before the economic crisis in 2008. Alternately, it represents the reduced cost if the energy price was subsidized/discounted for agriculture/irrigation. Power is given by: where: ρ is the density of the pumped fluid (kg/m 3 ); g is the gravitational acceleration (m/s 2 ); n p is the pumps' efficiency (dimensionless); NW is the total number of pumping wells (nominal); Q i is the flow-rate of pumping well i (in L/s); s i is the hydraulic head drawdown in the side of pumping well i (in m); δ is the distance between the initial horizontal level of the hydraulic head and the predefined reference level.
Since δ is the same everywhere, VB1 could be equal to the variable term of pumping cost, namely: where: A is a coefficient depending on the density of the pumped fluid, the pump's efficiency, the pumping duration, the gravitational acceleration, and the cost of a kWh, given: (11) where: T p is the study period in h (here T p = 1 yr = 12 × 30 × 24 h = 8640 h); ρ H2O = 997 kg/m 3 is the density of water; g = 9.81 m/s 2 ; n p = 0.8; for

Pipe Network Cost VB2
Pipe network optimization alone is an extremely complicated problem. It is a complex graph theory optimization problem, approached by meta-heuristics or machine learning methods. These cannot fit into the objective function of a meta-heuristic method (GAs) of an iterative nature itself. For this reason, the pipe network optimization process is reduced into a minimum spanning tree problem [37]: VB2 cost minimization is approximated by minimizing the total length of the pipe network alone. Pipes' diameters are only indirectly considered, as they are fixed to the respective flow rate. Therefore, VB2 practically represents the pipe network amortization cost, directly proportional to the initial network cost (e.g., [12,13,38,39]).
The initial network cost is assumed to be directly proportional to the total length of the pipes and their diameters. The pipes convey groundwater pumped by the AWs to the irrigation reservoir (X RES = 1312.5 m; Y RES = 1937.5 m; Figure 2). The pipe network's construction initial cost is assumed to be 45 €/m for small pipe diameters and 60 €/m for large pipe diameters. Practically, the pipe diameter for each section is selected depending on the respective flow rate. The flow-rate threshold that indirectly characterizes a pipe diameter as small or large is set at Q thres = 50 L/s. For an amortization period of 10 years and a 5% annual interest rate, the respective annual amortization cost is A a1 = 6 €/m and A a2 = 8 €/m, calculated by: where: i is the type of pipe (1 for small or 2 large diameter; nominal); C i is the initial network cost for pipe type i (C 1 = 45 €/m for Q < Q thres ; C 2 = 60 €/m for Q > Q thres ); r is the annual interest rate (here, r = 5%); n is the amortization period (here, n = 10 yr). This simplified network minimization process contributes to the more realistic simulation of the theoretical optimal groundwater resources management problem. It adds the pipe network cost item to the total management cost, which is antagonistic to the pumping cost. Nonetheless, it does not compromise computational efficiency. Finally, VB2 is given: (13) where: N adW is the number of AWs (here N adW = 2); L i is the length of the pipe section that carries water away from AW i (in m); A ai is the annual amortization cost (in €/m/yr; see Equation (12)).

Nitrogen Retrieval/Reuse Profit VB3
VB3 cost is the profit from the retrieval and reuse of N from pumped nitro-polluted groundwater by the AWs as fertilizer. The irrigation with water full of nitrate, depending Agronomy 2023, 13, 1534 10 of 31 on the irrigation type (drip irrigation, sprinkler irrigation etc.), the root uptake of the plants and other climatic factors, leads to a partial satisfaction of N fertilizing needs. The profit comes from the proportional decrease of N-based fertilizer. It is calculated as the market price of the N mass that is not applied in the fields/farmlands. At each of the 30 consecutive 10-day time steps, the algorithm sums up the volume of the pumped groundwater by the AWs. Provided that it does not exceed the irrigation needs and reservoir capacity, the retrieved N mass is calculated: where: the negative sign in front of the equation of VB3 implies that the VB3 cost is actually profit and not loss; R N-NO3 = 14/(14 + 3 × 16) = 0.22581 is the ratio of atomic weights of N over the nitrates' elements (AW N = 14; AW O = 16), used to convert nitrate mass to N mass so that the profit per kg of N can be calculated; N uptake is the ratio describing the average percentage of N mass that can be absorbed (root uptake) by the plants in the farmlands, which will be irrigated by the irrigation reservoir containing retrieved nitropolluted groundwater pumped by the AWs (here, values N uptake = 0.70 and 0.07 are tested); N Cost is the current (winter 2021) market price of N in Greece (e.g., 40 kg of Ammonium Nitrate with 34.5 units of N cost 45 €), hence N Cost is approximately 3 €/kg; investigating the impact of fertilizer price increase (current trend) or of a possible subsidization of N reuse (pump-and-fertilize strategies) on the optimal strategy produced, the value of N Cost = 30 €/kg is also tested; N adW is the number of AWs (here, N adW = 2); M NO3 i is the mass of total nitrate pumped by (AW i during the whole study period (in kg; calculated by Modflow; NS is the number of timesteps of the study period (here, NS = 30); ∆t is the duration of each timestep (in d; here, ∆t = 10 d; total study period = NS × ∆t = 300 d); C N,ij is the concentration of N pumped by EW i during timestep j (in mg/L; calculated by Modflow); Q i is the flow-rate of AW i (in m 3 /d).

Constraint Handling
Constraint handling techniques in GAs can be generally categorized into methods of (a) penalty imposition for the infeasible solutions, (b) repair of infeasible solutions, (c) preservation of feasibility, (d) hybrid techniques [40]. Here, constraints C1-C4 (Equations (3)-(6)) are handled as follows: C1: The aggregate of the flow rates of EWs (W 1 , W 2 ) must be constant and equal to ΣQ ex = 200 L/s. Thus, only Q 1 is inserted into the typical chromosome. Q 2 is calculated by the relevant difference 200-Q 1 . Q 1 should obtain a max physical value of marginally 200 L/s. However, after the binary chromosome's decoding (binary to decimal), it can reach up to 255 L/s (  (Table 1). For this constraint, a category (a) constraint handling technique is implemented: if any coordinate acquires a value equal to zero (0), it instantly converts to a random integer number N i > 80. After that, for each solution that includes several coordinates with a value N i > 80, a dynamic penalty (Pen1) is imposed. Pen1 is proportional to the number of C3 constraint violations: where: C PEN1 is a coefficient for adjusting Pen1 values to reach the order of magnitude of other penalty items in the total penalty function (see Equation (17); here, C PEN1 = 1,000,000); N adW is the number of AWs (nominal); X j and Y j are the abscissa (number of column in the 80 × 80 Modflow top view grid) and the ordinate (number of rows in the Modflow grid) of AW j (j = 3, 4), respectively (nominal); maxX and maxY are the max physical values of X j and Υ j , respectively (nominal). C4: The EWs (W 1 , W 2 ) must not pump pollution during the study period of 300 d. Thus, the pollutant concentration in the respective cell must remain equal to zero. This constraint is handled with a category (a) constraint handling technique: the algorithm sums up the N mass pumped by EWs wells during the study period for each EW. If it is non-zero, a dynamic penalty (Pen2) is imposed. Pen 2 is proportional to the magnitude of the constraint violation (pumped N mass by EWs). Pen2 includes a constant term (C con PEN2 ). It also consists of a variable term that derives from the product of the pumped by EWs' nitrate mass (M NO3 i ) and the constant coefficient C var PEN2 . Pen2 is also proportional to the number of constraint violations (number of EWs polluted), as it is the aggregate of all separate penalty items that are imposed on each EW: where: R N-NO3 is the ratio of atomic weights of N vs. the nitrates' elements (R N-NO3 = 0.22581); N exW is the number of EWs (here, N exW = 2); C con PEN2 is the constant term of Pen2 (here, C con PEN2 = 2000); C var PEN2 is the coefficient of the variable term of Pen2 (here, C var PEN2 = 200); M NO3 i is the total mass of nitrate pumped by (existing) well during the whole study period (in kg); NS is the number of timesteps (here, NS = 30); C N,ij is the concentration of N pumped by (existing) well i during timestep j (in mg/L; calculated by Modflow); ∆t is the duration of each timestep (in d); Q i is the flow-rate of (existing) well i (i = 1,2; in m 3 /d).
The total penalty of each solution violating at least one constraint of C3, C4, is given: where: Pen1 is the Penalty item concerning solutions that entail a violation of constraint C3; Pen2 is the Penalty item relating to solutions that involve a violation of constraint C4; C f1 and C f2 are weighting factors of Pen1 and Pen2 in the total penalty, respectively (here, C f1 = 50,000; C f2 = 2000). The investigation for the optimal values of the penalty function parameters (C PEN1 , C con PEN2 , C var PEN2 , C f1 , C f2 ) requires an extensive series of tests. The criteria that must be met are: 1.
Pen1 parameter (C PEN1 ) value must ensure that the penalty term that deals with solutions violating constraint C3 (no AW's coordinate out of range) is of a similar order of magnitude as the other penalty term (Pen2) 2.
Pen2 parameters' values (C con PEN2 , C var PEN2 ) must ensure that the penalty term for suppressing solutions violating constraint C4 (no pollution of EWs) is proportionally imposed on infeasible solutions, sorting infeasible solutions fairly.

3.
Parameters that affect the relevant weightings of Pen1 and Pen2 in the total Penalty function (C f1 , C f2 ), in combination with the penalty term-specific parameters (C PEN1 , C con PEN2 , C var PEN2 ) must be assigned values that prioritize the individual penalties, and thus, the constraints' importance. Coordinates out of range that are not linked to a feasible, realistic solution (solutions violating C3) are less welcome than EWs being polluted (solutions violating C4), which is unfavorable but exists as a real physical scenario.

4.
Parameters affecting the relevant weightings of Pen1 and Pen2 in the total Penalty function (C f1 , C f2 ), in combination with the penalty term-specific parameters (C PEN1 , C con PEN2 , C var PEN2 ) must be assigned values that also favor "bad" solutions (solutions with high Fitness Value, meaning expensive groundwater management resources' solutions) over infeasible solutions, that exhibit violations of any constraint and entail a penalty.

5.
Overall, the minimum penalty rule must apply; the fittest penalty function is the lower one that can consistently and as quickly as possible, lead to penalty-free optimal solutions (e.g., [41]).

Results and Discussion
The software application created and used, "Modflow-GA_Nitro" is written in Python. It exploits the "Flopy" package for Modflow handling. It is part of the "OptiManage" v4 suite for optimal groundwater resources management. Three scenarios are investigated, featuring different values (Table 2), in respect to high or low (hi/lo) (a) energy cost C kWh (hi = 0.28244 €/kWh; lo = 0.06641 €/kWh), (b) N root uptake N uptake (hi = 0.70; lo = 0.07), and (c) N mass cost N Cost (hi = 30 €/kg; lo = 3 €/kg). Scenario S1 can be seen as a high energy cost, high N uptake, and low N cost problem version (hi-hi-lo). It can also be seen as a high energy cost, low N uptake, and high N cost problem version (hi-lo-hi). In both cases, the product N uptake × N Cost (see Equation (14)) remains the same (=2.1). Scenario S2 is a high energy cost, high N uptake, and high N cost problem version (hihi-hi). Scenario S3 is a low energy cost, high N uptake, and high N cost problem version (lo-hi-hi). For each of the three scenarios, "Modflow-GA_Nitro" is implemented five times (runs). Each Modflow simulation needed 4-5 s. For a population of 60 chromosomes and 500 generations, the average computational time of each run was approximately 38 h (Intel Core i7 7700 @3.60 GHz; 16 GB RAM @1197 MHz). Table 2. Features of the three problem scenarios studied. S1: hi-hi-lo or hi-lo-hi energy cost-Nitrogen root uptake and Nitrogen mass cost, respectively; S2: hi-hi-hi; S3: lo-hi-hi.  (14)); 2 N Cost : N price (in VB3; Equation (14)); 3 N uptake × N Cost is a product in VB3 (Equation (14)).
3.1. Scenario S1: hi-hi-lo (S1a) or hi-lo-hi (S1b) Energy Cost-Uptake-Fertilizer Cost For Scenario S1 (hi-hi-lo or hi-lo-hi), the five runs produced 170 discrete acceptable (no penalty) solutions as far as their FV is concerned. All raw results are presented in Supplementary Materials SM4_1-SM4_5, while Supplementary Materials SM5_1 presents all results in a more readable fashion together with diagrams of FV, VB1, VB2, VB3 vs. the number of generations for each run. The 170 solutions are presented in Supplementary Materials SM6 (Sheet 1). 50 out of the 170 identified solutions exhibit up to a 10% increase in the FV value compared to the optimal solution FV S1 (1) = minFV S1 . The 50 best solutions of Scenario S1 are algebraically presented in Supplementary Materials SM7 (Sheet 1) and graphically presented in Supplementary Materials SM8_1. They are also graphically simultaneously shown in Figure 5 to survey the various concepts of proposed solutions better. The algebraically optimal solution for Scenario S1 (and S2 and S3) is presented in Table 3, together with statistics for all 50 "best" solutions. It is also graphically shown in Figure 6, with the concentration map of the last 300th day of simulation (EWs' protection period). The respective video of the predicted/simulated nitrate mass/pollution spread and pumping by AWs is provided as Supplementary Materials SM10_1.  . Algebraically optimal solution of the 50 best identified of Scenario S1 (S1-1), with the nitrate concentration map of the last simulation day (300). Radii of all wells are proportional to flow rates.
While the raw produced results include AWs' results in an enumerated fashion (Q3, Q4 meaning Qadd,1, Qadd,2), they are processed to filter possible duplicate results (identical solutions with a different enumeration of AWs). Another reason is for easier postprocessing for the identification of discrete solutions and categorization into different Figure 6. Algebraically optimal solution of the 50 best identified of Scenario S1 (S1-1), with the nitrate concentration map of the last simulation day (300). Radii of all wells are proportional to flow rates. Table 3. The algebraically optimal solution for all scenarios (S1, S2, S3), together with statistics for all "best" solutions (<10% of respective minFV) for each scenario (50, 54, 103 for S1, S2, S3, respectively). * metrics of the best solutions (50 for S1, 54 for S2, 103 for S3) of the 5 runs of each scenario (exhibiting up to 10% increase in FV vs. the optimal solution).
While the raw produced results include AWs' results in an enumerated fashion (Q 3 , Q 4 meaning Q add,1 , Q add,2 ), they are processed to filter possible duplicate results (identical solutions with a different enumeration of AWs). Another reason is for easier post-processing for the identification of discrete solutions and categorization into different strategies. Thus, the results are presented by classifying AWs as high (Q add,hi ) and low (Q add,lo ), regarding their flow rates. The high flow-rate AWs are highlighted in green, while the low flow-rate wells are highlighted in red. The radius of each well is proportional to the respective flow rate (except for really low values that were increased to be visible marginally).
This concept can be explained: a system of two AWs is positioned inside the North plume in its North-East boundary, with a total flow rate of ΣQ add ≈ 5307 m 3 /d. It is primarily hydrodynamically controlling and secondarily pumping the larger North nitropollution plume. The North plume is closer to an EW (North EW; W ex,1 ). That is why it is preferred to be directly controlled by the AWs instead of the South plume, which does not need direct control by AWs. A suitable distribution of the required flow rate by the EWs (ΣQ ex = 200 L/s or 17,280 m 3 /d) between them (Q ex,1 /Q ex,2 = 0.78) further protects the North EW. It allows lower AWs' flow rates, directing the bulk North plume towards the South EW. Finally, the South EW marginally does not pump nitrate pollution above the 1 ppm detection threshold during the 300-d study period. The high-Q AW (W add,hi ) is positioned between the North plume center and the North EW. Thus, it pumps the nitrate mass, mainly attracted by the North EW high flow rate. This is common practice in PAT remediation/pollution control strategies. W add,lo is positioned near W add,hi , to the North, with a lower but similar flow rate. Its existence and position are not random. It serves this dominantly pollution-control and less PAF management concept, as the system of the two AWs can pump and/or hydraulically control the same amount of polluted groundwater with lower hydraulic head drawdowns. Thus, less energy is consumed, and the dominant in Scenario S1 pumping cost VB1 is lowered. Moreover, even though the pipe network cost VB2 (connecting AWs to the reservoir) is trivial, the layout of the AWs (minimizing the total network length) indicates that the algorithm is sensitive to its impact in minimising the total FV.
Further studying Scenario S1 results, instead of various versions/alterations of the optimal solution, other management concepts are also discovered. For example, solution S1-5 proposes using only one AW with just an FV increase of +16,525 €/yr or +3.78% compared to S1-1. S1-5 is graphically presented in Supplementary Materials SM8_1 (page 5), while the respective video of mass transport is provided as Supplementary Materials SM10_2. Therefore, a systematic investigation of all "good" solutions is imperative. This way, different management strategies are given explicitly stated criteria that can be identified in various alternative algebraical versions. Such a detailed post-processing of all good solutions/results is presented in Section 3.4.

Scenario S2: hi-hi-hi Energy Cost-Uptake-Fertilizer Cost
For Scenario S2 (hi-hi-hi), the five runs produced 154 discrete acceptable (no penalty) solutions as far as their FV is concerned. All raw results are presented in Supplementary Materials SM4_6-SM4_10. Supplementary Materials SM5_2 presents all results in a more readable fashion together with diagrams of FV, VB1, VB2, and VB3 vs. the number of generations for each run. The 154 solutions are presented in Supplementary Materials SM6 (sheet 2). 54 out of the 154 identified solutions exhibit up to a 10% increase in the FV value compared to the optimal solution FV S2 (1) = minFV S2 . The 54 best solutions of Scenario S2 are algebraically presented in Supplementary Materials SM7 (sheet 2) and graphically in Supplementary Materials SM8_2. They are also graphically simultaneously presented in Figure 7 to survey the various concepts of proposed solutions better. The algebraically optimal solution for Scenario S2 (and S1 and S3) is presented in Table 3, together with statistics for all 54 "best" solutions. It is also graphically presented in Figure 8, with the concentration map of the last 300th day of simulation (EWs' protection period). The respective video of the predicted/simulated nitrate mass/pollution spread and pumping by AWs is provided as Supplementary Materials SM10_5. together with statistics for all 54 "best" solutions. It is also graphically presented in Figure  8, with the concentration map of the last 300th day of simulation (EWs' protection period).
The algebraically optimal solution (S2-1), proposing construction of two AWs, exhibits: FV S2 (1) = 225,895 €/yr, |FV| S2 (1) = 794,277, VB1 S2 (1) = 500,397 €/yr (63.00% of |FV| S2 (1)), VB2 S2 (1) = 9689 €/yr (1.22% of |FV| S2 (1)), and VB3 S2 (1) = -284,191 €/yr (35.78% of |FV| S2 (1)). Scenario S2 is more complex than S1. S2-1 exhibits min FV but not min VB1. It is ranked 10th based on VB1, 15th based on VB2, and 37th based on VB3 (ranking from low to high; see Supplementary Materials SM7-sheet 2. One can distinguish various versions of the algebraically optimal solution S2-1, e.g., solutions S2-2 to S2-11. These solutions generally represent the same management strategy as the three best solutions of Scenario S1 but with some differences: the AWs are located closer to the center of the North plume, exhibiting higher pumping flow rates. The variation of this strategy was expected in Scenario S2, as the equilibrium between (increased) pumping cost VB1 and the profit from pumped nitrate mass VB3 is now set less towards VB1 and more towards VB3. The proposed management strategy is similar but adjusted to the new reality. It promotes more PAF and less pollution control compared to S1. This is apparent by the total predicted pumped nitrate mass by AWs (see Table 3 "NO3 mass" column), which is approximately 13.5 t in S2-1 instead of 8.4 t in S2-1 (during the studied period of 300 d). It is also apparent in the maximum concentration in the aquifer. It is approximately 30 ppm for S2-1 (see Figure 8) instead of 45 ppm for S1-1 (see Figure 6). The system of two AWs now pumps a total of ΣQ add ≈ 6783 m 3 /d (in S2-1, instead of 4520 m 3 /d in S1-1). The distribution of the required EWs' flow rate is now in the order of Q ex,1 /Q ex,2 = 0.59 (instead of 0.76 in S1-1). The algorithm again considers even the trivial contribution of VB2 in the total FV, as proven by the positioning of the AWs relative to each other and the reservoir (minimum spanning tree).
Further studying Scenario S2 results, other management concepts are also discovered. For example, solution S2-25 positions the second AW, W add,lo in the South plume with an FV increase of +7885 €/yr or +3.5% compared to S2-1. S2-25 is graphically presented in Supplementary Materials SM8_2 (page 25). The respective video of mass transport is provided as Supplementary Materials SM10_7. The systematic investigation of all "good" solutions after a detailed post-processing of all good solutions/results is presented in Section 3.4.

Scenario S3: lo-hi-hi Energy Cost-Root Uptake-Fertilizer Cost
For Scenario S3 (lo-hi-hi), the five runs produced 221 discrete acceptable regarding FV.  Figure 9 to survey the various concepts of proposed solutions better. The algebraically optimal solution for Scenario S3 (and S1 and S2) is presented in Table 3, together with statistics for all 103 "best" solutions. It is also graphically shown in Figure 10, with the concentration map of the last 300th day of simulation (EWs' protection period). The respective video of the predicted/simulated nitrate mass/pollution spread and pumping by AWs is provided as Supplementary Materials SM10_10.
The algebraically optimal solution (S3-1), proposing the construction of two AWs, exhibits FV S3 (1) = −228,204 €/yr, which is the annual profit for all farmers in the area due to the cut-down on N-based fertilizer purchase. Moreover, |FV| S3 (1) = 623,256, VB1 S3 (1) = 183,243 €/yr (29.40% of |FV| S3 (1)), VB2 S3 (1) = 14,283 €/yr (2.29% of |FV| S3 (1)), and VB3 S3 (1) = −425,730 €/yr (68.31% of |FV| S3 (1)). Scenario S3 is more complex than 1 and 2. S3-1 exhibits minimum FV but not minimum VB1. It is ranked 41st based on VB1, 67th based on VB2, and 27th based on VB3 (ranking from low to high; see Supplementary Materials SM7-sheet 3.     There are various versions of the algebraically optimal solution S3-1, e.g., solutions S3-2 to S3-54. These solutions represent a different management strategy than that of the best solutions of Scenarios S1 and S2: (a) W add,hi is positioned even closer to the center (maximum concentration) of the North plume; (b) W add,lo is positioned inside the South plume; (c) both of them exhibit even higher flow-rates than S2. The variation of this strategy was expected in S3. The equilibrium between the (decreased) pumping cost VB1 and the high profit from pumped nitrate mass VB3 is leaning towards VB3 and not VB1 for the first time. This is a primarily PAF and, secondarily, pollution control strategy. This is portrayed in Table 3 concerning the total predicted pumped nitrate mass by AWs ("NO3 mass" column). It is approximately 20.3 t in S3-1 instead of 13.5 t and 8.4 t in S2-1 and S2-1, respectively. It is also evident in the maximum aquifer concentration, equal to 15 ppm for S3-1 (see Figure 10) instead of 30 (see Figure 8) and 45 ppm (see Figure 6) for S2-1 and S1-1, respectively. The system of two AWs now pumps a total of ΣQ add ≈ 13,248 m 3 /d (in S3-1, instead of 4520 in S1-1 and 6312 in S2-1). The required EWs' total flow rate is now evenly distributed between the two EWs (Q ex,1 /Q ex,2 = 1, (instead of 0.76 in S1-1 and 0.59 in S2-1). The algorithm again considers even the trivial contribution of VB2 in the total FV. This is proven by the positioning of the AWs relative to each other and the reservoir.
Other management concepts are also produced in Scenario S3. For example, solution S3-72 positions both AWs in the North plume with similar flow rates, with an FV increase of +22,049 €/yr or +9.7% compared to S3-1. S3-72 is graphically presented in Supplementary Materials SM8_3 (page 72). The respective video of mass transport is provided as Supplementary Materials SM10_12. The whole post-processing stage of "good" or "best" solutions' investigation and strategies' identification is featured in Section 3.4.

Systematic Investigation of Solutions
The solutions of the three scenarios identified as "good" or "best 10%" are investigated further and compared against each other so that the algorithm's robustness is scrutinized. Moreover, management strategies can be identified, solutions can be classified, and valuable conclusions can be drawn.
In an attempt to map the differences between the three scenarios and the impact each variable (AWs' coordinates and flow rates and EWs' flow rates) exhibits against each cost item (VB1, VB2, VB3) and the overall cost (FV), one can rely on the respective correlation coefficient (Cc) values. Table 4 presents the Cc values for all combinations of group A values (FV, VB1, VB2, VB3) and group B values (Q ex,1 /Q ex,2 , Q add,hi , Q add,lo , ΣQ add , VB1, VB2, VB3) for all scenarios and their respective "best" solutions. Interesting Cc values (>0.5) are highlighted in bold. Figure 11 also graphically presents the variation of critical variables/costs for the three scenarios (S1, S2, S3) in box plots (Supplementary Materials SM7-Sheet 4).
The complexity of the genuinely multi-objective problem is obvious when studying the correlation coefficient values of Table 4 and the differences between the plots in Figure 11. There are few clear connections between one decision variable and one cost item, but some interpretations and conclusions can be made. They conform with what empirically a specialist would generally expect. The distribution of flow rates of EWs (Q ex,1 /Q ex,2 ) does not univocally affect FV or any item cost, particularly in any scenario. Their role is auxiliary to the AWs. On the other hand, all VB1-VB3 correlations are significantly negatively high in all scenarios, confirming that the GA led to good solutions on a VB1-VB3 (the main and antagonizing cost items/objectives) Pareto front: increase of VB1 entails a decrease of VB3, and vice versa.
In S1, the dominance of VB1 in FV is portrayed by Cc(FV-VB1) = +0.995). W add,hi is the main reason for high VB1 so Cc(Q add,hi -VB1) = +0.645 and Cc(Q add,hi -FV) = +0.653. The Q add,lo /Q add,hi ratio is low (see Table 3), so Cc(Q add,lo -VB1) and Cc(Q add,lo -FV) are low. On the other hand, Cc(ΣQ add -FV) = +0.938, since ΣQ add 's main component is Q add,hi , and Cc(ΣQ add -VB1 = +0.941), since FV's main component is VB1. The primary use of W add,hi is pollution control and pump-and-fertilize (see locations of W add,hi in Figure 5), hence Cc(Q add,hi -VB3) = −0.545. W add,lo is mainly for pollution control (see Figure 5). Therefore, it exhibits low Cc(Q add,lo -VB3), but as Q add,hi is the main component of ΣQ add , it is Cc(ΣQ add -VB3) = −0.879. VB3 s increased importance in FV and its interconnection with VB1 are also portrayed in correlations: Cc(FV-VB3) = −0.896 and Cc(VB1-VB3) = −0.926.  The complexity of the genuinely multi-objective problem is obvious when studying the correlation coefficient values of Table 4 and the differences between the plots in Figure  11. There are few clear connections between one decision variable and one cost item, but some interpretations and conclusions can be made. They conform with what empirically a specialist would generally expect. The distribution of flow rates of EWs (Qex,1/Qex,2) does not univocally affect FV or any item cost, particularly in any scenario. Their role is auxiliary to the AWs. On the other hand, all VB1-VB3 correlations are significantly negatively high in all scenarios, confirming that the GA led to good solutions on a VB1- S2 is a more complex scenario, as VB1 is still the most important cost item (about 63% of |FV|), but VB3 s impact is significantly increased (about 35% of |FV|). The complexity of S2 is portrayed in low Cc(FV-VB1) and Cc(FV-VB3) values, and also low Cc(ΣQ add -FV) values, despite Cc(ΣQ add -VB1) = +0.936 and Cc(ΣQ add -VB3) = −0.967. This can be explained by the ambiguous use of W add,lo in S2, used for pollution control or/and PAF, depending on the strategy (see Figure 7). This can also explain the odd find of Cc(Q add,lo -FV) = −0.725, while correlations of Q add, lo with VB1 and VB3 are low. W add,hi is still the main reason for high VB1 (low Q add,lo /Q add,hi values; see Table 3), hence Cc(Q add,hi -VB1) = +0.661 and Cc(Q add,hi -FV) = +0.623. Some interesting correlations regarding VB2 are due to some solutions with ΣQ add > Q critical (=50 L/s), which increases the pipe cost per m and VB2.
S3 is also complex, like an inverse S2: VB3 is now the most important cost item (about 68% of |FV|), while VB1 s impact is not trivial but decreased (about 30% of |FV|), all graphically presented in Figure 11. The complexity of S3 is portrayed in low Cc(FV-VB1) while Cc(FV-VB3) = +0.556, and also low Cc(ΣQ add -FV) and values, despite Cc(ΣQ add -VB1) = +0.988 and Cc(ΣQ add -VB3) = −0.776. This can be explained by the ambiguous use of W add,lo in S3 for pollution control or/and pump-and-fertilize, depending on the strategy (see Figure 9). This can also explain the odd find of Cc(Q add,lo -FV) = +0.755, while Cc(Q add,lo -VB1) is low and Cc(Q add,lo -VB3) = 0.592. W add,hi is still the main reason for high VB1 so (low Q add,lo /Q add,hi values; see Table 3), hence Cc(Q add,hi -VB1) = +0.823, but that does not entail high Cc(Q add,hi -FV), as VB1 does not govern FV values in S3. It does correlate extremely well with VB3, though (Cc = −0.928). The nitrate mass it pumps in any solution seriously increases the dominant VB3 cost. W add,lo is used for the first time in S3 mainly for PAF in most solutions, exploiting the South plume's nitro-pollution for indirect profit (VB3), hence the interesting Cc(Q add,lo -VB3) = 0.592. The interesting correlations regarding VB2 here are due to most of the solutions exhibiting long networks (W add,lo is in the South plume) while also proposing ΣQ add > Q critical (=50 L/s). This causes increased VB2.

Management Strategies
The statistical processing of "good" solutions alone cannot give insight into better aquifer management or provide general directions to good management strategies. The systematic investigation of all solutions is a required process to step-by-step:

2.
Propose the criteria for the initial classification of solutions into strategy.

3.
Identify the management strategies.

4.
Re-evaluate criteria and classification if needed (e.g., if a strategy includes only a single solution or if two are practically very similar).

5.
Finally, classify solutions into the final strategy.
The four currently selected criteria are presented in Table 5. These criteria identify 12 scenario-specific and nine overall aquifer management strategies. The strategies, their classification and comments on their differences are presented in Table 6. Table 7 shows the best versions of the 12 scenario-specific and nine overall strategies. Figure 12 graphically presents sketches of all management strategies identified for all scenarios in their algebraically best versions. Figure 13 presents nitrate concentrations' temporal progression in AWs for the algebraically best versions of the 12 scenario-specific (9 overall) strategies, following the selected criteria of

C4
The ratio of AWs' flow rates (R=Q add,hi /Q add,lo ) R>1 R≈1 --- Table 6. The 12 scenario-specific and nine overall strategies, together with their classification based on the criteria of Table 5, and comments on differences of various strategies ("1" = yes, "0" = no).  Table 7. Best versions of all 12 scenario-specific and nine overall strategies' classifications. All these solutions are presented in Figure 12 and Supplementary Materials SM9_1-SM9_12.  Figure 12. All management strategies identified for all scenarios in their algebraically best ve Titles suggest the overall strategy (9 strategies, A to I), scenario-specific strategy, and best v of each scenario-specific strategy; see Supplementary Materials SM7). Radii are loosely propo to respective flow rates.  Commenting on the various of Figure 13, all concentration-time (C-T) diagrams corresponding to the best algebraical solution of each scenario-specific strategy start with a positive concentration value. This occurs as all AWs are positioned inside plumes. All strategies propose construction/operation of one AW inside the North plume (S1b, S2d) or two AWs inside the North plume (S1a, S1c, S1d, S2a1, S2a2, S2c, S3b) or two plumes, with Wadd,hi inside the North plume and Wadd,lo inside the South plume (S2b, S3a1, S3a2). In any case, all AWs are positioned inside one of the two plumes. Hence the initial concentration since day 1 is always greater than zero. In some cases, where Wadd,lo is positioned deeper in the plumes than Wadd,hi, the initial concentration is higher for them (S1d, S2b, S2c). This is why the respective Wadd,lo and Wadd,hi C-T diagrams intersect. In other cases, the intersection is circumstantial (S2a2, S3a2, S3b). Finally, AWs positioned in the leading dispersive front (considering the general flow) of the two initially Gaussian plumes generally exhibit a concave C-T function (e.g., Wadd,hi of Strategy F, S2-41; Figure  13h). On the other hand, AWs positioned in the tail (considering the general flow) of any of the two plumes generally exhibit a convex function (e.g., Wadd,lo of Strategy F, S2-41; Figure 13h).
To graphically present the solutions regarding their objectives' trade-offs, a VB1-VB2-VB3 3D diagram is not needed as VB2 is trivial in every scenario. Figure 14 presents  Commenting on the various of Figure 13, all concentration-time (C-T) diagrams corresponding to the best algebraical solution of each scenario-specific strategy start with a positive concentration value. This occurs as all AWs are positioned inside plumes. All strategies propose construction/operation of one AW inside the North plume (S1b, S2d) or two AWs inside the North plume (S1a, S1c, S1d, S2a1, S2a2, S2c, S3b) or two plumes, with W add,hi inside the North plume and W add,lo inside the South plume (S2b, S3a1, S3a2). In any case, all AWs are positioned inside one of the two plumes. Hence the initial concentration since day 1 is always greater than zero. In some cases, where W add,lo is positioned deeper in the plumes than W add,hi , the initial concentration is higher for them (S1d, S2b, S2c). This is why the respective W add,lo and W add,hi C-T diagrams intersect. In other cases, the intersection is circumstantial (S2a2, S3a2, S3b). Finally, AWs positioned in the leading dispersive front (considering the general flow) of the two initially Gaussian plumes generally exhibit a concave C-T function (e.g., W add,hi of Strategy F, S2-41; Figure 13h). On the other hand, AWs positioned in the tail (considering the general flow) of any of the two plumes generally exhibit a convex function (e.g., W add,lo of Strategy F, S2-41; Figure 13h).
To graphically present the solutions regarding their objectives' trade-offs, a VB1-VB2-VB3 3D diagram is not needed as VB2 is trivial in every scenario. Figure 14 presents the correlation diagram of the two dominant cost items VB1-VB3 for the "best 10%" solutions for (a) all scenarios, (b) Scenario S1, (c) S2, and (d) S3. The indicated spaces that the three scenarios graphically take over in the 2D projected policy space are completely different. This helps to comprehend the impact external to the scientific problem (decision variables and objectives) factors (energy cost, N uptake, N price) have on a constrained, non-linear, stochastic, multi-objective optimization problem as the one studied. The algebraically optimal solutions are also indicated. Many solutions are actually on the Pareto front, implying a robust methodology of optimization and results post-processing. Moreover, solutions of a specific strategy are generally grouped, implying a good selection of criteria. the correlation diagram of the two dominant cost items VB1-VB3 for the "best 10%" solutions for (a) all scenarios, (b) Scenario S1, (c) S2, and (d) S3. The indicated spaces that the three scenarios graphically take over in the 2D projected policy space are completely different. This helps to comprehend the impact external to the scientific problem (decision variables and objectives) factors (energy cost, N uptake, N price) have on a constrained, non-linear, stochastic, multi-objective optimization problem as the one studied. The algebraically optimal solutions are also indicated. Many solutions are actually on the Pareto front, implying a robust methodology of optimization and results post-processing. Moreover, solutions of a specific strategy are generally grouped, implying a good selection of criteria. Figure 14. Correlation diagrams of the two dominant cost items VB1-VB3 for the "best 10%" solutions for scenarios: (a) S1+S2+S3, (b) S1, (c) S2, (d) S3; each strategy is differently colored; optimal solutions are indicated. Many solutions are actually on the Pareto front, while solutions of a specific strategy are generally grouped, implying a successful optimization process and a good selection of criteria.

Conclusions and Future Research
A complex theoretical optimization problem is studied: managing a confined aquifer polluted by nitrate. Nitro-pollution may be caused by intensive agricultural production and fertilizer application, sewage or industrial disposal or artificial recharge via boreholes/recharge wells using treated wastewater. This pilot theoretical approach focuses mainly on the optimization methodology and the production of alternative management strategies rather than the detailed nitrate transport phenomenon and fate per-se. This framework used a simple advection-dispersion model for the mass transport simulation in groundwater. Following other similar research work assumptions and also meeting the demands dictated by the computationally greedy nature of the optimization method, a simplified 3D hydraulic model was formulated.
Moreover, the most complex and challenging optimization method scenario concerning the initial concentration conditions was selected after many tests. That is the Figure 14. Correlation diagrams of the two dominant cost items VB1-VB3 for the "best 10%" solutions for scenarios: (a) S1+S2+S3, (b) S1, (c) S2, (d) S3; each strategy is differently colored; optimal solutions are indicated. Many solutions are actually on the Pareto front, while solutions of a specific strategy are generally grouped, implying a successful optimization process and a good selection of criteria.

Conclusions and Future Research
A complex theoretical optimization problem is studied: managing a confined aquifer polluted by nitrate. Nitro-pollution may be caused by intensive agricultural production and fertilizer application, sewage or industrial disposal or artificial recharge via boreholes/recharge wells using treated wastewater. This pilot theoretical approach focuses mainly on the optimization methodology and the production of alternative management strategies rather than the detailed nitrate transport phenomenon and fate per-se. This framework used a simple advection-dispersion model for the mass transport simulation in groundwater. Following other similar research work assumptions and also meeting the demands dictated by the computationally greedy nature of the optimization method, a simplified 3D hydraulic model was formulated.
Moreover, the most complex and challenging optimization method scenario concerning the initial concentration conditions was selected after many tests. That is the assumption of two specifically engineered Gaussian plumes, assumedly caused by two boreholes/recharge wells injecting treated wastewater directly into the studied confined aquifer. The initial and or other management/policy-making authorities for every kg of N retrieved and reused through pump-and-fertilize strategies. Scenario S3, given the lower energy cost and high N-based fertilizer value, exhibits not only decreased costs but an actual profit in the order of 215,000 €/yr. S3 could represent a special subsidization/funding by policymakers for energy used for irrigation by water management authorities or individual farmers.
The search for just the algebraically optimal solutions is myopic in such complex multi-objective problems. Nevertheless, this can be deduced given the following: • Addressing such complex, constrained, nonlinear, stochastic, multi-objective water resources management optimization problems with metaheuristics provides no guarantees in finding the algebraically optimal solution(s).

•
The use of a computationally greedy optimization tool, like GAs, requires simplification of the conceptual/simulation models, inducing uncertainty and lowering models' accuracy.

•
Optimization results have proven highly volatile and sensitive regarding external variables (e.g., energy cost, N uptake, N price) governed by environmental, socioeconomic, and political factors (e.g., climate change, political unrest, conflicts, and poor management policies).
The real practical goal is the search for many alternative (sub)optimal solutions. This approach, implemented here, and the resulting methodology utterly provides a data-ready optimization/decision support tool for the management authorities in the design phase. A pool of good alternative pollution control and pump-and-fertilize management strategies is created in various (sub)optimal algebraical versions. Decision-makers can review different alternative strategies in numerous variations to select from new, previously unseen solutions or consider design features previously unconsidered. Moreover, even if some initial data are modified after the design phase (a-posteriori modified constraints, added constraints, budget constraints etc.), alternative strategy versions or even overall strategies that meet the updated problem parameters can be selected. As a result, the economic burden is low (certainly below 10%). There is also no need to initiate new bureaucratic procedures and re-assign new costly and time-consuming simulations and investigations.
Some examples regarding the currently studied problem could be: • In Scenario S1, constructing a pumping well proposed by the optimal solution's (S1-1) location (e.g., W add,hi ; see Figure 6) could be prohibited. The reasons could be problems in expropriating or using specific private land or false geological data of specific locations. In this situation, the managing authorities could easily select another solution with the W add,hi in a different area, like S1-4 (see Supplementary Materials SM8-Page 4). The additional cost would be only 1.28% (see Supplementary Materials SM7-Sheet 1).

•
In Scenario S1, constructing a second pumping well could not be favored. The reason could be elevated drilling costs, the strictest timetable, or difficulties in permits. Then, the decision-makers could select a one-well solution, like S1-5 (see Supplementary Materials SM8-Page 5). The added cost would only be 3.62% (see Supplementary Materials SM7-Sheet 1).

•
In Scenario S3, there could be a new constraint regarding the pipe network length. The reason could be the increased construction costs or shortages of pipes, or legal/permit issues concerning the route of the pipe in the South area of the study field. Then, the decision-makers could select another solution representing a different strategy, where no pumping occurs in the South plume, like S3-83 (see Supplementary Materials SM7-Page 83); the added cost would be 9.77% (see Supplementary Materials SM7-Sheet3).
First, Future research should include the investigation of the nitrogen root uptake (N uptake ) of various crops when irrigated with nitro-polluted groundwater in relation to nitrate concentration, irrigation type, climatic conditions, soil parameters, etc. Another goal should be the more accurate modeling of the nitrogen cycle and mass transport, including reactive transport, adsorption, and desorption. Furthermore, the denitrification process and stratification phenomenon should be considered, as well as the groundwater residence time or age and the impact on the nitrogen fate and, consequently, on the management strategies. Moreover, multi-objective optimization algorithms, like NSGA-II [42], should be tested as optimization tools to compare their efficiency with the proposed methodology in the complex aquifer management problem.