Fuzzy-Probabilistic Model for a Risk Assessment of Groundwater Contamination: Application to an Urban Zone in the City of Belém, Pará, Brazil

This study proposes a fuzzy-probabilistic modelling approach for groundwater contamination risk assessment (FPM-risks) regarding underground fuel storage tanks (UFST). Considering the subjective measures of hydrogeological parameters, a fuzzy inference system is proposed to assess the intrinsic vulnerability of aquifers (Y). Measurement of the UFST hazard degree (H) and natural groundwater quality (Q) is considered as a pattern framing issue, such that they were quantified by fuzzy-analytic hierarchy process (AHP) of the recognizing patterns model. Though the association among Q and the probability of using groundwater reserves (G), estimated by the Monte Carlo method, the consequences of contamination (C) were measured. Associating Y, H, and C, the basic and value-weighted risk assessment of groundwater contamination was performed in the urban zone of Belém city, Pará state, Brazil. The results showed that the majority of UFSTs concentrated in the more urbanized zone were classified by FPM-risks as high basic risk and very high value-weighted risk of groundwater contamination. Although the risk assessment should be updated regularly because of the dynamic characteristics of hazards from the USFTs, the FPM-risks was shown as a tool to be considered for managing groundwater resources, as these models overcome subjectivities and address uncertainties, thereby providing a higher level of accuracy than usual risk methods and possibly become a decision-making way.


Introduction
Risk is usually defined by the association of the terms origin-pathway-target [1,2]. The "origin" refers to the identification of polluting agents and contaminant loads associated with the source hazards (H). The "pathway" is associated with the intrinsic vulnerability (Y) of an aquifer for contamination. For association among H and Y, the authors of [2] suggested a basic risk (R b ) to express it by the following: R b = Y × H. The "target" is related to the possible consequences of contamination events (socioeconomic and environmental impacts). According to the authors of [2], the consequence of contamination (C) is associated with the product of the groundwater reserve (G) and the quality (Q) of the groundwater (i.e., C = G × Q). Thus, a type of risk known as the value-weighted risk (R) is calculated by the following equation: R = R b × C.      Figure 2. FIS for assessing the intrinsic vulnerability degree by Y GOD-Fuzzy *. GOD-groundwater hydraulic confinement, overlaying strata, and depth to groundwater table.
Water 2020, 12, 1437 5 of 27 The fuzzy index of intrinsic vulnerability (Y GOD-Fuzzy *) is obtained through the interaction among the lithology (x O ) and HAI ( Figure 2).
There are nine rules to define the relationships among x O and HAI, resulting in Y GOD-Fuzzy * (e.g., Rule: "IF" the hydraulic accessibility is "L", and the contaminant attenuation capacity is "L", "THAN" the intrinsic vulnerability is "M"). The defuzzification is processed using the centroid method, by calculating the HAI and Y GOD-Fuzzy *.
The application of the FIS is preceded by a performance analysis and subsequent validation. The performance of the GOD-fuzzy model is analyzed via error minimization. Therefore, the lower the relative errors, the better the performance of the GOD-fuzzy model, which is calculated as follows: where RE refers to the relative error (%), "N" is the number of data points, Y GOD is the output value of the GOD method, and Y GOD-Fuzzy * is the output value of the FIS.
The GOD-fuzzy model has better performance when used to classify a larger degree of vulnerability relative to the GOD-based method. Thus, more preventative interpretations of vulnerability by the FIS are suggested. The GOD-fuzzy model is validated by inserting samples of the primary contaminant tracer from a contamination source on the vulnerability map. If the samples indicate higher-than-normal values and are in areas in which the degree of vulnerability does not correspond to these values, the rules of the FIS are adjusted. The vulnerability map is produced by the spatial interpolation of the Y GOD-Fuzzy * index.

Hazard Pattern Definitions
The contamination sources considered here are fuel stations (FS), which handle and store BTEX (benzene-boluene-ethylbenzene-xylene), benzo(a)pyrene, and naphthalene compounds. Therefore, this source fits into five hazard classes (very low "VL", low "L", moderate "M", high "H", and very high "VH"). The classes are functions of the five criteria associated with the physico-chemical properties of the contaminant, including its toxicity (p 1 ), mobility (p 2 ), degradability (p 3 ), solubility (p 4 ), and contaminant volume (p 5 ).
• p 1 : refers to the adverse effects of the interaction of the chemicals with the human body. Therefore, the lower the permitted concentration in potable water, the more toxic the action of the compound. It is expressed in mg/L. • p 2 : quantified by the water-octanol partition coefficient (K OC ) and associated with the mobility of a given compound and its displacement in the direction of the aquifer. Therefore, the lower log(K OC ), the more recalcitrant the compound tends to be, that is, there will be negligible adsorption to soils and sediments and rapid migration to groundwater. • p 3 : an index produced using the EPI Suite™ (Estimation Program Interface) application, which was developed by the authors of [31]. The lower the index, the greater the resistance of the compound to chemical and biochemical degradation. • p 4 : associated with the saturated concentration of a compound in water, conditioned to a given temperature and pressure, and inversely proportional to the K OC . Thus, p 4 indicates that compounds with elevated solubility are more susceptible to biodegradation and negligible sorption by soils and sediments, and rapid migration potential to groundwater. • p 5 : refers to the amount of contaminant stored by the underground fuel storage tank (UFST). The authors of [32] indicated that an UFST has a capacity of 15 m 3 /fuel station (FS). Thus, to define the p 1 , p 2 , p 3 , and p 4 standards, the compounds most frequently used in the FS are considered by using the prescription referred to in [2], and they are expressed in Table 1. Using Table 1, an apparent hazard pattern matrix was developed (Table 2).  However, some contaminants are of greater importance in defining the value of p n -properties. Thus, the average weighted value for each p n (p n(m) ) is calculated, referring to the main contaminants of the FS activity. The weights are obtained using an AHP method (analytic hierarchy process) [33]. The Appendix A shows the acquisition of the referenced weights. The value of p n(m) is thus determined (product of w c -the weight of contaminant c, and the value of p n ).
Using this procedure, it is possible, for example, to define that compounds such as BP and B are most important for the value p 1 , because they are more toxic than the other compounds.

Natural Groundwater Quality Pattern Definitions
To classify the groundwater quality within various types of hydrogeological domains (sedimentary, metasedimentary, karst, and crystalline types), the base method proposed in [24] was used. This method weighs the use of the pH (q 1 ), chloride (q 2 ), total residuals (q 3 ), hardness (q 4 ), nitrate (q 5 ), and fluorine (q 6 ) parameters, evaluating the water quality using a Q-index, which is qualified by degrees. These degrees are optimal "Op" (80 < Q ≤ 100), good "Go" (52 < Q ≤ 79), acceptable "Ac" (37 < Q ≤ 51), and unfit "Un" (0 ≤ Q ≤ 36). This base method is restricted to classifying the natural quality of groundwater for human consumption. However, classifying the groundwater quality is considered a multi-criteria pattern classification problem, and the weighting of parameters is an important issue, even during the definition of the water use type. For potable purposes, the weighting of q n -parameters is preceded by a qualitative analysis, related to the natural origins of the impurities that can place the potability of groundwater at risk and contribute to the emergence of illnesses as follows: • q 1 : at pH~5.5, meteoric water interacts with iron sulfides in mining areas, and its acidity increases. More acidic waters attack most rocks and minerals, increasing the percentage of heavy metals with toxic potential in solution and creating toxic conditions [34].
• q 2 : during water treatment by chlorination and depending on the water quality, potentially carcinogenic chloride can be formed. The World Health Organization (WHO) recommends chloride concentrations of <250 mg/L. • q 3 : to classify the natural groundwater quality, a weight is conferred that is also referenced in the proposal in [35]. • q 4 : when ingesting hard and very hard waters, renal and cardiovascular diseases can develop [34]. CaCO 3 concentrations <50 mg/L are soft, 50 to 100 mg/L are slightly hard, 100 to 200 mg/L are hard, and >200 mg/L are very hard. • q 5 : although NO 3 at low concentrations is slightly toxic to humans, this substance can be reduced in the body (the stomach, intestines, or liver) into nitrite ions (NO 2 − ), causing the emergence of illnesses. In addition to methemoglobinemia, the authors of [34] noted that NO 2 − in the human body can react with amine substances, resulting in nitrosamine (carcinogenic agent). The maximum permitted concentration of nitrogen-nitrate (N-NO 3 ) in groundwater for human consumption was set by the WHO at ≤50 mg of N-NO 3 /L. • q 6 : fluorine concentrations >1.5 mg/L can cause dental fluorosis. Concentrations of ≥0.5 mg/L to ≤1.5 mg/L act as a natural tooth protectant, aiding in the avoidance of cavities. Concentrations of ≥4 mg/L can cause skeletal fluorosis (pains in the back and neck bones). Concentrations of >10 mg/L are associated with the appearance of deforming fluorosis, which causes bone hyperdensity, leading to disability, and excessive doses can lead to death [34].
Thus, as a function of these qualitative analyses, q 5 , q 2 , q 6 , q 1 , q 4 , and q 3 make up the hierarchy of parameters for classifying the natural quality of groundwater for potable uses, as evaluated by AHP (Appendix A). Thus, a multi-criteria model is structured when defining the parameters (criteria) and water uses (alternatives). Therefore, a multi-criteria fuzzy model is adopted for pattern recognition [36]. Initially, the relative pattern matrix is defined with regard to the potable water uses, which is conceived as a function of q n -parameter. It is categorized by four "h" classes relative to the degrees: 4 "Op", 3 "Go", 2 "Ac", or 1 "Un". Thus, the Z q matrix is defined, and its terms are z qn,h (Table 3). Table 3. Natural groundwater quality pattern matrix (Z q ). Un-unfit; Ac-acceptable; Go-good; Op-optimal.

Parameters-q n
Natural Groundwater Quality Degree-Q n 1 "Un" 2 "Ac" 3 "Go" 4 "Op" q 1 <6->9 6-9 7-8 7.5 q 2 (mg/L) The patterns indicated by the Brazilian Ministry of Health, U.S. EPA (United State of Environmental Protection Agency), National Council on the Environment (Conselho Nacional do Meio Ambiente-CONAMA), and World Health Organization (WHO) for determining the natural water quality were used. For values below the maximum permitted values, the designations "Op" and "Go" are given. For values equal to the maximum permitted, "Ac" is given, and those above the permitted maximum are designated "Un".

Fuzzy-AHP Multi-Criteria Model for Recognizing Patterns
The multi-criteria model for recognizing the fuzzy patterns proposed in [36] was used to fit the dataset to the previously established hazard degree and groundwater quality patterns. Both the data and the patterns must be placed in matrices. Therefore, the number of contaminating sources, associated p n , and water samples related to q n must be fit to the patterns defined by the matrices Z p ( Table 2) and Z q (Table 3).
To perform the fitting, the terms of the data matrixes X p and X q , and x p,j and x q,j , respectively, must be made uniform by normalization or harmonization, thus becoming x p,j ' and x q,j ', respectively. When the hazard and quality increase with increases in p and q, the normalization proceeds and harmonization decreases. This process proceeds in the same manner with the terms of the pattern matrices z p,h and z q,h . Thus, if z p,h ' and z q,h ' are equal to 1, a source or water sample belongs completely to class h, while if z p,h ' and z q,h ' are equal to 0, it is completely outside of class h.
The recognition model assumes the degree of belonging of each sample j to class h, considering a matrix U whose terms are designated by u h,j . To classify the hazard degree or the water quality, the value of h can be 1-5 or 1-4, considering the existence of five and four h-classes, respectively. The authors of [36] stated that, according to the fuzzy sets theory ( [37]), U is subject to the following restriction: Given this restriction, p n and q n are important for classifying their respective degrees. The fuzzy pattern recognition model considers a distance that separates j from its degree h equal to the following: in which w (p,q) is the weight of w pn and w qn that influences the hazard assessment and water quality classification, respectively. The attribution of weights was measured by AHP ( [33]) (see Appendix A). After w pn and w qn are defined from the perspective of the fuzzy theory, the term u h,j is considered as a weight for d h,j . Thus, a weighted distance (D h,j ) is defined; it separates j from its class h and is equal to u h,j *d h,j . To solve u h,j , the following objective function must be minimized: To obtain a Lagrangian function and optimize the classification process by restriction, in which λ j is a Lagrange multiplier. The resolution of the optimal value that fits a source or sample j to the degree of pertinence h by u h,j meets two conditions, namely, Using the partial derivatives of the Lagrange function, the pertinence function of u h,j is obtained for class h by using the following: Equation (8) is conditioned to d h,j 0. When d h,j = 0, the j sample belongs to level h, and thus u h,j is equal to 1. However, in using Equation (8), a classification can be performed only by the extremes of the h classes (i.e., the data matrixes X p ' and X q ', which could be grouped into two levels of h classes corresponding to the pattern matrices Z p ' and Z q ', that is, u h,j = 0 or u h,j = 1). The data matrices X p ' or X q ' can be fit to the interval between 0 and 1, such that they can correspond to the categories of h classes and are associated with their respective values (from category 1 to 5 and from 1 to 4). Therefore, u h,j *, as terms of the matrix U*, which represents a matrix of pertinence of each sample j belonging to each class h, will vary from 1 to 5 and from 1 to 4. Classification into one of the h categories cannot be effected directly by matrix U*, and thus the u j * terms that are equal to (u 1,j , u 2,j , u 3,j , u 4,j , and u 5,j ) are considered. Therefore, when vector u j * is associated with category h, a hazard index is obtained, H j(P) , as an index of the natural quality of groundwater for potable uses, H j(Q) .
where H j(P,Q) refers to the indices H j(P) and H j(Q) . To classify water quality, h varies from 1 to 4.

Probabilistic Estimate of Groundwater Reserve Use
The total groundwater reserve (G) can be calculated by summing up the G P and G R reserves. The G P reserve corresponds to the volume of groundwater in the saturated zone, which is below the minimum position of the seasonal variation of the aquifer. The G R reserve corresponds to the infiltrated water volume or annual recharging of the aquifers. The estimate for groundwater reserves meeting the demand refers to the probability that the annual consumption for multiple uses is met such that it does not present a risk to the G R reserves, and the G P reserves are conserved. However, the lack of data for calculating the variables that make up the calculation of G P and G R is an aspect that denotes uncertainties.
A strategy that considers the random nature of the variables η e and ∆h is used with the adoption of a probabilistic method. The Monte Carlo simulation (MCS) method requires a random sampling method, in addition to information about the shapes of the frequency distributions of a variable and its statistical attributes (e.g., its mean and standard deviation). Simple random sampling (SRS) is used to simulate η e and ∆h. In this research, 10,000 iterations are considered. The random value is sampled from each distribution of the random variable under study. The probability density function of the regulation of groundwater reserves, f (G R ), is determined as follows: (i) generation of pseudo-random numbers for η e and ∆h and (ii) acquisition of f (G R ) by the following: The term A (area of occurrence of the aquifer) was considered to be a deterministic variable. To generate pseudo-random numbers, the supplement NtRand ® was used, which was installed in Microsoft Excel ® (version 2016, development by Microsoft Corporation, Redmond, Washington, USA). A conservative concept was adopted in which the annual average G R must be fully exploited, given that part of this reserve must ensure river flows. A number W of inhabitants is considered that requires a quantity of water such that it is possible to estimate the average annual consumption.
If this volume of water is removed from G R (it can be used without harming the aquifer), there is a probability that the exploited volume is less than that required for annual consumption. Therefore, the larger this probability is, the greater the chances are that G R cannot be used safely, thus requiring the use of G P . In these terms, the following was adopted when the probability percentage is as given below: • ≤10%: very low probability "VL" of G R being used without the need to use G P . • >10% and ≤30%: low probability "L" of G R being used such that part of G P can be used. • >30% and ≤50%: moderate probability "M" of G R being used, with partial use of G P . • >50% and ≤70%: high probability "H" of G R being used, requiring the use of G P . • >70%: very high probability "VH" of G R being used, requiring a large use of G P .

Case Study: Belém Urban Zone
The application of the fuzzy-probabilistic approach was conditioned by data acquisition strategies for a domain that corresponds to the urban zone of Belém city, in the northern region of Brazil, which is an integral part of the metropolitan region of Belém (MRB). In addition to belonging to the hydrogeological province of the Amazon River, which is one of the largest in Brazil with a volume of greater than 30 thousand km 3 , the selection of this spatial domain is justified by the fact that evidence of hydrocarbon contamination from leaking UFST has been found.
To apply the risk model, the geological and hydrogeological data were initially extracted to evaluate the intrinsic vulnerability. Thus, data were extracted from the Groundwater Information System (Sistema de Informação de Águas Subterrâneas-SIAGAS), which was generated by the Research and Mineral Resources Company (Companhia de Pesquisa e Recursos Minerais-CPRM, Serviço Geológico do Brasil) belonging to the Brazilian Geological Service. A total of 762 data points from wells were catalogued. The information referred to the hydrostatic level and the lithological and construction profiles of the unsaturated zone (to a depth of less than 2 m), the spatial distribution of which is illustrated in Figure 3A,B.
The vulnerability of the referenced urban zone was measured in shallow aquifers in which a large portion of contamination events began in the most superficial layers of the underground environment. The authors of [38][39][40] discussed the hydrogeological aspects of the metropolitan region of Belém (MRB).
In this vulnerable environment, UFST leaks have already occurred. In northern Brazil, there are currently 2798 registered fuel stations. Their average monthly sales volume as recorded through the end of 2013 was greater than 160 m 3 /station. Data from the National Petroleum, Natural Gas, and Biofuel Agency of Brazil (Agência Nacional do Petróleo, Gás Natural e Biocombustível do Brasil-ANP) show that, in the MRB, there are more than 200 registered fuel stations. In the urban zone of the city of Belém and its islands, there were 91 FS, and their storage capacity is illustrated in Figure 4. Volumetric data are available on the ANP platform.
With regard to groundwater quality evaluations, recent studies have restricted the analysis to tracers of anthropogenic contamination (e.g., pH, NO 3 , and fecal coliforms) [41][42][43]. This information is pertinent because the shallow aquifers are more exposed to contamination than confined aquifers (which are "naturally" protected). However, the natural quality of groundwater evaluated in studies in [44] showed that, in shallow wells (depths ≤ 50 m), the pH varied from 3.4 to 5.2 and directly corresponded to low levels of electrical conductivity, hardness, and bicarbonate alkalinity, classifying them as chlorinated sodium, sulfated calcium, and sodium bicarbonate when using the Piper diagram.
The samples made available by the authors of [44] showed the physicochemical parameters that make up the natural quality analyses of groundwater, which can also be tabulated. In turn, 14 data points were inventoried from water analyses of wells (from shallow aquifer) in the urban and periurban zones of the city of Belém and its islands, extracted from the work of [44]. The authors of [41] analyzed 20 samples from wells, which captured groundwater from the phreatic aquifer, in June 2002. Nevertheless, data from an experimental water sampling campaign were collected from 13 wells in [42] from June 2000 to March 2002, and the average values were used. Data catalogued in [43] were also used, and they were from water collected from 24 wells in the urban and periurban portion of the island of Mosqueiro during two seasons, namely, the dry season in October 2010 and the rainy season in March 2011 (average values were used). Owing to the lack of data from the island of Outeiro, the natural quality of the groundwater there was analyzed as being equal to the portion that corresponds to the urban zone of the island of Mosqueiro. Figure 4 shows all of the sampling used.  The vulnerability of the referenced urban zone was measured in shallow aquifers in which a large portion of contamination events began in the most superficial layers of the underground environment. The authors of [38][39][40] discussed the hydrogeological aspects of the metropolitan region of Belém (MRB).
In this vulnerable environment, UFST leaks have already occurred. In northern Brazil, there are currently 2798 registered fuel stations. Their average monthly sales volume as recorded through the end of 2013 was greater than 160 m 3 /station. Data from the National Petroleum, Natural Gas, and Biofuel Agency of Brazil (Agência Nacional do Petróleo, Gás Natural e Biocombustível do Brasil-ANP) show that, in the MRB, there are more than 200 registered fuel stations. In the urban zone of the city The demand for groundwater was characterized in a study conducted by [38], which indicated a water quantity of slightly greater than 430,000 m 3 /day, of which approximately 20% comes from aquifers (around 86 thousand m 3 /day). The growing urban demand for treated water for industrial and domestic uses combined with the current inefficiency and insufficiency of the water supply network in the urban and periurban zones of Belém and its islands is reflected in the increased capture of groundwater from phreatic wells. A per capita demand of 200 liters/day was adopted, as projected for the year 2033 [45], considering that 50% of the potable and non-potable demand by the city of Belém comes from wells, that is, for approximately 750,000 inhabitants (the population of the city of Belém is close to 1.5 million inhabitants, according to data in [46]), the daily volume is close to 150,000 m 3 (54.75 million m 3 /year). of Belém and its islands, there were 91 FS, and their storage capacity is illustrated in Figure 4. Volumetric data are available on the ANP platform. With regard to groundwater quality evaluations, recent studies have restricted the analysis to tracers of anthropogenic contamination (e.g., pH, NO3, and fecal coliforms) [41][42][43]. This information is pertinent because the shallow aquifers are more exposed to contamination than confined aquifers (which are "naturally" protected). However, the natural quality of groundwater evaluated in studies in [44] showed that, in shallow wells (depths ≤ 50 m), the pH varied from 3.4 to 5.2 and directly corresponded to low levels of electrical conductivity, hardness, and bicarbonate alkalinity, classifying them as chlorinated sodium, sulfated calcium, and sodium bicarbonate when using the Piper diagram.
The samples made available by the authors of [44] showed the physicochemical parameters that make up the natural quality analyses of groundwater, which can also be tabulated. In turn, 14 data points were inventoried from water analyses of wells (from shallow aquifer) in the urban and periurban zones of the city of Belém and its islands, extracted from the work of [44]. The authors of [41] analyzed 20 samples from wells, which captured groundwater from the phreatic aquifer, in June  The island zones are dependent on capturing groundwater to supply demand for multiple uses. With an estimated population of close to 50,000 inhabitants, the island of Mosqueiro demands a volume of close to 10,000 m 3 /day (3.65 million m 3 /year). For the island of Outeiro, with a population of close to 90,000 inhabitants, the daily demand is close to 18,000 m 3 (6.57 million m 3 /year). Therefore, from these annual demands for potable and non-potable uses, the probability that the G R reserve will be used safely is estimated, without the need to use the G P reserve of the confined aquifers.

Assessment of Fuzzy Intrinsic Vulnerability
Although the GOD-fuzzy model shows a relative error (RE) of 10.92% compared with the GOD method, the FIS was used to evaluate the intrinsic vulnerability, because the GOD method (i) could lead to vague or diffuse conclusions and (ii) has no alternative to classifying the intrinsic vulnerability in an Aristotelian way, for example, the lithology of the aquifer is either sandy-clay (true) or not (false). These interpretations are not considered by the fuzzy principle model, as the lithology of the aquifer, for example, can be more or less clay-sandy and then become clay, sandy, or even sandy-clayey. Thus, using a non-fuzzy approach could compromise the assessment of the aquifer intrinsic vulnerability. Therefore, an interpolation of the Y GOD-Fuzzy * index was performed, such that your spatial distribution is isotropic and was adjusted by a power model (POW), identifying the nugget effect "C 0 " and the constant term "α" and power "p" (Figure 5A) (the details of which are described in Appendix B). The fuzzy vulnerability map in a large portion of the spatial domain was classified as moderate, or "M", in degree. The mapping of Outeiro Island indicated that more than 50% of the urban and periurban zone has high "H" intrinsic vulnerability to contamination and that slightly more than 40% of the referenced zone was classified as having moderate "M" intrinsic vulnerability ( Figure 5B). Figure 5C shows the high "H" intrinsic vulnerability coverage of the shallow aquifer that surfaces in the urban and periurban zone of the island of Mosqueiro, which was less than 4%, compared with the 90% classified as "M". Less than 30% of the urban and periurban zone of Belém was classified as having high "H" intrinsic vulnerability for the contamination of the shallow aquifer. Slightly more than 70% of this zone was classified as having moderate "M" intrinsic vulnerability ( Figure 5D).
The results obtained by the interpolation of the Y GOD-Fuzzy * index partially aligned with the results obtained by the authors of [40]. Lastly, the fuzzy intrinsic vulnerability map was validated ( Figure 5E). For this purpose, well water samples were included, which, after analysis, demonstrated contamination by BTEX [46] in a zone classified by the GOD-fuzzy method as having high intrinsic vulnerability.

Fuzzy-AHP Classification of the Hazard Degree of UFST
Before the application, the configured fuzzy-AHP model was measured to group the fuel stations by patterned hazard degree. Figure 5F shows the results of physicochemical properties weighted by the AHP method. Because the volumetric capacity of the UFST was the only property that varied between the stations, a linear simple regression was fitted for this measurement, which can be used to classify the hazard ( Figure 5G). The application of the fuzzy-AHP model quantified the hazard or potential to cause damage to groundwater from 91 FS (see Figure 4). The fuzzy-AHP model classified four UFST/FS as a very high "VH" hazard. However, 60 UFST/FS were classified as having moderate risk "M", and more than 29% were classified as having high risk "H" ( Figure 5H). The stations catalogued in the spatial domain that showed the most potential to cause damage are located in the city of Belém. Only one station was catalogued on the island of Mosqueiro, which was classified as having moderate "M" hazard ( Figure 5I).

Fuzzy-AHP Classification of the Natural Quality of Groundwater
The fuzzy-AHP model was configured to classify the natural quality of groundwater for potable uses. Figure 6A shows the results of physicochemical parameters weighted by the AHP method. Figure 6D shows the mapping of the natural quality of groundwater in the urban zone of Belém city by geostatistical interpolation of the H j(Q) -index, adjusted for exponential model (EXP), identifying the "C 0 " and sill (C) ( Figure 6B) (detailed in Appendix B). This way, 81% of the mapped area was framework acceptable "Ac" quality, and the remaining 19% was framework good "Go" quality ( Figure 6C).

Probability of Groundwater Reserve Use
For the variable ηe, the data obtained by the authors of [38] showed a probable value of 10%. Given the lack of data, a triangular distribution was adopted by considering an empirical criterion of ±10% of the probable value, showing a minimum of 0.05% and a maximum of 15%. Likewise, the lack of a monitoring network for the water column of wells justified the adoption of a probable value for ∆h equal to 1.8 m, as referenced in [38] for aquifers that surface in the city of Belém.
Thus, a triangular distribution was also adopted, with the minimum value under consideration of 0.8 m and the maximum value of 2.6 m, which are justifiable variations for Amazon water levels. On the basis of these considerations, pseudo-random numbers were generated for the physical variables η e and ∆h by simple sampling, which resulted in the triangular probability distribution functions (t-PDFs) illustrated in Figure 6E,F, respectively.     Thus, considering a surfacing area of the shallow aquifer systems of approximately 1200 km 2 , which corresponds to the saturation water [38,44], combined with the PDF of the variables η e and ∆h, f (G R ) was obtained. In general, through the PDF of the G R reserve, f (G R ), the probability of safely using these reserves without the need to use the G P reserves was obtained. The following approximations were adopted:

GR (m³)
• A per capita water demand of 200 L/day, as projected in the spatial domain for the year 2033 and referenced in [45].

•
Fifty percent of the consumption in the city of Belém comes from wells, that is, approximately 750,000 inhabitants demand a daily volume of 150,000 m 3 (54.75 million m 3 /year).

•
On the island of Mosqueiro, which depends on groundwater to supply an estimated population of approximately 50,000 inhabitants, the water volume demand is close to 12,500 m 3 /day (4.5 million m 3 /year). On the island of Outeiro, which also depends on groundwater, the population of close to 90,000 inhabitants demands approximately 22,500 m 3 per day (8.2 million m 3 /year).
According to these considerations, the probability of using the regulating reserve (G R ) in the city of Belém was estimated at almost 11%. Thus, the probability that G P would be used was classified as low "L", so the G R reserve can be safely used. On the island zones, the probability that the G P would be used was classified as very low "VL" (5.09% and 5.45%), so the G R reserve can be used with complete safety.
Although the Municipal Plan for Basic Sanitation, Water Supply, and Sanitary Sewer of Belém-Pará [45] projected a time horizon of 20 years to universalize the coverage of the water supply network, there are no studies suggesting that the demand for groundwater will decrease or increase. Nevertheless, the more than 5000 wells that have been drilled for potable and non-potable uses will not abruptly cease to operate. Thus, considering the low recharge of the shallow aquifer, an accumulated volume of five years in the city of Belém, island of Mosqueiro, and island of Outeiro was projected such that the demand increased to 273.75, 22.5, and 41 million m 3 , respectively.
In this scenario, the probability of using the G R reserve in Belém was estimated at slightly more than 61%. Thus, the probability that the G P reserve would be used was classified as high "H", so the use of the G R reserve was considered as unsafe. However, in the island zones, the very low classification "VL" for the probability of G P being used (6.94% and 9.15%) was not changed, meaning that G R can be safely used ( Figure 6G). Thus, to assess the consequence of contamination, the relationship among the natural quality degree of groundwater and probability of use groundwater was adopted ( Figure 6I).

Contamination Risk Assessment
The basic contamination risk was assessed (R b ) from the point sources located in the urban zone of the city of Belém and its islands using a risk matrix. This matrix relates Y GOD-Fuzzy * and H j(P) , and the product is R b . The results from this relationship were obtained by inputting the hazard quantification for each UFST/FS into the vulnerability map. The degrees of basic groundwater contamination risk by point sources are indicated in Figure 7A.
To assess the value-weighted risk (R), the degree of consequence from groundwater contamination was initially quantified by finding a relationship among the degrees of natural groundwater quality and the probability of G P of the phreatic aquifer being used. The contamination of optimal-quality water, for which the probability of the G R reserve being used is very high with the activation of a small portion of the G P reserve, shows the consequence of contamination to be very high or "VH". The contamination of groundwater with an unfit natural quality extracted from phreatic aquifers, for which the probability of the G R being used is very low, thus requiring a large portion of the G P reserve, has a very low or "VL" consequence.  In a large portion of the urban zone of Belém, the impact of contamination in zones classified as having good natural water quality was considered very high or "VH". In the urban and periurban zones of the islands of Mosqueiro and Outeiro, the contamination impact was classified as moderate or "M". Thus, the degrees of value-weighted risk relating to the contamination of the groundwater from point sources are indicated, and the results are distributed in the urban zone of Belém ( Figure 7B).

Conclusions and Recommendations
The groundwater contamination risk assessment method proposed here is a preventive alternative to be considered during the management of groundwater resources, as it overcomes subjectivities and addresses uncertainties in the assessment of intrinsic vulnerability, risk, and water quality and the estimation of groundwater reserve use. The GOD-fuzzy model was proposed to overcome subjective and ambiguous limitations, which are generally associated with the GOD method, and it was used for an intrinsic vulnerability assessment. The GOD-fuzzy model was processed by employing a systematic and simultaneous relationship between "fuzzifier" hydrogeological parameters, which was shown as an alternative with a more preventative nature than the GOD method through an intelligent controller.
For specific hydrogeological configurations, the knowledge base of the GOD-fuzzy model can be readapted such that the interval for the ranges of values for the linguistic variables that defined the fuzzifier of the hydrogeological parameters under consideration can be changed. Thus, the set of rules can also be adjusted, functioning as a type of calibrator for the referenced fuzzy model. In turn, depending on the type of aquifer, other base methods can be used, for example, SINTACS or AVI methods, to measure the fuzzy intrinsic vulnerability.
Within the scope of applying the Mamdani-type fuzzy controller, the use of "defuzzifier" methods other than the centroid method used in this study is suggested. The use of other "defuzzifier" methods can cause significant changes in the final results, and consequently, in the intrinsic vulnerability map. Therefore, the Y GOD-Fuzzy * index produced by the GOD-fuzzy model was interpolated in the environment within the urban zone of Belém, and the shallow aquifer intrinsic vulnerability map was produced.
The point sources were located in this environment, namely, fuel stations (associated with the UFST), whose degrees of hazard were quantified with an H j(P) index generated by a multi-criteria fuzzy pattern recognition model integrated with AHP. The judgment of the importance of the physical and chemical properties (toxicity, degradability, mobility, solubility, and quantity) considered to quantify the degree of hazard can be evaluated. Thus, the classification of a specific degree of hazard can differ from that obtained in this study. Likewise, the inclusion of other properties, for example, the half-life, can lead to hazard classifications different from those obtained in this study. The fuzzy-AHP model designed to classify the hazard is flexible to the point that it performs the following: • adapts to types of patterns established by any regulation; • performs convenient judgments of importance for each type of contamination source, including for scenarios with vehicles; • inserts n physical and chemical contaminant properties.
To quantify the groundwater consequence of contamination, a concept was used that measures the socioeconomic and environmental value of groundwater, which is a function of the association between the degrees of natural water quality and the estimated use of the aquifer reserves. For this purpose, the natural quality of the groundwater for potable uses was classified using a multi-criteria fuzzy-AHP model. This model was configured to classify the natural quality of groundwater from the weighting of the pH, chloride, total residue, hardness, nitrate, and fluorine parameters in the same manner as mentioned previously.
Considering the uncertainties and stochastic nature of the physical variables of an aquifer, a probabilistic proposal was adopted using the MCS method. Thus, the degree of probability of using the G R reserves without activating the G P reserves was estimated. These results can undergo changes if a less conservative concept is adopted, namely, the probability of using the G P reserves while also requiring the G R reserves. For future studies, the inclusion of a stochastic assessment of other physical variables is recommended, as is the use of other algorithms for the generation of pseudo-random numbers, for example, a Latin hypercube.
Finally, using the proposed FPM-risks, strategies can be developed to manage contamination sources to protect underground springs when considering a contamination event in a dense urban environment. The implementation of this model must not be the sole strategy for groundwater protection. Thus, strategies that can be integrated and coordinated with the characteristics of the proposed risk model must be adopted, such as the following: • Imposing compliance with special requirements for handling chemical and toxic substances and persistent wastewater in any industrial zone located in zones of low or moderate vulnerability, which includes the adoption of devices that can prevent and monitor flows; • Widespread weighted risk management must be implemented.
Thus, the fuzzy and probabilistic models proposed for risk assessment permitted the establishment of a connection between groundwater users and decision-makers through the following: • The synthesis and simplification of the hydrogeological conditions by the GOD-fuzzy model; • The quantification of the hazard posed by contamination sources and the classification of the natural groundwater quality for potable uses using the fuzzy-AHP model; • The probabilistic estimation of the use of G R reserves, without using G P reserves; • The generation of thematic maps, in particular.
Although, technically, the risk assessment should be updated regularly because of the dynamic characteristics of hazards from the USFT, the FPM-risks was shown as a tool to be considered for managing groundwater resources. On the order hand, as these models overcome subjectivities and address uncertainties, thereby providing a higher level of accuracy than usual risk methods, a better way for decision making became available. In practice, this study is presented as a strategy to expand the public's consciousness and understanding about the need to protect groundwater.

Conflicts of Interest:
The authors declare no conflict of interest.
Condition. C R ≤ 20% (is consistent with the judgment of importance between the criteria)

Appendix B. (Spatial Interpolation by Geostatistical Method)
To Map the Intrinsic Vulnerability Fuzzy.

. Preliminary Analysis of Drift Formation
Although the Y GOD-Fuzzy index presented a certain slope of the regression line, it was considered that this slope is not accentuated to the point of not varying in average in a well-defined behavior, in accordance with the latitude and longitude. Therefore, it was assumed that the respective indices did not present a drift phenomenon.
To Map the Intrinsic Vulnerability Fuzzy.

Preliminary Analysis of Drift Formation
Although the YGOD-Fuzzy index presented a certain slope of the regression line, it was considered that this slope is not accentuated to the point of not varying in average in a well-defined behavior, in accordance with the latitude and longitude. Therefore, it was assumed that the respective indices did not present a drift phenomenon.  The GOD-Fuzzy index estimates that assessing the intrinsic vulnerability of the urban zone of Belém and islands is valid and consistent without the need for corrective (e.g., trend kriging) filtering. It should be noted that the preliminary analysis, statistical and spatial description, pointed to a situation of homogeneity (only one population).

Appendix B.2. Variography
After several empirical attempts, through a preliminary approach of the geostatistical proposal, lag was adopted in the x and y direction, equal to 700 and 1600 m, respectively; and a lag number equal to 15, obtaining the variographic surface. Due to these conditions, a situation of slight anisotropy in the N-W and S-E directions occurred, because within the search circle, the pixels of blue tones below the variance were predominant ( Figure A2A). The GOD-Fuzzy index estimates that assessing the intrinsic vulnerability of the urban zone of Belém and islands is valid and consistent without the need for corrective (e.g., trend kriging) filtering. It should be noted that the preliminary analysis, statistical and spatial description, pointed to a situation of homogeneity (only one population).

Variography
After several empirical attempts, through a preliminary approach of the geostatistical proposal, lag was adopted in the x and y direction, equal to 700 and 1600 m, respectively; and a lag number equal to 15, obtaining the variographic surface. Due to these conditions, a situation of slight anisotropy in the N-W and S-E directions occurred, because within the search circle, the pixels of blue tones below the variance were predominant ( Figure B2A). On the other hand, an analysis of the variable surface suggests an isotropic tendency that can be modeled via the omnidirectional variogram model. To confirm this finding, after several empirical attempts, the tolerance, lag, and maximum bandwidth equaled to 550, 1950, and 650 m, respectively.
The omnidirectional variogram confirmed the nonexistence of the trend of YGOD-Fuzzy data. In these terms, it has been observed that the graph illustrated by Figure B2B points to a pure "white noise" nugget effect so that, in this case, the spatial correlation is tenuous, so as to call into question the advantage of a geostatistical approach to map the intrinsic vulnerability. Thus, the distribution of the isotropic YGOD-Fuzzy index was considered, adjusting the omnidirectional variogram via the power (POW) theoretical model ( Figure B3). On the other hand, an analysis of the variable surface suggests an isotropic tendency that can be modeled via the omnidirectional variogram model. To confirm this finding, after several empirical attempts, the tolerance, lag, and maximum bandwidth equaled to 550, 1950, and 650 m, respectively.
The omnidirectional variogram confirmed the nonexistence of the trend of Y GOD-Fuzzy data. In these terms, it has been observed that the graph illustrated by Figure A2B points to a pure "white noise" nugget effect so that, in this case, the spatial correlation is tenuous, so as to call into question the advantage of a geostatistical approach to map the intrinsic vulnerability. Thus, the distribution of the isotropic Y GOD-Fuzzy index was considered, adjusting the omnidirectional variogram via the power (POW) theoretical model ( Figure A3).

zone.
On the other hand, an analysis of the variable surface suggests an isotropic tendency that can be modeled via the omnidirectional variogram model. To confirm this finding, after several empirical attempts, the tolerance, lag, and maximum bandwidth equaled to 550, 1950, and 650 m, respectively.
The omnidirectional variogram confirmed the nonexistence of the trend of YGOD-Fuzzy data. In these terms, it has been observed that the graph illustrated by Figure B2B points to a pure "white noise" nugget effect so that, in this case, the spatial correlation is tenuous, so as to call into question the advantage of a geostatistical approach to map the intrinsic vulnerability. Thus, the distribution of the isotropic YGOD-Fuzzy index was considered, adjusting the omnidirectional variogram via the power (POW) theoretical model ( Figure B3). Figure B3. Omnidirectional and directional variogram-Belém urban zone. Therefore, it was not possible to make adjustments to the theoretical models with sill, because it is not possible to gauge the priority of a better fit for the points relative to the smallest distances (the Figure A3. Omnidirectional and directional variogram-Belém urban zone. Therefore, it was not possible to make adjustments to the theoretical models with sill, because it is not possible to gauge the priority of a better fit for the points relative to the smallest distances (the first points). Thus, the application of the linear model POW, not sill, was the one that best fit the pairs of points of the omnidirectional variogram.
Thus, to compose the studies on an intrinsic vulnerability analysis, the map obtained via spatial interpolation of the Y GOD-Fuzzy index, adjusted by the POW model, was adopted. Thus, cross-validation of the variographic model revealed a root mean square error (RMSE) of around 0.086414.
To Map the Groundwater Nature Quality fuzzy.

Appendix B.3. Preliminary Analysis of Drift Formation
Although the H j(Q) index presented a certain slope of the regression line, it was considered that this slope is not accentuated to the point of not varying in average in a well-defined behavior, in accordance with the latitude and longitude. Therefore, it was assumed that the respective indices did not present a drift phenomenon ( Figure A4).
Water 2020, 11, x FOR PEER REVIEW 23 of 27 first points). Thus, the application of the linear model POW, not sill, was the one that best fit the pairs of points of the omnidirectional variogram. Thus, to compose the studies on an intrinsic vulnerability analysis, the map obtained via spatial interpolation of the YGOD-Fuzzy index, adjusted by the POW model, was adopted. Thus, cross-validation of the variographic model revealed a root mean square error (RMSE) of around 0.086414.
To Map the Groundwater Nature Quality fuzzy.

Preliminary Analysis of Drift Formation
Although the Hj(Q) index presented a certain slope of the regression line, it was considered that this slope is not accentuated to the point of not varying in average in a well-defined behavior, in accordance with the latitude and longitude. Therefore, it was assumed that the respective indices did not present a drift phenomenon ( Figure B4).

Variography
In order to measure the spatial distribution of the fuzzy index Hj(Q) via a geostatistical approach, after several empirical attempts, a lag in the x and y direction was adopted at 8000 and 5600 m, respectively; and a lag number equal to 5, obtaining the variographic map. Under these conditions,

Appendix B.4. Variography
In order to measure the spatial distribution of the fuzzy index H j(Q) via a geostatistical approach, after several empirical attempts, a lag in the x and y direction was adopted at 8000 and 5600 m, respectively; and a lag number equal to 5, obtaining the variographic map. Under these conditions, little anisotropy in the N-S and E-W directions was observed; therefore, within a small search radius, the pixels of blue tones whose hue was below the variance were detected based on the map of the variographic surface ( Figure A5A).
After several empirical attempts, the tolerance, lag, and maximum bandwidth were found to equal to 2100, 4100, and 2000 m, respectively. Therefore, through the omnidirectional variogram, we confirmed the nonexistence of a trend in the index data set H j(Q) , ratified by the preliminary spatial drift analysis. The omnidirectional variogram reached the sill for a distance of about 8000 m.
It was also observed that the directional variograms for 0º, 45º, and 90 • , with an angular tolerance of 45 • , presented similarities, entering a state of randomness after an 8000 m distance ( Figure A5B). Although anisotropy was observed on the variographic surface when the directional variogram was observed at 135 • , it was considered that the omnidirectional variogram behaved as an average in relation to the directional ones. Thus, without using any "normalization", by adopting a small search radius, the distribution o e Hj(Q) index was considered as isotropic, adjusting the omnidirectional variogram via the spheric PH), exponential (EXP), and Gaussian (GAU) models ( Figure B6). Thus, without using any "normalization", by adopting a small search radius, the distribution of the H j(Q) index was considered as isotropic, adjusting the omnidirectional variogram via the spherical (SPH), exponential (EXP), and Gaussian (GAU) models ( Figure A6).
The best theoretical adjustment was associated with the EXP model. Adjusted via the aforementioned EXP model, the variability of the H j(Q) index was estimated via normal kriging. In addition to being used to choose the EXP model, among others, the cross-validation process of the variographic model revealed a statistical indicative RMSE which was equal to 0.08715.