A Hybrid Neural Network–Particle Swarm Optimization Informed Spatial Interpolation Technique for Groundwater Quality Mapping in a Small Island Province of the Philippines

Water quality monitoring demands the use of spatial interpolation techniques due to on-ground challenges. The implementation of various spatial interpolation methods results in significant variations from the true spatial distribution of water quality in a specific location. The aim of this research is to improve mapping prediction capabilities of spatial interpolation algorithms by using a neural network with the particle swarm optimization (NN-PSO) technique. Hybrid interpolation approaches were evaluated and compared by cross-validation using mean absolute error (MAE) and Pearson’s correlation coefficient (R). The governing interpolation techniques for the physicochemical parameters of groundwater (GW) and heavy metal concentrations were the geostatistical approaches combined with NN-PSO. The best methods for physicochemical characteristics and heavy metal concentrations were observed to have the least MAE and R values, ranging from 1.7 to 4.3 times and 1.2 to 5.6 times higher than the interpolation technique without the NN-PSO for the dry and wet season, respectively. The hybrid interpolation methods exhibit an improved performance as compared to the non-hybrid methods. The application of NN-PSO technique to spatial interpolation methods was found to be a promising approach for improving the accuracy of spatial maps for GW quality.


Introduction
Acid mine drainage (AMD) is a natural or man-made environmental occurrence that transpires when sulfide minerals are exposed to weathering conditions or as a result of specific mining activities. It is a discharge with low pH, high heavy metal, and deadly component concentrations. It is also typically formed when sulfide-abundant wastes have been introduced to the environment. This condition has been viewed as a severe environmental issue encountered by mineral extraction enterprises around the world [1,2]. AMD is caused by oxidation of pyrite and other sulfate metals when mining sources are exposed to air, microbial activity, and water. Among the potentially hazardous dissolved metals found in high quantities, iron (II) is the most prevalent and frequent in the majority of AMD locations. Iron (II) in AMD interacts with dissolved oxygen to form iron oxide precipitates often referred to as yellow boy and may kill life all along river or stream banks [3,4]. This problem has been identified as a significant environmental concern for

Materials and Methods
The study uses NN-PSO methodology relative to a range of interpolation types, including deterministic techniques, geostatistical methods, and interpolation with barriers, to improve its GW quality prediction performance. The subsequent sections detail the study's focus and the combined strategies used to create enhanced spatial maps.

The Area of Study
The Province of Marinduque, the smallest island province in the MIMAROPA Region (or Region IV-B), is the area of research. Marinduque Island, about 200 km south of Manila, is a province known as Philippines' heart due to its geometric shape and geographical location [55]. With a total land area of 96,000 hectares, Marinduque is a 4th income class island province made up of six municipalities: Boac, Buenavista, Mogpog, Gasan, Santa Cruz, and Torrijos. The province's topography is mostly mountainous, with continuous and severe slope areas. Considering the province's overall land area, 77 percent, or 737.2 square kilometers, is classified as alienable land, while the remaining 23 percent, or 222.05 square kilometers, is classified as forest land [56]. The island province has a Climate Type III climate, with the dry season lasting from November to April and the wet season covering Toxics 2021, 9,273 4 of 33 the rest of the year [57]. The annual rainfall in the province of Marinduque range from 1700 to 2500 mm [58].
Several bodies of water, primarily rivers and their tributaries, make up the province of Marinduque. The province's major rivers and tributaries have a total length of 178 km which is composed of the municipalities of Boac (20.11 percent), Buenavista (10.96 percent), Gasan (19.49 percent), Mogpog (20.06 percent), Santa Cruz (15.22 percent), and Torrijos (15.22 percent) (14.16 percent). The province has a total of 614.1003 km 2 considering the drainage area. The Municipality of Boac has 34.80 percent of the total surface water body drainage area, followed by 7.47 percent in the Municipality of Buenavista, 11.30 percent in the Municipality of Gasan, 14.61 percent in the Municipality of Mogpog, 19.58 percent in the Municipality of Santa Cruz, and 12.24 percent in Torrijos [59].  Table 1.
continuous and severe slope areas. Considering the province's overall land area, 77 percent, or 737.2 square kilometers, is classified as alienable land, while the remaining 23 percent, or 222.05 square kilometers, is classified as forest land [56]. The island province has a Climate Type III climate, with the dry season lasting from November to April and the wet season covering the rest of the year [57]. The annual rainfall in the province of Marinduque range from 1700 to 2500 mm [58].
Several bodies of water, primarily rivers and their tributaries, make up the province of Marinduque. The province's major rivers and tributaries have a total length of 178 kilometers which is composed of the municipalities of Boac (20.11 percent), Buenavista (10.96 percent), Gasan (19.49 percent), Mogpog (20.06 percent), Santa Cruz (15.22 percent), and Torrijos (15.22 percent) (14.16 percent). The province has a total of 614.1003 km 2 considering the drainage area. The Municipality of Boac has 34.80 percent of the total surface water body drainage area, followed by 7.47 percent in the Municipality of Buenavista, 11.30 percent in the Municipality of Gasan, 14.61 percent in the Municipality of Mogpog, 19.58 percent in the Municipality of Santa Cruz, and 12.24 percent in Torrijos [59]. Figure 1 depicts 34 watersheds that make up the province of Marinduque. With an extent of 195.94 km 2 , the Boac Watershed is the province's largest watershed. The list of watersheds of the province of Marinduque is shown in Table 1.
Marinduque Island is home to one of the Philippines' largest copper deposits. Since 1969, copper mining activities have been carried out on the island. Mine tailings from these activities began to be deposited in Calancan Bay in 1975. Since then and until 1997, around 200 million tons have been dumped [60].   Marinduque Island is home to one of the Philippines' largest copper deposits. Since 1969, copper mining activities have been carried out on the island. Mine tailings from these activities began to be deposited in Calancan Bay in 1975. Since then and until 1997, around 200 million tons have been dumped [60].
In the 1990s, the island was hit by two mining catastrophes. The first incident occurred in 1993, when the Maguilaguila siltation dam in San Antonio, Sta. Cruz, collapsed. This caused property and agricultural damage and adverse effects to public health in downstream communities [61]. Three years later, in 1996, the Tapian Pit collapsed, releasing between 180,000 and 260,000 cubic meters of mine tailings into the Boac River, causing environmental and community damages [62]. GW samples were gathered from different wells in six municipalities in Marinduque and stored in plastic bottles (1 L). Data were collected and stored in compliance with EPA No. SESDPROC-301-R3, which is the GW sampling operational procedure [63]. The Hanna HI 9811-5 handheld multi-parameter sampler was used to collect all field measurements, including in situ physicochemical parameters for GW samples, such as temperature (in Celsius), pH, total dissolved solids (TDS) in milligrams per liter, and electrical conductivity (EC) in microsiemens per centimeter. Two separate geographical maps were created to depict dry and wet seasons in a year [64]. The map and details of the sampling locations are presented in Figure 2 and Table 2.
Toxics 2021, 9, x FOR PEER REVIEW  5 of 33   9  Murallon-Boac  26  Torrijos  10  Ihatub-Boac  27  Tambangan-Santa Cruz  11 Caganhao-Boac 28 Tawiran-Tagum  12 Maybo-Boac 29 Tagum In the 1990s, the island was hit by two mining catastrophes. The first incident occurred in 1993, when the Maguilaguila siltation dam in San Antonio, Sta. Cruz, collapsed. This caused property and agricultural damage and adverse effects to public health in downstream communities [61]. Three years later, in 1996, the Tapian Pit collapsed, releasing between 180,000 and 260,000 cubic meters of mine tailings into the Boac River, causing environmental and community damages [62].

Sampling, Storage, and Collection of GW Samples
GW samples were gathered from different wells in six municipalities in Marinduque and stored in plastic bottles (1 L). Data were collected and stored in compliance with EPA No. SESDPROC-301-R3, which is the GW sampling operational procedure [63]. The Hanna HI 9811-5 handheld multi-parameter sampler was used to collect all field measurements, including in situ physicochemical parameters for GW samples, such as temperature (in Celsius), pH, total dissolved solids (TDS) in milligrams per liter, and electrical conductivity (EC) in microsiemens per centimeter. Two separate geographical maps were created to depict dry and wet seasons in a year [64]. The map and details of the sampling locations are presented in Figure 2 and Table 2.

Elemental Analysis of Groundwater Samples
The measurement of total metals, which includes suspended and dissolved components as well as soluble metals, is required when looking for HMs in GW samples. The EPA Method 3005A was employed for Inductively Coupled Plasma spectroscopy using water acid digestion for total dissolved and recoverable metals as the reference guidelines for GW sample digestion [65,66].

Descriptive and Multivariate Statistical Analysis
The IBM Statistical Package for the Social Sciences (SPSS) was utilized to evaluate descriptive statistics linked to GW physicochemical parameters and HM intensities. Skewness and kurtosis were applied to assess the asymmetry of physicochemical characteristics and HM concentrations in GW. The skewness of the GW quality parameters shows the Toxics 2021, 9, 273 7 of 33 relative locations of the median and mean, whereas the kurtosis represents the form of the distribution [67,68]. The most important element in characterizing the variability of GW physicochemical parameters and HM content was the coefficient of variation (CV). The coefficient of variability was used to analyze the dataset's variability as follows: CV ≤ 15%, low; 15% < CV ≤ 35%, intermediate; and CV ≥ 35%, high [69]. Moreover, a Kolmogorov-Smirnov (K-S) test was employed to test the normality of the datasets and to examine if the GW quality parameters have a normal distribution [70].
The GW quality data are frequently identified and evaluated using multivariate statistical analysis. Multivariate statistical approaches enable the extraction of meaningful meaning from data by simplifying, organizing, and classifying it [71]. Using a correlation matrix utilizing MATLAB 2021a and R studio, the relationship between physicochemical parameters and HM intensities in GW in the research area was observed. The correlation matrix established occurrence, HM associations, and potential source of contaminants in the area of study [72]. The R value of −1 signifies that the parameter shifts inversely with respect to the other. A very strong correlation is exhibited by 0.90 < r < 1.00, strong correlation by 0.70 < r < 0.89, moderate correlation by 0.40 < r < 0.69, weak correlation by 0.10 < r < 0.39, and negligible correlation by 0 < r < 0.10. There is no association between the two variables if the correlation is zero [73].

Machine Learning: Hybrid Neuro-Particle Swarm Optimization Modelling
Machine learning is an application of artificial intelligence that allows software applications and uses statistical models and algorithms to analyze and draw interferences from patterns in data. This study uses MATLAB 2021a to enhance the prediction capability of the spatial interpolation maps of GW quality using a Particle Swarm Optimization (PSO) trained Artificial Neural Network model. Subsequent sections below elaborate how machine learning was used as technique in creating GW quality mapping. Hence, pages 7-9 discuss and illustrate how machine learning has been used in the study and which stage of the technique development been used.

Backpropagation Neural Network (BP-NN)
The ANN is an approach inspired by a real biological neuron that has been used in prediction and forecasting, particularly for complicated and non-linear systems such as environmental problems such as water quality modeling [74]. The Artificial Neural Network learns by training the connectivity between the neurons, which is performed using known input and output values provided in an organized manner so that the network can extract the relationship and patterns in the dataset [75]. MATLAB R2021a was used to create the neural network model, which included 70 percent, 15 percent, and 15 percent data partitioning for the data sets utilized in the training, validation, and testing stage, respectively [76]. The Levenberg-Marquardt algorithm was employed as the model's training algorithm since it is the quickest method to train a moderate-sized feed forward neural network with several hundred weights [77]. The model uses a hyperbolic tangent sigmoid (tansig) transfer function as the driving component for the interaction between a neuron's weights (W) and the input element. It also has a significant impact on the network's complexity and performance, and it was chosen because it provides ideal decision biases (b). The tansig transfer function can understand the complex non-linear connection between the input and output parameters by producing values ranging from −1 to +1 [78]. Figure 3 depicts the design and architecture of the ANN model development.
neural network with several hundred weights [77]. The model uses a hyperbolic tangent sigmoid (tansig) transfer function as the driving component for the interaction between a neuron's weights (W) and the input element. It also has a significant impact on the network's complexity and performance, and it was chosen because it provides ideal decision biases (b) . The tansig transfer function can understand the complex non-linear connection between the input and output parameters by producing values ranging from −1 to +1 [78]. Figure 3 depicts the design and architecture of the ANN model development.  One of the most frequently used ANN models is the BP-NN, which is a multiple layer feed-forward ANN trained using the error BP method. BP-NN, on the other hand, has weaknesses, which may be fixed by combining it with the Particle Swarm Optimization (PSO) approach, which improves the model's accuracy and efficiency [79]. The PSO approach is utilized in this research in order to optimize the link weights of the ANN.
The PSO is a type of swarm intelligence approach used in evolutionary computing. These techniques were influenced by natural bio-social phenomena such as flocks of birds, schools of fish, and other natural bio-social phenomena. The PSO is particularly applicable for non-linear generalization capabilities with discontinuities because of its quick convergence and robustness. PSO is a promising choice for optimization modelling when compared to other evolutionary algorithms [80,81]. Due to its faster learning speed and lower memory demand, the PSO is preferred over other optimization algorithms including the Genetic Algorithm (GA) and Imperialist Competition Algorithm (ICA) [82].

Hybrid NN-PSO Model
The Neural Network (NN) model's connection weights are optimized via Particle Swarm Optimization (PSO). The PSO is utilized because it can determine the best solution while also reducing the ANN's errors. The PSO calculates the positions of the particles and transmits it to the learning process. The optimal weights and biases for the training method of the ANN were found using the PSO-generated particle population [83]. The framework of the application of the NN-PSO method in spatial interpolation techniques is presented in Figure 4.

Performance Evaluation
The correlation coefficient (R) and mean squared error (MSE) were utilized to assess the performance of the NN-PSO model. A complete positive correlation is implied by a R value of 1. The R value indicates how closely two variables were linked [84]. The R values were observed and utilized as performance indicators throughout the validation and testing phases. The R value was utilized in the validation phase to assess network generalization, terminate the simulation when generalization ceased to improve, and determine the optimum architecture, while the R value in the testing phase serves as an additional independent measure of network performance during and after the simulation [85,86]. In a simulation using the NN-PSO method, the MSE was minimized which includes the overall MSE in the validation and testing phases. The MSE is a helpful tool for assessing model predictions since it reflects the sum of squared bias and variance. For the model, zero is the optimum value for the MSE [87]. The equations for the R and MSE are shown in Equations (1) and (2) where "N" is the number of data sets, y 0 is the predicted value, y m is the measured value, y and y m were the mean values, and e i is the contrast between the measured values and the predicted values [88].
The Neural Network (NN) model's connection weights are optimized via Particle Swarm Optimization (PSO). The PSO is utilized because it can determine the best solution while also reducing the ANN's errors. The PSO calculates the positions of the particles and transmits it to the learning process. The optimal weights and biases for the training method of the ANN were found using the PSO-generated particle population [83]. The framework of the application of the NN-PSO method in spatial interpolation techniques is presented in Figure 4.

Spatial Interpolation Methods for Heavy Metals
The measured physicochemical parameter and detected HM concentrations were mapped by utilizing the ArcGIS platform. The sampling sites' precise locations were recorded using a GARMIN Montana 650 GPS, which was integrated onto the Geographical Information System (GIS) platform. Moreover, data collected from the GW samples were applied to create maps operating the Geostatistical Analyst Tool and Geostatistical Wizard in the ArcGIS software. The different interpolation methods are utilized to apply spatial analysis which are included in the ArcGIS spatial analysis extension tool. Deterministic techniques, geostatistical methods, and interpolation with barriers are among the three (3) interpolation types utilized in the study. The deterministic techniques include Inverse Distance Weighting (IDW), Global Polynomial Interpolation (GPI), Radial Basis Functions (RBF), and Local Polynomial Interpolation (LPI). The IDW is a deterministic spatial interpolation approach that determines the data in an unsampled location by using the data from a distributed collection of sampled locations. The data in an unidentified site are based on the weighted sum of the values of the recognized locations which is based on the distance of the unidentified site to the sampled locations [89]. GPI is a deterministic and approximate trend surface analysis wherein a smooth two-dimensional polynomial function of first, second, or higher degree is used to describe a surface. It computes the target point's value by using all nearby points [90]. RBF achieves its accuracy by using a large number of accurate interpolators that reduce the overall curvature of the surface depending on the space between the interpolated and sampled locations [91]. The LPI method includes fitting the weighted least squares to store the data inside the search ellipse of the grid node wherein the projected value is used to calculate the surface value of the neighboring points that can be used to build surface that account for short-range variation [92].
The geostatistical techniques include Ordinary Kriging (OK), Universal Kriging (UK), and Empirical Bayesian Kriging (EBK). OK is a linear geostatistical process that depends less on stationary mean assumptions by using the search radius. The OK method approximates values in unsampled regions by averaging nearby data and visualizing the correlations between surrounding values as a function of the geographic distance between the sites in the area of study using a weighted average of neighboring data and a variogram [93]. UK uses a trend surface that may include factors that account for variation in the global component, and it more likely to provide residuals that are more closely related to a stationary mean with identical distribution [94]. EBK automates the most time-consuming and challenging stages in creating a realistic kriging model. EBK automatically optimizes by subsetting and simulating several semi-variogram models instead of a single semivariogram. EBK creates a semi-variogram model from existing data and then simulates the new value at each input data point until the final calculation of the new semi-variogram model based on the simulated data [95].
The interpolation with barriers includes Kernel Smoothing (KS) and Diffusion Kernel (DK). The DK employs a complex distance metric specified by the cost surface, which is a widely used raster function that estimates the cost of traveling from one cell of a raster to the next, and then generates forecasts on automatically chosen grids. The KS method is a variation of first order LPI that avoids computing uncertainty by using a technique similar to that used to estimate regression coefficients in ridge regressions. KS utilizes the shortest distance between locations to connect places on opposing sides of an absolute barrier using a sequence of straight lines [96]. These interpolation techniques were implemented to the data arrays in order to distinguish the concentration map that best describes the HM pollution in the province of Marinduque [97].

Cross Validation
A frequently utilized technique for comparing the interpolation methods is cross validation. Due to the small sample size, cross validation was used. A cross validation procedure consists in removing data points one at a time, interpolating a value from the remaining observations, and comparing that value to the real value [98].
The Mean Absolute Error (MAE) and Pearson's Correlation Coefficient (R) were used to determine the predictive accuracy of distinct methods, with the least MAE representing the most exact predictions. Equation (3) shown below is used to calculate MAE: where Z i is the predicted value, Z is the observed value, and n is the number of observations [99].

Results
The major findings of the study are presented in this section. This section contains all the maps created with the integrated NN-PSO and spatial interpolation algorithms. This part also includes the data's descriptive statistics as well as the maps' prediction performance. Table 3 shows descriptive data for GW physicochemical characteristics and HM concentrations. During the dry season, all GW physicochemical parameters and HM concentrations were observed to be highly variable, except for pH and temperature, which were found to be low and moderately variable, respectively. The Kolmogorov-Smirnov test revealed that the dry season physicochemical characteristics and HM concentrations of GW in the research region were not uniformly distributed, with all parameters having p < 0.05.

Heavy Metals in Groundwater
The GW physicochemical characteristics were compared to Philippine National Standards for Drinking Water (PNSDW) 2017 and World Health Organization (WHO) Drinking Water Standards. Table 3 shows that the pH of the GW during the dry season ranged from 6.10 to 7.90, with an average pH of 7.02, which is within the PNSDW 2017 and WHO standards. The EC of the GW during the dry season varies from 80 to 2350 µS/cm with mean EC of 935.17 µS/cm, which is below the maximum value set by the WHO. The asymmetries of the physicochemical characteristics and HM concentration of the GW during the dry season were measured by skewness and kurtosis. With the exception of pH, all physicochemical characteristics and HM concentrations in GW exhibit a positive skewness, meaning that the right side is longer than the left. This indicates that an asymmetry distribution with a positive skewness tends to be less concentrated. TDS, Cr, Fe, Mn, and Zn are the only elements with positive kurtosis values, indicating a steeper distribution than normal.
A comparative assessment of HM concentrations of GW during the dry season was likewise performed relative to the PNSDW 2017 and WHO Standards for Drinking Water. Mean concentrations of Ni ranging from 0.000110 to 0.125310 ppm and Cu ranging from 0.000868 to 0.260497 ppm were below the limit of PNSDW 2017 and WHO standards for drinking water. Zn concentration ranged between 0.000985 and 56.96133 ppm, with mean concentration exceeding the WHO limit. Average concentrations of Cr, Cd, Fe, and Pb exceeded the allowable limits by PNSDW 2017 and WHO Standards. Moreover, some locations had concentrations way above the limit.
During the wet season, all GW physicochemical parameters except for pH and temperature and HM concentrations showed considerable variability. The Kolmogorov-Smirnov test revealed that the wet season physicochemical characteristics and HM concentrations of GW in the research region were not normally distributed, with p < 0.05 for all parameters.
Wet season GW physicochemical characteristics were also compared to PNSDW 2017 and WHO Drinking Water Standards. Table 4 indicates that the pH of the GW throughout the wet season ranged from 6.00 to 9.55, with a mean pH of 7.43, which is within the PNSDW 2017 and WHO standards. The EC of the GW during the wet season varies from 70 to 2640 µS/cm with average 780.61 µS/cm which is inside the acceptable value set by the WHO. The asymmetries of the physicochemical characteristics and HM concentration in GW during wet season were measured by skewness and kurtosis. With the exception of Cr, Cd, Pb, and Cu, all physicochemical characteristics and HM concentrations in GW exhibited positive skewness, meaning the right side is longer than the left side. This indicates that an asymmetry distribution with a positive skewness tends to be less concentrated. All parameters have a negative kurtosis value, which suggests that the distribution of the datasets was flatter than a normal distribution.
The HM concentrations in GW during the wet season were assessed in comparison to the PNSDW 2017 and WHO Standards for Drinking Water. Only the mean Ni content was below the PNSDW 2017 and WHO drinking water requirements. Cr, Cd, Fe, Mn, Pb, Zn, and Cu values were above the permissible limits set by PNSDW 2017 and WHO. Furthermore, some sites have concentrations that were far higher than the permissible level.
A Pearson's Correlation Matrix (PCM) was utilized to establish the level of correlation between GW HMs and physicochemical properties in the island province of Marinduque, with the goal of identifying a potential source of the HMs. Table 5 shows the metals correlation matrix that was generated during the dry season. Cd (r = 0.72), Ni (r = 0.69), and Pb (r = 0.81) were all highly linked with chromium. Cadmium and Ni (r = 0.78), Cd and Pb (r = 0.83), and Ni and Pb (r = 0.77) showed substantial positive associations. Positive significant correlations imply that these metals have a shared origin, were mutually dependent, and behaved similarly throughout transport [100].
The correlation analysis for water quality parameters during the wet season is shown in Table 6 and illustrated in Figure 5. For the physicochemical characteristics, the EC was observed to have a positive correlation to TDS which agreed to the findings of Manikandan et al. in 2020 [101]. A significant connection was found between Cr and Cd, as well as Cu and Zn, for the HMs in GW during the wet season. This indicated a possible shared source for these HMs. Moreover, these correlations were in agreement with the findings of Kumar    1 . 0 0 ** Correlation is significant at the 0.01 level (two-tailed); * Correlation is significant at the 0.05 level (two-tailed).

NN-PSO Simulation Results
The hybrid NN-PSO was utilized to improve the efficiency and robustness of spatial interpolation mapping of GW quality. The proposed method was evaluated by considering different internal characteristics of the network. The Levenberg-Marquardt method was used as the training algorithm, and the hyperbolic tangent sigmoid was used as the transfer function for the input layer (IL) to hidden layer (HL) and HL to output layer (OL) in the hybrid NN-PSO informed spatial interpolation approaches. The results of the NN-PSO simulation for the dry and wet season GW quality parameters are presented in Tables 7 and 8.

NN-PSO Simulation Results
The hybrid NN-PSO was utilized to improve the efficiency and robustness of spatial interpolation mapping of GW quality. The proposed method was evaluated by considering different internal characteristics of the network. The Levenberg-Marquardt method was used as the training algorithm, and the hyperbolic tangent sigmoid was used as the transfer function for the input layer (IL) to hidden layer (HL) and HL to output layer (OL) in the hybrid NN-PSO informed spatial interpolation approaches. The results of the NN-PSO simulation for the dry and wet season GW quality parameters are presented in Tables 7 and 8. The results and efficiency of the NN-PSO simulation for dry and wet season GW quality parameters showed excellent results, as evidenced by the extremely low MSE (ideal value is zero) and extremely high R values for internal validation and testing of the NN-PSO models (ideal value is one). These NN-PSO models informed the spatial interpolation techniques which improved the accuracy and performance of the GW quality maps. The R plots of the NN-PSO simulation for the physicochemical parameters and HM concentrations during the dry and wet season are presented in Figures A1-A6 in Appendix A.

NN-PSO Informed Spatial Interpolation Techniques for GW Quality Mapping
The performance of the NN-PSO informed spatial interpolation approaches for GW quality mapping was evaluated via cross validation, with MAE and R values utilized in order to assess the prediction capability of the various interpolation techniques. The approach that produced the lowest MAE value and the greatest R value was regarded the best. The accuracy and effectiveness of the techniques employed to interpolate GW quality during the dry season are displayed in Figures 6 and 7. The complete values for the performance of the different interpolation techniques during the dry season are shown in Table A1 of Appendix B. Various interpolation techniques were applied in order to evaluate the spatial variability of GW quality in the province of Marinduque during the dry season, including IDW, GPI, RBF, LPI, OK, UK, EBK, DK, and KS. Moreover, these interpolation techniques were likewise assessed after being informed and integrated with NN-PSO.
Physicochemical properties and HM concentrations of GW were mapped and analyzed using various interpolation approaches throughout the dry season. Except for Nickel (Ni), which had the best mapping prediction performance using NN-PSO informed radial basis functions, NN-PSO informed geostatistical techniques and OK and EBK methods, in particular, were the best approaches for mapping the GW quality during the dry season. The performance of the best interpolation technique was manifested through its lowest MAE and highest R values. Figure 6 illustrates that the NN-PSO informed OK model had the best mapping prediction for GW temperature and EC during dry season. The OK-NN-PSO method had the least MAE and highest R value which is significantly higher than the best method observed without NN-PSO. The EBK-NN-PSO method was observed to have the best prediction performance for groundwater pH and TDS as evident to its lowest MAE and highest R value. Similarly, these validation criteria have significant improvement than compared to the interpolation methods that were not NN-PSO informed.
During the dry season, HM concentrations in GW were also mapped and evaluated by using several spatial interpolations approaches as shown in Figure 7. Using the NN-PSO informed Empirical Bayesian Kriging technique (EBK-NN-PSO), the best mapping prediction performance was found for Cr, Cd, Fe, Pb, Zn, and Cu. Various interpolation techniques were applied in order to evaluate the spatial variability of GW quality in the province of Marinduque during the dry season, including IDW, GPI, RBF, LPI, OK, UK, EBK, DK, and KS. Moreover, these interpolation techniques were likewise assessed after being informed and integrated with NN-PSO.
Physicochemical properties and HM concentrations of GW were mapped and analyzed using various interpolation approaches throughout the dry season. Except for Nickel (Ni), which had the best mapping prediction performance using NN-PSO informed radial basis functions, NN-PSO informed geostatistical techniques and OK and EBK methods, in particular, were the best approaches for mapping the GW quality during the dry season. The performance of the best interpolation technique was manifested through its lowest MAE and highest R values. Figure 6 illustrates that the NN-PSO informed OK model had the best mapping prediction for GW temperature and EC during dry season. The OK-NN-PSO method had the least MAE and highest R value which is significantly higher than the best method observed without NN-PSO. The EBK-NN-PSO method was observed to have the best prediction performance for groundwater pH and TDS as evident to its lowest MAE and highest R value. Similarly, these validation criteria have significant improvement than compared to the interpolation methods that were not NN-PSO informed.
During the dry season, HM concentrations in GW were also mapped and evaluated by using several spatial interpolations approaches as shown in Figure 7. Using the NN-PSO informed Empirical Bayesian Kriging technique (EBK-NN-PSO), the best mapping prediction performance was found for Cr, Cd, Fe, Pb, Zn, and Cu. Toxics 2021, 9,   Among the interpolations employed, the EBK+NN-PSO technique for mapping Cr, Cd, Fe, Pb, Zn, and Cu has the lowest MAE. The R value of the NN-PSO informed EBK method used for Cr increased by 35.61%, Cd = 17.91%, and Pb = 21.65%. Furthermore, R values of Fe, Zn, and Cu improved by 5.6, 2.4, and 2.9 times, respectively. Manganese (Mn) and Nickel (Ni) were observed to have the best prediction performance using NN-PSO informed OK and RBF method, respectively. Correspondingly, the performance of these models for Mn and Ni was observed to have the lowest MAE among the interpolation techniques utilized and 3.1 and 1.2 times higher than compared to interpolation techniques without NN-PSO. The NN-PSO informed dry season GW physicochemical characteristics and HM concentration map of Marinduque are presented in Figures A7 and A8 of Appendix C.
The NN-PSO informed geostatistical approaches, comprising OK and EBK, were shown to be the optimum method for the wet season GW physicochemical properties, with the lowest MAE and highest R. For groundwater temperature and pH, OK was the best technique, whereas EBK was the best method for EC and TDS. The NN-PSO informed OK was the best among the interpolation techniques integrated with NN-PSO and had an R observed to be 3.9 times greater than compared to the best interpolation method without NN-PSO. Additionally, the other physicochemical parameters including GW pH, EC, and TDS were observed to have an R value of 2.7, 4.3, and 3.4 times higher than the spatial interpolation methods with NN-PSO integration. The performance of the different spatial interpolation techniques for the parameters observed during the wet season is presented in Figures 8 and 9. As illustrated in Figure 8, the OK+NN-PSO model provided the most accurate mapping forecast for GW temperature and pH during wet season. Among the observed spatial interpolation techniques for temperature and pH, the OK+NN-PSO approach achieved the highest R value and the lowest MAE. The EBK+NN-PSO approach had shown the greatest prediction performance for GW EC and TDS during wet season.  Figures A7 and A8 of Appendix C. The NN-PSO informed geostatistical approaches, comprising OK and EBK, were shown to be the optimum method for the wet season GW physicochemical properties, with the lowest MAE and highest R. For groundwater temperature and pH, OK was the best technique, whereas EBK was the best method for EC and TDS. The NN-PSO informed OK was the best among the interpolation techniques integrated with NN-PSO and had an R observed to be 3.9 times greater than compared to the best interpolation method without NN-PSO. Additionally, the other physicochemical parameters including GW pH, EC, and TDS were observed to have an R value of 2.7, 4.3, and 3.4 times higher than the spatial interpolation methods with NN-PSO integration. The performance of the different spatial interpolation techniques for the parameters observed during the wet season is presented in Figures 8 and 9. As illustrated in Figure 8, the OK+NN-PSO model provided the most accurate mapping forecast for GW temperature and pH during wet season. Among the observed spatial interpolation techniques for temperature and pH, the OK+NN-PSO approach achieved the highest R value and the lowest MAE. The EBK+NN-PSO approach had shown the greatest prediction performance for GW EC and TDS during wet season.  As illustrated by Figure 9, EBK+NN-PSO provided the best mapping prediction performance for Fe, Mn, Ni, Pb, and Zn. OK+NN-PSO provided the best mapping prediction performance for Cr, Cd, and Cu. RBF+NN-PSO provided the best mapping prediction performance for Zn.  As illustrated by Figure 9, EBK+NN-PSO provided the best mapping prediction performance for Fe, Mn, Ni, Pb, and Zn. OK+NN-PSO provided the best mapping prediction performance for Cr, Cd, and Cu. RBF+NN-PSO provided the best mapping prediction performance for Zn. Table A2 in Appendix B elaborates the complete values for the performance of the various interpolation methods for the wet season.
The HM concentration during the wet season was mapped using the different interpolation techniques with and without the integration of NN-PSO. Chromium (Cr), Cadmium (Cr), and Copper (Cu) concentrations were mapped using the best method observed which is the integration of Ordinary Kriging (OK) and NN-PSO. The OK+NN-PSO technique had the lowest MAE, and its R value was 4.7, 1.2, and 1.3 times greater than the maximum R value found for the interpolation method without NN-PSO integration. The metals Fe, Mn, Ni, and Pb have the best prediction performance based on MAE and R utilizing Empirical Bayesian kriging (EBK) combined with NN-PSO. When compared to interpolation approaches without NN-PSO integration, the EBK-NN-PSO method provided the lowest MAE and R values, which were 2.5, 1.5, 1.4, and 2.4 times higher. The spatial interpolation maps with the best prediction performance for the physicochemical characteristics and HM concentrations are presented in Figures A9-A11 of Appendix C.
The summary of the cross-validation performance of the models for the physicochemical parameters and the HM concentrations both for dry and wet season is exhibited in Table 9.

Discussion
Based on the results, the use of the hybrid neural network-particle swarm optimization method in spatial interpolation was able to address the sampling density issue experienced in water quality monitoring. It provides solutions on data gaps when the spatial distribution map of GW quality becomes available by having simple water parameter such as pH [104]. The impact of this research is emphasizes on prediction improvement for mapping specific features and water quality parameters and ease in GW quality monitoring. Furthermore, it contributes to solutions in data gaps when processed data are necessary.
The application of a range of machine learning techniques, as well as the combination and comparison of these approaches, resulted in a larger pool of potential environmental monitoring systems. The integration of machine learning techniques to spatial interpolation techniques have been implemented in several studies including Least Squares Support Vector Machine (LSSVM) and Population-based Incremental Learning to Ordinary Kriging (OK) [105], Random Forest (RF) to IDW and OK [106], Deep Neural Network (DNN) to Ordinary Kriging (OK) [107], Decision Tree to Kriging and Inverse Distance Weighting [108], and Non-linear AutoRegressive eXogenous (NARX) model to Geographic Information System [109].
A different method was integrated relative to spatial interpolation techniques in this research study. The use of the NN-PSO methodology with the capability of both prediction and optimization was utilized to enhance the prediction capacity of spatial interpolation methods. The findings of the cross validation of the HM concentration and physicochemical parameters showed that the spatial interpolation methods integrated with NN-PSO were the best method as manifested to its lowest MAE and highest R value. The governing interpolation methods were mostly under the geostatistical methods which integrated with NN-PSO including OK (OK+NN-PSO) and EBK (EBK+NN-PSO) method. Moreover, some HMs such as Ni (dry season) and Zn (wet season) have best performance using NN-PSO informed radial basis functions (RBF-NN-PSO). The findings of this study in reference to the performance of the NN-PSO informed spatial interpolation techniques agreed with the study of Li et al. [110], wherein they confirm that the integration of machine learning techniques produces more superior performance in the spatial interpolation method than compared to spatial interpolation methods without machine learning integration. The NN-PSO integration relative to spatial integration techniques addressed the issue in sampling density and was able to improve the performance of the spatial interpolation methods.

Conclusions
The objective of this research is to improve the mapping prediction capability of spatial interpolation algorithms by using an NN-PSO technique. It analyzed interpolation methods and mapping for physicochemical parameters such as temperature, pH, TDS, and EC, as well as mapping of HM concentrations in GW. This technique comprises three spatial interpolation methods such as deterministic, geostatistical, and interpolation with barriers with neuro-particle swarm optimization guided interpolation approaches. The measurement criteria for the best method were the least MAE and the highest R value.
The results recorded the governing interpolation techniques during the dry and wet season. The OK+NN-PSO method was recorded as best performing for temperature and EC during dry season, while EBK+NN-PSO for pH and TDS. The best methods during wet season were OK+NN-PSO for temperature and pH and EBK+NN-PSO for EC and TDS. These methods have the highest R and lowest MAE among the spatial interpolation techniques observed. The best methods for mapping physico-chemical characteristics were found to have the least MAE and R values ranging from 1.7 to 5.6 times higher than the interpolation techniques without NN-PSO integration.
The HM concentration maps during the dry season were observed to have the best performance using EBK+NN-PSO for Cr, Cd, Fe, Pb, Zn, and Cu while RBF+NN-PSO was the best method for Ni mapping. The best method during the wet season was found to be OK+NN-PSO for Cr, Cd, and Cu; EBK+NN-PSO for Fe, Mn, Ni, and Pb; and RBF+NN-PSO for Zn. The best method for mapping the GW HM concentrations during the wet season was observed to have the least MAE and R values ranging from 1.2 times to 5.6 times greater than the best interpolation method without NN-PSO integration.
Hybrid methods in general showed better performance than compared to the nonhybrid methods. The development of these hybrid methods using NN-PSO and geostatistics provides a promising innovative approach for environmental quality monitoring as it improves the accuracy of predictive mapping and modelling of GW quality in an area. The integration of NN-PSO into spatial interpolation methods addresses the challenge of sample density and its effect on the spatial interpolation method's performance. It opens a new avenue for enhancing the predictive capability of spatial interpolation algorithms. On the basis of the results of this research, it can be stated that the employment of models such as NN-PSO is suitable for overcoming the challenges in spatial data processing and mapping, as well as for improving the model's predictive capabilities. The findings of the study suggest that the integration of ML techniques such as NN can be utilized in mapping GW quality as well as its application in spatio-temporal maps.  Acknowledgments: This is to acknowledge the support-in-kind of Mapua University, Marinduque Local Government Units, and Marinduque State College.

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

Appendix A
OK+NN-PSO for Cr, Cd, and Cu; EBK+NN-PSO for Fe, Mn, Ni, and Pb; and RBF+NN-PSO for Zn. The best method for mapping the GW HM concentrations during the wet season was observed to have the least MAE and R values ranging from 1.2 times to 5.6 times greater than the best interpolation method without NN-PSO integration.
Hybrid methods in general showed better performance than compared to the nonhybrid methods. The development of these hybrid methods using NN-PSO and geo-statistics provides a promising innovative approach for environmental quality monitoring as it improves the accuracy of predictive mapping and modelling of GW quality in an area. The integration of NN-PSO into spatial interpolation methods addresses the challenge of sample density and its effect on the spatial interpolation method's performance. It opens a new avenue for enhancing the predictive capability of spatial interpolation algorithms. On the basis of the results of this research, it can be stated that the employment of models such as NN-PSO is suitable for overcoming the challenges in spatial data processing and mapping, as well as for improving the model's predictive capabilities. The findings of the study suggest that the integration of ML techniques such as NN can be utilized in mapping GW quality as well as its application in spatio-temporal maps.     Appendix B