Optimization of the Groundwater Remediation Process Using a Coupled Genetic Algorithm-Finite Difference Method

: In situ chemical oxidation using permanganate as an oxidant is a remediation technique often used to treat contaminated groundwater. In this paper, groundwater ﬂow with a full hydraulic conductivity tensor and remediation process through in situ chemical oxidation are simulated. The numerical approach was veriﬁed with a physical sandbox experiment and analytical solution for 2D advection-diffusion with a ﬁrst-order decay rate constant. The numerical results were in good agreement with the results of physical sandbox model and the analytical solution. The developed model was applied to two different studies, using multi-objective genetic algorithm to optimise remediation design. In order to reach the optimised design, three objectives considering three constraints were deﬁned. The time to reach the desired concentration and remediation cost regarding the number of required oxidant sources in the optimised design was less than any arbitrary design.


Introduction
In recent decades, groundwater contamination has raised concerns of the risk to human health and ecological systems. These concerns are increased because of the growing number of contaminated sites [1,2], and sources of contamination, combined with an increased wordwide water demand from domestic, industrial and agricultural sectors. Since groundwater is a well-defined pathway for contaminant transport and transfer to receptors, the importance of groundwater management is increasing. This management approach consists of identifying and conserving uncontaminated sources, and finding and using suitable approaches to manage and remediate contaminated aquifers. In situ chemical oxidation (ISCO) is one of the most popular and practical technologies for remediation because it can be used for many site types and contaminants. ISCO is based on the injection of an oxidant, such as permanganate, into the groundwater flow. Following the injection, the oxidant concentration is high but it dramatically declines due to its reaction with natural and contaminant organic matter, necessitating further frequent injections often resulting in higher costs [3]. A possible approach to handling these issues is the use of controlled-release oxidant. Mathematical modelling of groundwater flow and contaminant transport plays an essential role in many groundwater management projects, including prediction, assessment [4] and optimization [5][6][7] because it provides useful scientific information, such as predicting the effects of hydrological changes like groundwater abstraction for industrial purpose or irrigation development on the behaviour of the aquifer, for policy makers and resource planners which improve their decision-making process. The groundwater flow is governed by a partial differential equation based on Darcy's law and the conservation of mass. The groundwater velocity or Darcy velocity, which is calculated from it, couples groundwater equation to solute transport equation. To solve these equations, there are two different approaches, namely the analytical and the numerical approaches. Different numerical approaches, such as finite difference methods (FDM) [8], finite element methods (FEM) [9][10][11][12][13][14], finite volume methods (FVM) [15] and Meshfree methods including radial point collocation method (RPCM) [16,17], and radial base functions (RBFS) [18,19] can be utilised to solve the governing equation of the flow and solute transport in groundwater. The selection of the best method among these depends on site-specific needs, availability of input data, expectation and use of the model results, and the analytical solution can be used as a benchmark to show the accuracy of the numerical approach used, considering different conditions such as the aquifer's heterogeneity [11]. Based on different numerical approaches, there are various codes for groundwater modeling, including MODFLOW [20], FEFLOW [21], HYDRUS [22], HydroGeoSystem [23], Parflow [24], each with its own capabilities, operational features, and limitations. These codes can model groundwater flows as well as unsaturated flow above the water table and stream flow. All of these codes solve the classic three-dimensional groundwater flow equation, subject to initial and boundary conditions. However, the codes differ in the algorithms used to solve the equation, they differ in their representation of the initial and boundary conditions. MODFLOW requires model layers to be continuous throughout the entire model domain, making it a challenge within a larger body to model smaller material formations, without adding a large number of unnecessary cells. Moreover, it is difficult to simulate complex geological features, such as angled faults and simulate steep hydraulic gradients. FEFLOW is similar to MODFLOW in the sense that model layers must be continuous across the entire model domain. All these codes and in house software needs to be validated. Besides analytical solutions, physical models using sandboxes not only benchmarks, but provides relevant and useful data to validate fluid flow and solute transport [25,26] but can also be used for an estimation procedure of contaminant transport parameters [27].
Mathematical modelling of groundwater and contaminant fate together with the optimization of remediation processes can offer more appropriate scenarios to tackle contaminated groundwater problems to achieve the effective and efficient design of remediation [28]. New optimisation methods such as adaptive and natural computing methods have been utilised because in many engineering areas, conventional mathematical optimisation methods, such as linear programming (LP), quadratic programming (QP) and nonlinear pro-programming (NLP) are no longer effective for considering the increasing complexity of problems. Natural Computing methods are a general term which refers to computing inspired by nature and natural systems. These methods belong to a relatively large family of analysis methods which utilizes biomimicry to translate and manipulate data. These methods include neural networks which mimics the mechanisms of the nervous system, fuzzy systems based on an extension of traditional logic in order to represent uncertainty and qualitative reasoning, machine learning approaches and general optimization techniques such as Evolutionary Computation on simulation of biological evolution, swarm intelligence based on simulation of social behaviour of animals, immunocomputing inspired by the biological immune system, and simulated annealing derived by Statistical Mechanics. Some of the natural computing methods such as artificial neural network [29][30][31], fuzzy based [32,33], ant colony algorithms [34,35], genetic algorithms (GA) [29] and particle swarm optimization [36] are often used in groundwater management issues.
The GA, which were introduced by John Henry Holland [37] based on the Darwinian evolution theory, is one of the advanced methods in artificial intelligence and heuristic search algorithms that mimic the mechanics of natural selection and natural genetics which use the principles of evolution known as genetic operators (population, selection, crossover and mutation) to discover a better solution or optimization to a particular problem. The concept of GA is based on the initial selection of a relatively small population. Each individual in the population represents a possible solution in the parameter space. The fitness of each individual is determined by the value of the objective function, calculated from a set of parameters which can be used to impose the constraints on the solution. The natural evolutional processes of reproduction, selection, crossover, and mutation are applied using probability rules to evolve new and better generations. The probabilistic rules allow some less fit individuals to survive. GA has an advantage over conventional optimization strategies because it does not require derivative knowledge to solve complex and discontinuous optimization problems.
In this paper, the coupled groundwater and reactive transport of oxidant and contaminant considering the full tensor of hydraulic and The FDM discretization scheme are introduced. The physical sandbox experiment and analytical solution of 2D advectiondiffusion were used as benchmarks to show the developed model's accuracy. The verified model is coupled with the GA approach in order to find the optimum remediation design.

Groundwater Flow
The 2D groundwater flow equation through a saturated anisotropic porous aquifer can be written as [38,39] S ∂h where S is the storage coefficient, h(x, y, t) is the piezometric head [L], t is the time, ∇ is the divergence operator, Q w is the source or sink term L 3 T −1 L −2 and q denotes the recharge rate LT −1 , K = K(x) is the hydraulic-conductivity symmetric full tensor LT −1 expressed as The following initial and boundary conditions, considering Ω as an aquifer domain Γ as its Lipschitz continuous boundary comprising ∂Ω = Γ D ⊕ Γ N , are imposed in the groundwater flow equation. Where Γ D and Γ N interpret the portions of Γ for Dirichlet and Neumann boundary conditions c.f. Figure 1, h(x, y, t) = h 1 (x, y) (x, y) ∈ Γ D . ∂h Water 2021, 1, 0 3 of 20 crossover and mutation) to discover a better solution or optimization to a particular problem. The concept of GA is based on the initial selection of a relatively small population. Each individual in the population represents a possible solution in the parameter space. The fitness of each individual is determined by the value of the objective function, calculated from a set of parameters which can be used to impose the constraints on the solution. The natural evolutional processes of reproduction, selection, crossover, and mutation are applied using probability rules to evolve new and better generations. The probabilistic rules allow some less fit individuals to survive. GA has an advantage over conventional optimization strategies because it does not require derivative knowledge to solve complex and discontinuous optimization problems. In this paper, the coupled groundwater and reactive transport of oxidant and contaminant considering the full tensor of hydraulic and The FDM discretization scheme are introduced. The physical sandbox experiment and analytical solution of 2D advectiondiffusion were used as benchmarks to show the developed model's accuracy. The verified model is coupled with the GA approach in order to find the optimum remediation design.

Groundwater Flow
The 2D groundwater flow equation through a saturated anisotropic porous aquifer can be written as [38,39] where S is the storage coefficient, h(x, y, t) is the piezometric head [L], t is the time, ∇ is the divergence operator, Q w is the source or sink term L 3 T −1 L −2 and q denotes the recharge rate LT −1 , K = K(x) is the hydraulic-conductivity symmetric full tensor LT −1 expressed as The following initial and boundary conditions, considering Ω as an aquifer domain Γ as its Lipschitz continuous boundary comprising ∂Ω = Γ D ⊕ Γ N , are imposed in the groundwater flow equation. Where Γ D and Γ N interpret the portions of Γ for Dirichlet and Neumann boundary conditions c.f. Figure 1,

Reactive Transport
Reactive transport of the contaminant and oxidant in groundwater is given by the following coupled advection-dispersion equations [40,41] R where R = 1 + ae b k d /θ is the retardation factor and and considered to be equal to 2. It describes sorption and considered the same for the contaminant and the oxidant, in which ae b is the soil bulk density, k d is the sorption coefficient, and θ denotes the porosity of the aquifer, t denotes time, D is the dispersion L 2 T −1 , C 1 and C 2 are concentrations of organic matter as a contaminant and permanganate as an oxidant respectively ML −3 , k sor is the second-order reaction constant T −1 and F Release is the release function of permangante. v is the velocity vector LT −1 . The components of the velocity vector, namely v x and v y , are evaluated from the piezometric head equations using the following relations [38,39] v where K xx and K yy are the hydraulic conductivities in x and y directions respectively. The components of the dispersion coefficient tensor, D = D(x), are evaluated by the following relations: where α L and α T are longitudinal and transverse dispersivity and D * is the effective molecular diffusion coefficient. v x and v y in the Equations (7) and (8) are evaluated from the flow equation, and these two equations couple the groundwater flow and reactive transport. In our study, the following initial and boundary conditions are imposed on the reactive transport equations

FDM Scheme
The FDM is a suitable and well-known numerical approach to solve coupled differential equations of groundwater flow and reactive transports of contaminant and oxidant, in which differential quotients are replaced by difference quotients.

Time Discretization
Considering time span [0, T] is divided into n t , a time period with ∆t = T n t . The temporal discrete of Equations (1), (5) and (6) based on the implicit FDM, are given by where Φ denotes the piezometric head and concentration of contaminant and oxidant in Equations (1), (5) and (6), and κ indicates the prefactors in these equations.

Spatial Discretisation
The aquifer domain [0, L x ] × 0, L y with an N x + 1 × N y + 1-point grid, defined by where ∆x and ∆y are the grid spacings in each direction. The gradients of the piezometric head and the concentration of contaminant and oxidant at the interface of cells are determined using the following equations ∂Φ Herein again, Φ denotes the unkown field parameter, piezometric head, concenteration of contaminant and oxidant.
By substituing Equations (14)- (16) in Equation (1) and multiplying both sides of the equation by ∆t and ∆x∆y, the following FDM scheme is found for the piezometric head equation Appendix A presents A 1 i,j A 9 i,j which is used in this scheme and the same approach was followed for the reactive transport equation of contaminant and oxidant.

Model Verification
Although numerical simulation models are valuable tools for analyzing groundwater systems, their predictive accuracy is limited. Increased sophistication of groundwater models has created a discrepancy between model predictions and the ability to verify or develop trust in predictions. Since we developed a model based on the introduced FD formulations in MATLAB to study contaminant fate in an anisotropic porous aquifer, we used the physical sandbox model and the analytical solution to reduce the uncertainty of the model and narrow the possible outcomes range.

Numerical Model Verification with the Physical Sandbox Model
In this study, sandbox experiments were performed to validate our simulation for mapping a changing groundwater plume distribution over time. A plexiglass tank (dimensions: 1.50 × 0.38 × 0.10 [m]) was constructed to conduct these experiments. The ends of the tanks consist of constant head tanks which are separated from the rest of the box by one impermeable wall and one perforated steel mesh filter to separate the sand from the head tanks. A peristaltic pump (Watson Marlow) was used to circulate water through the system. The pump is capable of delivering a maximum of 42 [L/h]. The characteristics of sand used are given in Table 1. Hydraulic conductivity was measured by using the constant-head method [42,43]. To mitigate the creation of preferential pathways and air bubbles, the tank was filled with a layer of a few centimeters of dry sand, after which tap water was added to saturate and cover the sand. More dry sand was layered over this now saturated sand, and the added sand also covered with tap water. This process was repeated until the tank was full. A 0.5% [w/v] potassium permanganate solution was selected and used for the changing plume distribution simulation over time. Because it is a commonly used oxidant for remediation and because of its intense purple color we were able to gain visual observations of the plume's development. The potassium permanganate was injected into the tank through the end of a 3 mm diameter tube located in the centre of the perforated wall of the tank as illustrated in Figure 2, using specialist tubing that helped control the flow from peristaltic pumps. The injection rate was 16 [L/h]. 5.5 liters permanganate were injected from the point source and after completing the injection, the experiment continued until the whole permanganate solution was diluted and washed from the sand by the water flow.

Model Verification with the Analytical Solution
To verify two-dimensional transport problems, the following combined advectiondiffusion with the first-order decay rate constant equation is considered as follows where C is the contaminant concentration ML 3 , ρ b is the Bulk density of soil ML 3 , θ is soil porosity, S is Sorbed concentration MM −1 , D xx and D yy are diffusion coefficients in x and y directions L 2 T −1 , λ is the first-order decay rate constant T −1 , α is the firstorder desorption rate constant T −1 and k d is the sorption distribution coefficient and is dimensionless. The initial and boundary conditions c.f. Figure 3 are Water 2021, 1, 0 7 of 20

Model Verification with the Analytical Solution
To verify two-dimensional transport problems, the following combined advectiondiffusion with the first-order decay rate constant equation is considered as follows where C is the contaminant concentration ML 3 , ρ b is the Bulk density of soil ML 3 , θ is soil porosity, S is Sorbed concentration MM −1 , D xx and D yy are diffusion coefficients in x and y directions L 2 T −1 , λ is the first-order decay rate constant T −1 , α is the firstorder desorption rate constant T −1 and k d is the sorption distribution coefficient and is dimensionless.
The initial and boundary conditions c.f. Figure 3 are The analytical solution to this problem in the Laplace domain was given by [44] as follows where x and y are the coordinates in x and y directions respectively, s is the Laplace complex variable, C 0 is the concentration at the contaminant source, b is the aquifer width andĈ(n) is presented in Appendix B. The physical constants and parameters for the corresponding analytical solution are are summarized Table 2. The analytical solution to this problem in the Laplace domain was given by [44] as follows where x and y are the coordinates in x and y directions respectively, s is the Laplace complex variable, C 0 is the concentration at the contaminant source, b is the aquifer width andĈ(n) is presented in Appendix B. The physical constants and parameters for the corresponding analytical solution are are summarized Table 2.

Remediation Design Optimization Using the GA Approach
A multiobjective optimisation algorithm is applied to the coupled groundwater flow and reactive transport problem considering ISCO remediation by slow release of permanganate. This problem involves finding the set of optimal remediation design regarding the y coordinate of oxidant sources whitin a contaminant plume. The GA is applied to simultaneously maximize the effectiveness of the remediation by minimising the average concentration of the contaminant in the geometry at a shorter time and maximise the region where the contaminant concentration is less than or equal to the final desired concentration of the contaminant and to minimise the cost of the remediation process by decreasing the number of necessary oxidant sources. To achieve this goal, we have defined the following functions: is the average concentration of the contaminant, N is the number of nodes, C 1 is the concentration of the contaminant at each node, C 1 FDCC is the final desired concentration of the contaminant, TC is the total cost of the remediation process, n OS is the number of oxidant sources and COS denotes the cost of each oxidant source. To achieve the objective of the study, we wish to minimise GA1 and GA3 and maximise GA2 min GA1 max GA2 min GA3 (23) by considering the following constraints: 1. The coordinates of oxidants in the X direction are fixed x| oxidant sources = x fixed (24) 2. The vertical distance between oxidant sources where d is |y ox1 − y ox2 | and d c is the critical distance between oxidant sources. The critical distance is the vertical distance between two oxidant sources whose influence domains overlap more than 75%. The influence domain is defined as a region in where the contaminant concentration is less than 85% of its initial concentration after 50 days if it implements alone. The following function defines the influence domain which is used to define critical distance where C 1 0 is initial contaminant concentration. 3. The number of oxidant sources where n OS is the actual number of oxidant sources. Figure 4 compares the results of the sandbox test and the numerical approach for the different time. To compare the sandbox test and the numerical approach, we have considered three sampling points at the tank, located at (36, 18.5) [cm], (56, 18.5) [cm], and (94, 18.5) [cm] with respect to the origin which is located at the bottom corner of the input end of the tank, and the permanganate concentration was measured every 5 min after the injection. The percentage difference between the measured and predicted concentration shows that the results are in good agreement cf. Table 3. During the experimental procedure, the concentrations recorded were similar between the different observation points. There was a peak concentration of permanganate in observation point no. 3 after 30 min compared with the other observation points and initial concentrations. This peak occurred by the movement of the mass and the compression in specific points, and the dilution effect caused by the hydrodynamic dispersivity is smaller on this scale than in real conditions.

Model Verification with the Analytical Solution
The contaminant of concentration C = 200[mg/L] is injected at x = 0 4 y 6 requires to find its concentration at the downstream for 100 minutes. Figure 5 compa FD and analytical solutions for two-dimensional transport from a continuous line so a confined aquifer. The simulation was done for t = 100 minutes.

Model Verification with the Analytical Solution
The contaminant of concentration C = 200 [mg/L] is injected at x = 0 4 y 6, and it requires to find its concentration at the downstream for 100 min. Figure 5 compares the FD and analytical solutions for two-dimensional transport from a continuous line source in a confined aquifer. The simulation was done for t = 100 min.

Model Verification with the Analytical Solution
The contaminant of concentration C = 200 [mg/L] is injected at x = 0 4 y 6, and it requires to find its concentration at the downstream for 100 min. Figure 5 compares the FD and analytical solutions for two-dimensional transport from a continuous line source in a confined aquifer. The simulation was done for t = 100 min.   Table 4 shows the RMSE error measures from Figure 6a,b. It can be seen that there is a very good agreement between the FD concentration values and the analytical results.   Table 4 shows the RMSE error measures from Figure 6a,b. It can be seen that there is a very good agreement between the FD concentration values and the analytical results.   Table 4 shows the RMSE error measures from Figure 6a,b. It can be seen that there i very good agreement between the FD concentration values and the analytical results.    (   Table 4 shows the RMSE error measures from Figure 6a,b. It can be seen that there very good agreement between the FD concentration values and the analytical results.

Optimized Remediation Design
The aquifer domain is made of rectangles of 100 m by 20 m which are heterogeneous regarding hydraulic conductivity. The reliability of groundwater modelling and contaminant transport depends on the quality of the characterisation of hydraulic conductivities.

Optimized Remediation Design
The aquifer domain is made of rectangles of 100 m by 20 m which are heterogeneous regarding hydraulic conductivity. The reliability of groundwater modelling and contaminant transport depends on the quality of the characterisation of hydraulic conductivities. The sequential Gaussian simulation can be used to model conditional stochastic continuous variables such as hydraulic conductivity of an aquifer [45]. The Sequential Gaussian begins with defining the univariate distribution of values, performing a normal score transform of the original values to a standard normal distribution. Figure 7 demonstrates the distribution of hydraulic implemented in the simulation. We have carried out two different studies, in both of which two oxidant sources have been considered to remediate the contaminant. The second order reaction rate controls the rate of contaminant oxidation by permanganate. The GA approach is used to find the optimum location of oxidant sources with respect to the discussed criteria. The functions presented by [19,46] were modified to simulate the oxidant release.

Study 1
In this study, two oxidant sources have been considered, whose x coordinates are x = 20 m and the optimum y coordinates have resulted from GA. The initial contaminant concentration is C 1 t0 = 10[mg/L] in the whole geometry. Figure 8a demonstrates the iso piezometric head contours. The optimum Y coordinate of the oxidant sources was 7.9 [m] and 12.4 [m]. Figure 8b,c show the contaminant and oxidant concentration at different times. The simulations were performed for 200 days and the time step is 0.001 day. Figure 9 illustrates the contaminant concentration at three different observation points located at the (36,10) [m], (60,15) [m] and (75, 10) [m] at the downside of the stream. As can be seen in Figure 8a, the piezometric head decreases from left to right overall. Contours of piezometric heads have peaks between the oxidant sources, as expected. As time passes, the piezometric head increases in value and its contours have higher peaks between the oxidant sources. Also, from Figure 8b

Study 2
In the second study a pumping borehole which is located at (75,5) m is considered for piezometric head which causes the contaminant and oxidant to flow to the bottom-left of the geometry. Figure 10a

Study 2
In the second study a pumping borehole which is located at (75,5) m is considered for piezometric head which causes the contaminant and oxidant to flow to the bottom-left of the geometry. Figure 10a

Discussion
The purpose of our numerical study is to better understand the remediation process of groundwater through the release of permanganate considering full tensor of hydraulic conductivity which resembles the real condition of the aquifer. To find the critical distance between two sources, we carried out prior simulations to decrease the simulation time. We found that the critical distance in our simulation is d c = 3.6 [m]. The results of prior simulations revealed that if the vertical distance between two sources is less than the critical distance, in some cases the reaching time to the desire average concentration in the geometry is more than implementing only one oxidant source located at (20,5) [m], and it can translate as reducing the performance of the remediation design. The effect of considering the full tensor of hydraulic conductivity in the flow equation is evident in the iso-head contours, and it is especially important for contaminant transport predictions because the velocity Equation (8) is affected by this . Overall, the water head contours decreased from left to right, and it mounded between two oxidant sources as expected. As we expected, the concentration of oxidant closer to the sources behave similar to the sources, but with increased distance from sources it changes rapidly. Farther away from sources, e.g. three times more than the distance between oxidant sources, there is less impact on the contaminant concentration. In the first study after 65 days, the contaminant concentration reached the 70 % of initial contaminant concentration, while in the second study after 68 days. From Figure 9 it is obvious that the contaminant concentration in optimized remediation design is much less than arbitrary design, which is randomly generated. In Study 1 the maximum difference between contaminant concentration at observation points located at (36,10) m, (60,15) m and (75,10) are 32.8%, 33.7% and 41.2% respectively and it occures at day 7, 14 and 36 after remediation, but after that the difference decreases reaching 14%, 8% and 6.5% at first, second and third observation point respectively. In other

Discussion
The purpose of our numerical study is to better understand the remediation process of groundwater through the release of permanganate considering full tensor of hydraulic conductivity which resembles the real condition of the aquifer. To find the critical distance between two sources, we carried out prior simulations to decrease the simulation time. We found that the critical distance in our simulation is d c = 3.6 [m]. The results of prior simulations revealed that if the vertical distance between two sources is less than the critical distance, in some cases the reaching time to the desire average concentration in the geometry is more than implementing only one oxidant source located at (20,5) [m], and it can translate as reducing the performance of the remediation design. The effect of considering the full tensor of hydraulic conductivity in the flow equation is evident in the iso-head contours, and it is especially important for contaminant transport predictions because the velocity Equation (8) is affected by this. Overall, the water head contours decreased from left to right, and it mounded between two oxidant sources as expected. As we expected, the concentration of oxidant closer to the sources behave similar to the sources, but with increased distance from sources it changes rapidly. Farther away from sources, e.g. three times more than the distance between oxidant sources, there is less impact on the contaminant concentration. In the first study after 65 days, the contaminant concentration reached the 70% of initial contaminant concentration, while in the second study after 68 days. From Figure 9 it is obvious that the contaminant concentration in optimized remediation design is much less than arbitrary design, which is randomly generated. In Study 1 the maximum difference between contaminant concentration at observation points located at (36,10) m, (60,15) m and (75,10) are 32.8%, 33.7% and 41.2% respectively and it occures at day 7, 14 and 36 after remediation, but after that the difference decreases reaching 14%, 8% and 6.5% at first, second and third observation point respectively. In other words, to reach the same contaminant concentration in the arbitrary design, either initial permanganate concentration at the sources must be higher, or the number of sources has to increase leading to higher cost. It can be concluded that the optimised process improves oxidant concentration and influence farther from the sources early on in treatment, though it has greater impact in the region close to the sources in later stages when compared to arbitrary design.

Conclusions
In this study, the simulation-optimisation model of groundwater and reactive transport of contaminant and slow release of permanganate, using the modification of the presented function by [19,46] for our study, in an anisotropic porous aquifer have been developed. The groundwater flow has been simulated considering heterogeneity regarding hydraulic conductivity using its full tensor leading to the simulation resembles a real aquifer. The developed model was verified using the sandbox experiment and analytical solution and has been applied to two different studies. In both studies, the optimum Y coordinates of oxidant sources using GA have been found. To find these optimum locations, three criteria were considered, including minimising the average concentration of contaminant and remediation costs in terms of the number of remediation sources and maximising the influence domain of the oxidant. The results revealed that optimisation reduces not only the time to reach the desired concentration but also the number of oxidant sources necessary, hence reducing the overall remediation cost. In addition, in the beginning of treatment with the optimised design, areas farther from oxidant sources were greater positively with higher oxidant concentrations, later on, however, the regions closer to the oxidant sources were better affected when compared to those treated with the arbitrary design. We will further extend the proposed model to deal with 3D models and will also find the optimum design of oxidant sources in the complex geometries. The next step will be to compare the effectiveness of the FD, FEM and Meshless methods to simulate the reactive transport of solutes in groundwater. Funding: This work was completed as part of the REMEDIATE (Improved decision-making in contaminated land site investigation and risk assessment) Marie-Curie Innovation Training Network. The network has received funding from the European Union's Horizon 2020 Programme for research, technological development and demonstration under grant agreement n. 643087. REMEDIATE is coordinated by the QUESTOR Centre at Queen's University Belfast http://questor.qub.ac.uk/ REMEDIATE/. In addition, this work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation). Project numbers: -327154368 (SFB 1313) -390740016 (Germany's Excellence Strategy EXC 2075/1).

Conflicts of Interest:
The authors declare no conflict of interest.