In Situ Measurements of Domestic Water Quality and Health Risks by Elevated Concentration of Heavy Metals and Metalloids Using Monte Carlo and MLGI Methods

The domestic water (DW) quality of an island province in the Philippines that experienced two major mining disasters in the 1990s was assessed and evaluated in 2021 utilizing the heavy metals pollution index (MPI), Nemerow’s pollution index (NPI), and the total carcinogenic risk (TCR) index. The island province sources its DW supply from groundwater (GW), surface water (SW), tap water (TP), and water refilling stations (WRS). This DW supply is used for drinking and cooking by the population. In situ analyses were carried out using an Olympus Vanta X-ray fluorescence spectrometer (XRF) and Accusensing Metals Analysis System (MAS) G1 and the target heavy metals and metalloids (HMM) were arsenic (As), barium (Ba), copper (Cu), iron (Fe), lead (Pb), manganese (Mn), nickel (Ni), and zinc (Zn). The carcinogenic risk was evaluated using the Monte Carlo (MC) method while a machine learning geostatistical interpolation (MLGI) technique was employed to create spatial maps of the metal concentrations and health risk indices. The MPI values calculated at all sampling locations for all water samples indicated a high pollution. Additionally, the NPI values computed at all sampling locations for all DW samples were categorized as “highly polluted”. The results showed that the health quotient indices (HQI) for As and Pb were significantly greater than 1 in all water sources, indicating a probable significant health risk (HR) to the population of the island province. Additionally, As exhibited the highest carcinogenic risk (CR), which was observed in TW samples. This accounted for 89.7% of the total CR observed in TW. Furthermore, all sampling locations exceeded the recommended maximum threshold level of 1.0 × 10−4 by the USEPA. Spatial distribution maps of the contaminant concentrations and health risks provide valuable information to households and guide local government units as well as regional and national agencies in developing strategic interventions to improve DW quality in the island province.


Introduction
The fundamental requirement for human growth and development is water. Around 71% of the world's population relies on a clean drinking water supply that is readily The Monte Carlo simulation (MCS) has been recognized by the United States Environmental Protection Agency (USEPA) as a practical approach for uncertainty analysis of an RA, and it is frequently used for HMM RAs. RAs for such environmental compartments have become more significant as the levels of HMM in water supplies have increased [19,20]. Marinduque, which hosted catastrophic mining disaster events, has a limited study on human HR analysis to exposure to HMM in the DW. Most of the population relies on SW and GW for DW supply. Moreover, some water refilling stations (WRS) used SW and GW as the source of water for their bottled and refilling business. Hence, an RA considering uncertainty analysis and probabilistic HR evaluation methodologies such as MCS is essential.
Regarding the stated difficulties of HMM contamination of water resources, in situ, precise, and real-time detection of HMMs in water are needed. Laboratory-based methods such as inductively coupled plasma optical emission spectroscopy (ICP-OES) and atomic absorption spectroscopy (AAS) have drawbacks, as these methods need several days to complete the analyses. These laboratory analytical techniques are also costly, involve complicated sample preparation, and are inapplicable in field environments [21]. The utilization of portable X-ray fluorescence (pXRF) and metals analysis systems (MASs) for in situ metals detection addresses these limitations by providing a speedy and reliable analysis. This in situ HMM detection and analysis, coupled with the MCS and mapping, produces a rapid and reliable RA. Hence, the primary purpose of this research was to provide new knowledge to the population of the quality of their DW supply, the available onsite detection, and analysis technology. The specific objectives of this study were to: (1) assess HMM pollution at specific sites using various indices such as MPI and NPI; (2) evaluate the carcinogenic and noncarcinogenic HRs with the application of MCS to analyze the uncertainty in the results; and (3) create a machine learning-geostatistical interpolation (MLGI)-based spatial interpolation maps of pollution and health risk index of DW in Marinduque. The significance of this study is the provision of benchmark data for local DW quality monitoring specifically the HMM concentrations and the health risks to go along with the DW supply quality. In addition, this information provides the groundwork for the development, utilization, and protection of water resources.

Study Area
The DW samples were gathered across the island province of Marinduque, Philippines, as shown in Figure 1. Marinduque is located about 200 km south of Manila, the capital city of the Philippines, covering about 96,000 hectares of land [22]. It comprises six municipalities: Boac, Buenavista, Gasan, Mogpog, Sta. Cruz, and Torrijos. It is a tropical and warm island with an annual mean temperature of 23 to 28 • C [23]. The geology of the island is characterized by permeable volcanic and sedimentary rocks, which allow groundwater to flow through across the island [14,24].

Water Samples Collection, Preparation, and Analysis
One hundred grab samples in a 500 mL volume were collected using the suggested standard procedure [25]. Twenty-five water samples from WRS, twenty-six groundwater (GW) samples, and twenty-one tap water (TW) samples were collected and placed in polyethylene (PE) bottles. All water samples were either used for drinking and/or used for cooking by households. The PE bottles were washed with distilled water or Type 1 water [26] to remove contaminants. Physicochemical parameters such as temperature, pH, total dissolved solids (TDS), and electrical conductivity (EC) were evaluated and recorded in situ using a Hannah Multiparameter HI 9811-5 portable meter [27]. The portable handheld Olympus Vanta X-ray fluorescence spectrometer and Accusensing Metals Analysis System (MAS) G1 were utilized for the HMM concentration analysis. Both analyzers are highperformance and on-site elemental analyzers that may be used with various environmental media, including water. The pXRF spectrometer was calibrated using the Olympus Vanta blank and set to Geochem mode before analysis [21]. The Accusensing MAS G1 was used for analyzing HMMs which registered "limits of detection" (LOD) in the pXRF such as As, Cu, Ni, and Pb. Figure 2 shows the metals analyzed using the two analyzers.

HM Pollution and HR Assessment
The detected elements in the water samples were assessed and compared to the Philippine National Standards for Drinking Water (PNSDW) 2017 [28]. A HMM pollution assessment utilizing several indices such as the MPI, NPI, non-CR Index, and CR Index was implemented in this study to comprehensively assess the degree of pollution in the DW considering HMMs.

HMM Pollution Index (MPI)
The MPI describes the WQ state and determines whether it is suitable for drinking concerning HMMs. The MPI is established on the weighted arithmetic mean quality procedure as presented in Equation (1).
where W i is denoted as the weight unit, which can be computed as 1/S i where S i is the suggested level of the pertinent HMM, "n" is the value of the evaluated HMMs, and Q i is the distinctive quality rating of the "ith" metal and specified by Equation (2).
C i is the detected value of the ith metals in micrograms/liter. The standard allowable value (S i ) for each constraint was acquired from the Philippine water quality standards [29]. The categorization of the WQ using MPI is demonstrated in Table 1. The Single Factor Pollution Index (SFPI) was employed to calculate and describe the effect of an individual HMM as it contaminates the DW at a specific sampling location. The SFPI is computed using Equation (3) presented below: where C i denotes the observed intensity of the pollutant "i" in the groundwater (in mg/L) and S i is the evaluation standard of the contaminant "i" in the surface and groundwater (in mg/L). An SFPI value greater than 1 signifies that the heavy-metal pollutant exceeds the standard [31].

Nemerow's Pollution Index (NPI)
The NPI was utilized to calculate and explain the influence of several HMMs as they pollute the water resource at a particular sampling location. The NPI is the most usual method and thorough pollution assessment approach that reflects the single factor index P i , highlights the impact of elevated levels of HMMs on the environmental quality, and Toxics 2022, 10, 342 6 of 32 removes the insufficiency of the mean value on assessment. The NPI is calculated using Equation (4), shown below: where SFPI max indicates the highest SFPI of a pollutant and SFPI ave denotes the mean SFPI of the pollutants considered [32,33]. The classification of NPI values is shown in Table 2. HMMs are introduced into the body by various pathways, with ingestion-the most frequent route of water consumption-as a common pathway resulting from oral human exposure. The CDI quantifies the extent of pollution absorbed by the human body and specifies the pollutant dosage per kilogram of body weight per day as received by direct eating, dermal absorption, or inhalation, as suggested by the USEPA. The CDI of water ingestion can be calculated using Equation (5) [35]: where CDI in is the exposure doses from ingestion of water in gram/kilogram-day and C i is the mean concentration of the ith HMM in water (micrograms/liter) [36]. Additional amounts and units of other parameters in the computation of CDI are presented in Table 3. The relevant RfD was compared to the exposure or mean intake of hazardous elements to determine the probability on noncarcinogenic substances. The noncancer risk was quantified using the hazard quotient (HQ) for a single chemical or the hazard index (HI) for multiple substances and exposure routes. Concerns about possible noncarcinogenic consequences may arise if exposure to a chemical exceeds the corresponding RfD, i.e., if HQ exceeds 1 [43] as shown as Equation (6). Table 4 presents the toxicity reactions of HMMs for RfD values for both oral and dermal exposure route as well as the oral slope factor (SF). The total potential for non-CRs influenced by calculating chemicals can be evaluated by the HI, which is the sum of each computed HQ. Equation (7) displays the formula for calculating the hazard index.
It is recommended to have even greater probabilities of harmful health effects when the HI is greater than 1. At the same time, no chronic risks were expected to transpire at the site when HI was less than 1 [48].
The USEPA defined CRs as the cumulative risk of a person developing cancer because of exposure to a probable cancer-causing agent. The cancer risk was calculated by Equation (8): The total cancer risk (TCR) is the sum of the cancer risks due to the ingestion exposures to multiple HMMs of concern. The TCR can be evaluated through Equation (9), and the risk value levels are shown in Table 5.

MCS and SA
The sensitivity analysis (SA) approach was used to establish the impact of changing the values of the input variables on the CR estimate under a set of assumptions [50]. The SA was implemented using Oracle Crystal Ball ® version 11.1.

Statistical Analysis
Using Excel software, descriptive statistics linked to the physicochemical parameters and HMM concentrations in the DW were evaluated to calculate the mean, standard deviation (SD), and coefficient of variance (CV). The CV was utilized to evaluate the dataset's variability as follows: CV ≤ 15%, low; 15% < CV ≤ 35%, intermediate; and CV ≥ 35%, high [51][52][53]. A Pearson rank-order correlation analysis coupled with correlation matrix plots was also implemented utilizing IBM SPSS Statistics Version 26.0 and R Studio. Data on the metals concentrations in the domestic water were subjected to a hierarchical cluster analysis (CA) to find homogeneous groups. A dendrogram was likewise produced to analyze the cohesion of the clusters produced, in which relationships among components may easily be identified [54,55].

Spatial Concentration Mapping of Risk Indices Using Machine Learning Geostatistical Interpolation (MLGI) Method
A hybrid neuro-particle swarm optimization (NN-PSO) technique coupled with the empirical Bayesian kriging (EBK) method was utilized to generate spatial maps of the risk indices of domestic water samples from the risk indices calculated at each sampling location [21]. Table 6 shows the average mean, SD, and CV of the DW's physicochemical parameters and HMM concentrations. These were assessed with WQ standards for drinking water regulated by WHO [1] and PNSDW [56]. The mean pH and TDS of all water samples, except for the mean TDS of TW, were below or within the standards. High TDS may affect the taste and palatability of drinking water [57]. The CV for the TDS and EC of all water samples were higher than 35%, indicating a greater dispersion around the mean and a relatively high variability in the data sets [51]. The measured HMM concentrations in the DW shown in Figure 3 varied between 0 and 19.0 mg/L for WRS, 0 and 13.4 mg/L for GW, and 0 and 9.31 mg/L for TW. The mean concentrations of As, Pb, and Ni in WRS, GW, and TW were higher than the acceptable limits of the WHO and PNSDW standards. These elevated concentrations of HMM were attributed to the existing two abandoned open mine pits that are located at a higher elevation and sit on permeable volcanic and sedimentary rocks [14,24] that allow groundwater to pass through. Toxics 2022, 10, x FOR PEER REVIEW 2 of 16   The observed HMM concentrations trend, shown in Table 7, for WRS, GW and TW were As > Pb > Fe > Ni > Cu > Zn > Ba > Mn, Pb > Fe > As > Ni > Zn > Ba > Cu > Mn, and As > Pb > Ni >Fe > Zn > Cu > Ba > Mn, respectively. It was recorded that As had the highest concentration detected in both the WRS and TW, while Mn had the lowest concentration observed in all DW sources. Continuous subsurface flow of HMM into inland waters contributes to the increased metal concentration in water resources. In addition, the weathering of rocks that leached HMMs may contribute to the concentration of these HMM.

MPI and NPI Results
MPI and NPI were broadly utilized to evaluate the total HMM contamination in water resources [4]. The MPI calculated in all DW sampling locations was observed to have a high degree of pollution. The average MPI value for the TW samples was 37.5 times more than the minimum MPI value. This is classified as having a high degree of pollution, while the average MPI values for the GW and WRS samples were 33.6 times and 22.3 times greater, respectively. Moreover, the NPI values observed in the DW samples were 8.4 times to 13.6 times higher than the minimum NPI value. These NPI values are categorized as heavily polluted. Table 8 presents the mean CDI of metals through the oral route from DW. The CDI values calculated for adults ranged from 0.0003 to 0.0386 for all DW sources. At the same time, the CDI values for children were observed to range from 0.0004 to 0.0491. The highest mean CDI, illustrated in Figure 4, for adults and children was observed for As for TW and WRS and Pb for GW. The smallest mean CDI for adults was observed for Mn for all DW samples. The study of de Jesus et al. (2021) in Marinduque also revealed high concentrations of Pb in the GW samples [58].    The HQ indices of metals through water intake in the study area are elaborated in Table 9. The average HQ indices of As for adults and children in TW and WRS were highest, while the mean HQ index value for Pb was the highest for GW samples both for adults and children. The HQ index trend of HMM for adult and children in TW and WRS were observed to be As > Pb > Ni > Cu > Fe > Ba > Zn > Mn, while the HQ index in GW for adult is Pb > As > Ni > Fe > Cu > Ba > Mn > Zn. The HQ index trend of GW for children was Pb > As > Ni > Fe > Cu > Ba > Zn > Mn. Furthermore, the HQ indices of As and Pb for all water sources were greater than 1, indicating potential health risks to the human population. It must be emphasized that exposure to HMM may also come from various sources and through other pathways of exposure [59]. The TCR to residents through water intake from TW, GW, and WRS is summarized in Table 10. Among the studied target contaminants, only As, Pb, and Ni are categorized as  [60,61]. The adult carcinogenic risk for As ranged from 1.41 × 10 −3 to 4.96 × 10 −2 ; Pb from 8.92 × 10 −5 to 4.16 × 10 −4 ; and Ni from 2.07 × 10 −3 to 5.50 × 10 −3 . The child carcinogenic risk for As ranged from 1.80 × 10 −3 to 1.50 × 10 −1 ; Pb from 1.14 × 10 −4 to 5.29 × 10 −4 ; and Ni from 2.35 × 10 −3 to 2.76 × 10 −3 . The mean TCR was 2.51 × 10 −2 and 5.95 × 10 −2 for adults and children, respectively. All these risks were greater than the threshold value proposed by USEPA, which is 1 × 10 −4 [62,63]. Having recorded these values, certain interventions and control measures are required to reduce the level of concentrations of HMM in the province and limit the population's exposure. Adequate remediation and prompt onsite treatment to safeguard people's health are necessary [64]. The TCR of the carcinogenic metals in all water sources were seen in the order As > Ni > Pb, and the TCR was in the order TW > WRS > GW both for adult and children, as shown in Figure 5.

Monte Carlo Simulation and Sensitivity Analysis
The carcinogenic risk (adults and children) related to As, Pb, and Ni in all residential water sources in Marinduque, Philippines was evaluated using the Monte Carlo method. The likelihood of lifetime cancer risk (adults) for As in all water samples is shown in Figure 6. The mean TCR for As was 3.18 × 10 −2 , while the risks of 5% and 95% were as high as 2.38 × 10 −2 and 4.10 × 10 −2 , respectively. This risk is very high compared to the indicated maximum acceptable risk of 1.00 × 10 −4 for As.   The factors included in the lifetime carcinogenic risk (LTCR) estimate was then determined using a sensitivity analysis (SA). The As in water demonstrated that the two components involving HMM content and body weight (BW) had the most significant effect on the LTCR values. Compared to other factors, an As concentration of 36.4 percent and a BW concentration of 13.0 percent had the highest positive and negative impacts on the LTCR calculation ( Figure A1).
The average likelihood of the LTCR for Ni was 6.86 × 10 −3 , while the risks of 5% and 95% were equal to 5.18 × 10 −3 and 8.82 × 10 −3 , respectively (Figure 7). The SA for the LTCR estimate regarding Ni indicated that the two parameters including a Ni concentration of 36.8% and an AT of −12.6%, respectively, had the highest positive and negative influences on the carcinogenic hazard value, as shown in Figure A2.   Moreover, the average likelihood of lifetime CR for Pb was 1.67 × 10 −4 , while the risks of 5% and 95% were equal to 1.25 × 10 −4 and 2.15 × 10 −4 , respectively, as shown in Figure 8. The SA for the LTCR computation concerning Pb demonstrated that the two parameters consisting of a Pb concentration of 36.0% and a BW of −13.5% had the highest positive and negative impacts on the carcinogenic hazard value as shown in Figure A3. For children, the average likelihood of LTCR of As was observed to be 4.06 × 10 −2 , while the risks of 5% and 95% were equivalent to 3.06 × 10 −2 and 5.24 × 10 −2 , respectively, as shown in Figure 9. The SA for the LTCR (children) computation on As demonstrated that the two parameters comprising an As concentration of 36.8% and a BW of −13.6% had the highest positive and negative impacts on carcinogenic hazard computation as shown in Figure A4. For children, the average likelihood of LTCR of As was observed to be 4.06 × 10 −2 , while the risks of 5% and 95% were equivalent to 3.06 × 10 −2 and 5.24 × 10 −2 , respectively, as shown in Figure 9. The SA for the LTCR (children) computation on As demonstrated that the two parameters comprising an As concentration of 36.8% and a BW of −13.6% had the highest positive and negative impacts on carcinogenic hazard computation as shown in Figure A4.
For children, the average likelihood of LTCR of As was observed to be 4.06 × 10 −2 , while the risks of 5% and 95% were equivalent to 3.06 × 10 −2 and 5.24 × 10 −2 , respectively, as shown in Figure 9. The SA for the LTCR (children) computation on As demonstrated that the two parameters comprising an As concentration of 36.8% and a BW of −13.6% had the highest positive and negative impacts on carcinogenic hazard computation as shown in Figure A4. The average likelihood of LTCR of Ni for children was observed to be 8.72 × 10 −3 , while the risks of 5% and 95% were equivalent to 6.56 × 10 −3 and 1.13 × 10 −2 , respectively as shown in Figure 10. The SA for the LTCR (children) computation involving Ni explained that the two parameters containing a Ni concentration of 36.0% and a BW of −13.3% had the highest positive and negative impacts on the carcinogenic hazard calculation, as shown in Figure A5. The average likelihood of LTCR of Ni for children was observed to be 8.72 × 10 −3 , while the risks of 5% and 95% were equivalent to 6.56 × 10 −3 and 1.13 × 10 −2 , respectively as shown in Figure 10. The SA for the LTCR (children) computation involving Ni explained that the two parameters containing a Ni concentration of 36.0% and a BW of −13.3% had the highest positive and negative impacts on the carcinogenic hazard calculation, as shown in Figure A5. Furthermore, the LTCR of Pb for children was observed to be 2.10 × 10 −4 , while the risks of 5% and 95% were equivalent to 1.60 × 10 −4 and 2.70 × 10 −4 , respectively as shown in Figure 11. The SA for the LTCR (children) computation involving Pb revealed that the two parameters consisting of a Pb concentration of 37.5% and an AT of −12.9% had the highest positive and negative impacts on the carcinogenic hazard calculation, as shown in Figure A6. Furthermore, the LTCR of Pb for children was observed to be 2.10 × 10 −4 , while the risks of 5% and 95% were equivalent to 1.60 × 10 −4 and 2.70 × 10 −4 , respectively as shown in Figure 11. The SA for the LTCR (children) computation involving Pb revealed that the two parameters consisting of a Pb concentration of 37.5% and an AT of −12.9% had the highest positive and negative impacts on the carcinogenic hazard calculation, as shown in Figure A6.
Furthermore, the LTCR of Pb for children was observed to be 2.10 × 10 −4 , while the risks of 5% and 95% were equivalent to 1.60 × 10 −4 and 2.70 × 10 −4 , respectively as shown in Figure 11. The SA for the LTCR (children) computation involving Pb revealed that the two parameters consisting of a Pb concentration of 37.5% and an AT of −12.9% had the highest positive and negative impacts on the carcinogenic hazard calculation, as shown in Figure A6.

Relationship of WQ Parameters in the Domestic Water Samples
The relationships between the selected WQ parameters, especially metals, in the water provided interesting information on their potential sources and pathways. Figure 12a presents the correlation plot for the WQ parameters obtained from the TW samples. It can be seen that there are significant positive correlations between EC and TDS as well as Pb and Cu. For the GW samples, the correlation matrix plot of each WQ parameter is presented in Figure 12b. It is observed that significant positive correlations exist between EC and TDS and Fe and Zn. Moreover, the correlation matrix plot of WQ parameters obtained for WRS is exhibited in Figure 12c. It is observed that significant positive correlations exist between EC and TDS, Fe and Cu, Pb and Cu, and Pb and Fe.

Relationship of WQ Parameters in the Domestic Water Samples
The relationships between the selected WQ parameters, especially metals, in the water provided interesting information on their potential sources and pathways. Figure 12a presents the correlation plot for the WQ parameters obtained from the TW samples. It can be seen that there are significant positive correlations between EC and TDS as well as Pb and Cu. For the GW samples, the correlation matrix plot of each WQ parameter is presented in Figure 12b. It is observed that significant positive correlations exist between EC and TDS and Fe and Zn. Moreover, the correlation matrix plot of WQ parameters obtained for WRS is exhibited in Figure 12c. It is observed that significant positive correlations exist between EC and TDS, Fe and Cu, Pb and Cu, and Pb and Fe.  The observed relationships of metals in all water sources were further supported by a hierarchical cluster analysis (CA) dendrogram. The primary clusters found in WRS were (1) Mn-Zn-Ba-Ni-Cu-Fe, (2) Pb, and (3) As (Figure 13a). The primary clusters for GW were (1) Ba-Mn-Zn-Cu-Ni-As, (2) Fe, and (3) Pb, as shown in Figure 13b. Lastly, the primary clusters found in TW were (1) Ba-Mn-Zn-Cu-Fe-Ni, (2) Pb, and (3) As (Figure 13c).

Pollution and Health Risk Mapping Using the MLGI Method
The results of the simulation for the MPI mapping of the DW samples using the NN-PSO technique integrated with the EBK method are presented in Table 11. The validation and testing plots of the MPI for the DW samples are shown in Figure A7 in Appendix B. As shown in Figure 14a, the highest MPI for TW was observed in Brgy. Market Site, Municipality of Mogpog. The highest MPI was observed in Brgy. Bicas-Bicas, Municipality of Buenavista for the GW samples ( Figure 14b). Moreover, it was observed that the highest MPI for WRS was in Brgy. Janagdong, Municipality of Mogpog as shown in Figure 14c.

Pollution and Health Risk Mapping Using the MLGI Method
The results of the simulation for the MPI mapping of the DW samples using the NN-PSO technique integrated with the EBK method are presented in Table 11. The validation and testing plots of the MPI for the DW samples are shown in Figure A7 in Appendix B. As shown in Figure 14a, the highest MPI for TW was observed in Brgy. Market Site, Municipality of Mogpog. The highest MPI was observed in Brgy. Bicas-Bicas, Municipality of Buenavista for the GW samples ( Figure 14b). Moreover, it was observed that the highest MPI for WRS was in Brgy. Janagdong, Municipality of Mogpog as shown in Figure 14c.    Table 12 shows the simulation results for the NPI mapping of DW samples using the NN-PSO approach combined with the EBK method. The validation and testing plots for the NPI maps are illustrated in Figure A8 in Appendix B. As seen in Figure 15a, Brgy. Buangan, Torrijos reported the highest NPI for TW. The highest NPI values for GW were found in Brgy. Bicas-Bicas, Municipality of Buenavista as shown in Figure 15b. Additionally, the highest NPI for WRS was recorded in Brgy. Janagdong, Mogpog Municipality, as shown in Figure 15c.  The findings of a simulation for the HI (adults) mapping of domestic water samples using the NN-PSO approach with the EBK technique are shown in Table 13. The validation and testing plots for the HI (adults) maps are illustrated in Figure A9 in Appendix B. As shown in Figure 16a, Brgy. Buangan, Torrijos reported the highest HI (adults) for tap water. As seen in Figure 16b, the highest HI (adults) values for GW were observed in Brgy. Bicas-Bicas, Municipality of Buenavista. Additionally, the highest HI (adults) for WRS was recorded in Brgy. Janagdong, Mogpog Municipality as shown in Figure 16c.  The findings of a simulation for the HI (adults) mapping of domestic water samples using the NN-PSO approach with the EBK technique are shown in Table 13. The validation and testing plots for the HI (adults) maps are illustrated in Figure A9 in Appendix B. As shown in Figure 16a, Brgy. Buangan, Torrijos reported the highest HI (adults) for tap water. As seen in Figure 16b, the highest HI (adults) values for GW were observed in Brgy. Bicas-Bicas, Municipality of Buenavista. Additionally, the highest HI (adults) for WRS was recorded in Brgy. Janagdong, Mogpog Municipality as shown in Figure 16c.   Table 14 summarizes the results of a simulation for the HI (children) mapping of residential water samples using the NN-PSO methodology integrated with the EBK method. Figure A10 in Appendix B illustrates the validation and testing plots for the HI (children) maps. Brgy. Buangan, Torrijos recorded the highest HI (children) for tap water, as seen in Figure 17a. As seen in Figure 17b, the highest HI (children) values for GW were found in Brgy. Bicas-Bicas, Buenavista Municipality. Additionally, as shown in Figure 17c, the highest HI (children) for WRS was obtained in Brgy. Janagdong, Mogpog Municipality.  The findings of a simulation for the TCR (adults) mapping of domestic water samples using the NN-PSO approach with the EBK method are shown in Table 15. The validation and testing plots for the TCR (adults) maps are illustrated in Figure A11 in Appendix B. As shown in Figure 18a, Brgy. Buangan, Torrijos reported the highest TCR (adults) for TW. As illustrated in Figure 18b, the highest TCR (adults) levels for GW were found in Brgy. Bagacay, Buenavista. Additionally, the highest TCR (adults) for WRS was recorded in Brgy. Janagdong, Mogpog Municipality as shown in Figure 18c.    Table 16 summarizes the simulation results for the TCR (children) mapping of domestic water samples using the NN-PSO methodology integrated with the EBK method. Figure A12 in Appendix B illustrates the validation and testing plots for the TCR (children) maps. Brgy. Buangan, Torrijos reported the highest TCR (children) for TW, as seen in Figure 19a. As illustrated in Figure 19b, the highest TCR (children) levels for GW were found in Brgy. Bagacay, Buenavista. Additionally, the highest TCR (children) for WRS was recorded in Brgy. Janagdong, Mogpog Municipality as shown in Figure 19c.

Discussion
Identifying the parameters influencing the GW chemistry is vital for the sustainability of water resources [65]. In this study, it was observed that the mean concentrations of As, Pb, and Ni in domestic water resources in the province of Marinduque were above the limit set by the WHO. The relationship of water quality parameters was likewise determined, and it was recorded that there was a substantial positive correlation between EC and TDS in tap water samples, which agrees with the results of the study by Qureshi et al. [66]. Moreover, a high positive correlation was also seen between Pb and Cu in TW and Pb and Fe in WRS, which is similar to the findings of Varghese and Jaya [67]. For WRS, Fe and Cu had a significant positive correlation that agreed with the results of Kuisi and Abdel-Fattah [68]. A positive correlation among these metals suggests a possible shared origin [69].
The MPI and NPI were broadly utilized in the evaluation of the total HMM contamination in GW. The MPI covers the weight of different HMMs in the computation of the overall quality of drinking water [4,70]. In the present study, HMM including As, Ba, Cu, Fe, Pb, Mn, Ni, and Zn were considered. The findings showed that all sampling locations recorded MPI values classified as having a high degree of pollution. Similarly, all sampling points were heavily polluted based on the calculated NPI values. As, Pb, and Ni in all water sources exceeded both the water quality standards of WHO and PNSDW. Similar findings were observed by Agarin et al. [14], who investigated the concentration of metals in the surface and groundwater of the island province in 2021.
The total carcinogenic risk computed through the ingestion route exceeded the maximum threshold level of 1.00 × 10 −4 [71] for all DW samples. The results indicated that children were more vulnerable to CR than adults, as shown by the observation that the TCR values for children were all greater than the TCR values calculated for adults, consistent with the results of Pervez et al. [35]. Therefore, it is essential to determine the As levels because exposure may cause either acute or chronic poisoning. Acute As toxicity is often characterized by nausea, vomiting, abdominal discomfort, and severe watery diarrhea [72]. Chronic arsenic effects may arise due to prolonged exposure to lower arsenic levels, although latent toxicity, such as cancer, can persist even after exposure has ended. Chronic arsenic toxicity may gradually develop, making it more difficult to identify. Skin manifestations and peripheral neurologic complaints are often more apparent than gastrointestinal symptoms, with chronic exposure and concerns also centering on an increased future risk of cancer [73]. Early skin indicators include hyperpigmentation or hypopigmentation of the skin. Hyperkeratosis and scaling, especially on the palms and soles, highly indicate arsenic exposure [74]. Carcinomas of the skin and Bowen's disease are connected with latent arsenic toxicity. Additionally, a peripheral vascular condition known as a Blackfoot disease accompanying gangrene has been related to chronic arsenic exposure [75,76]. Numerous research and case reports have also shown a link between arsenic exposure and cancer. Arsenic exposure has been linked to skin [77], lung [78], liver [79], kidney [80], bladder [81], and prostate cancers [82].
The different pollution and health risk indices were mapped using the MLGI approach, which combines NN-PSO with EBK [21]. The governing MLGI models were found to be established from the AIC measure values that were the lowest among the observable hidden neurons. As the AIC value approached its minimal value, it was revealed that increasing the quantity of hidden neurons resulted in an increased AIC value. This indicated that the network had reached a state of generalization [58].
The concentrations of HMM do not necessarily reflect the actual pollution state of a water resource since it only assesses each heavy metal separately and with equal severity in terms of its biological effects. The application of pollution and health risk indices could provide a more conclusive assessment of the pollution and health risk status of a region since it compares the concentrations of each heavy metal to their acceptable levels [102]. It also indicates the total amount of pollution a region is facing. These indices were applied to an island province that was hit by two mine tailing disasters and it could be adopted in other locations across the world that experienced similar mining disasters.
The utilization of pollution and health risk indices is critical for ensuring the quality of water resources consumed and used by residents. These indices can be employed to determine the current risk to which the community may be exposed and recommend potential mitigation measures to limit the level of exposure of community members [103]. The incorporation of the MLGI approach enables the creation of spatial maps of various indices to provide remediation strategies and decision-making processes.

Conclusions
This study examined the spatial variability of the MPI, the NPI, the HI and the TCR (for children and adults) by HMM such as As, Ba, Cu, Fe, Pb, Mn, Ni, and Zn in an island province's DW sources, i.e., TW, GW, and WRS. Moreover, it highlighted the use of portable devices for in situ detection of HMMs in water, providing accurate and realtime results. It was recorded that As, Pb, and Ni concentrations were much higher than the WHO's permissible levels. The correlation analysis and dendrogram revealed that the link between the HMMs detected indicated a common origin. The computed MPI values in all sampling sites for all domestic water samples were characterized as highly polluted. Additionally, the NPI values estimated at all sampling stations for all DW samples were categorized as heavily polluted. The HQ indices for As and Pb were significantly greater than 1 for all water sources, indicating a potentially significant health risk to the human population. Moreover, As had the highest CR observed in the tap water samples and accounted for 89.7% of the total CR in TW. Additionally, all sample sites exceeded the USEPA's suggested maximum threshold level of 1.00 × 10 −4 , which indicates a high carcinogenic risk. The sensitivity analysis results using MCS showed that the most influential variable in the LTCR was the contaminant concentration. Further, reducing the metal concentration reduced the carcinogenic risk. The calculated indices and MLGI-based maps can be utilized as a benchmark for future research by local government units. It would be helpful in the creation and implementation of remediation and mitigation strategies to promote sustainable development in the protection of domestic water resources. Prompt intervention is required, as well as a regular monitoring of the quality of domestic water sources in the island province. This is to reduce the negative health impacts of these elevated HMM to the local population. Priority and special attention shall be given to the HMM that have elevated concentration and exceed the permissible limits by WHO and PNSDW. In addition, a regular monitoring of water resources quality for domestic supply shall be carried out, the development of policies based on these results shall be pursued, and the enforcement of existing policies shall be effectively conducted to reduce adverse health effects to the population of the island province.        Toxics 2022, 10, x FOR PEER REVIEW 12 of 16 Appendix A Figure A1. SA of LTCR model (adults) for As in water. Figure A2. SA of LTCR model (adults) for Ni in water. Figure A3. SA of LTCR model (adults) for Pb in water. Figure A3. SA of LTCR model (adults) for Pb in water.