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

Abstract

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.

1. 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.

However, the commonly used methods to evaluate intrinsic vulnerability, to quantify the hazard degree of a contamination source, to classify the quality, and to estimate the use of groundwater reserves are subjective (shallow or deeper aquifers). Additionally, the lack of available geological, hydro-geochemical, and hydro-geological data is a source of uncertainty that can lead to errors in risks assessments. Therefore, proposing a risk model that can handle these limitations, that is, the subjectivities and uncertainties, is valuable for risk assessment.

1.1. Theoretical Background and Limitations

1.1.1. To Assessment Intrinsic Vulnerability

The typical way of assessing vulnerability involves methods structured by a Y-index, which results in the association of input data (X) (x₁, x₂, ...x_i), that is, of hydrogeological parameters. In general, this format imposes subjectivity on the methods, denoting them as conceptual (usual methods seek a simplified interpretation to explain complex interactions that occur in a simultaneous and systemic manner, such as vertical water flow phenomena and contaminant transport, as well as biophysical, chemical, and physical processes) and ambiguous (they define the universal set of X-parameters in imprecise terms, which are used to classify the Y-index into specific degrees of vulnerability).

Although recent work has been developed to overcome the assessment vulnerability limitations [3], this subjective, conceptual, and ambiguous nature has been circumvented using “intelligent models”, for example, the application of fuzzy models. The works of [4,5,6,7] applied a multi-criteria fuzzy model to assess vulnerability based on the X-parameters of the DRASTIC method [8]. Other works [9,10,11,12,13,14] proposed a fuzzy inference system (FIS) that permits systematic and simultaneous interactions of X-parameters to evaluate the Y-index.

Constructing an FIS requires the following three steps: “fuzzification”, inference, and “defuzzification”. The “fuzzification” is defined by attributing linguistic variables to the X-parameters and to the Y-index and by selecting a type of membership function (mf). Although this function can assume many different forms, the usual forms are triangular (mftri), trapezoidal (mftrap), and gaussian functions (mfgau). In the inference, all the applications of one type of FIS to assess the intrinsic vulnerability follow the Mamdani model [15], as they do in the “defuzzification” using the centroid method [16]. The FIS was not validated, nor were the performances measured in studies that apply a fuzzy model to assess vulnerability. In these studies, the FIS was only compared to the base-method.

1.1.2. Assessment of the Hazard Degree from Contamination Sources

It is necessary to identify the polluting agents and quantify the contaminant loads. This identification refers to the registration of sources [17,18]. The pollutant origin and its surcharge hydraulically (POSH) method, which classifies sources into levels of contaminant load generation, was proposed [19]. This approach considers the type of activity, contaminant-generating capacity, and associated hydraulic head. Recently, some studies proposed hazard indices (H), a function of the physico-chemical properties of contaminants [1,2,20,21,22]. In general, the studies reviewed that classify the hazard degree are based on the importance of the relationship among physico-chemical properties. Thus, they represent a pattern framing issue with an approach by the multi-criteria method, the objective of which is to classify the hazard degree. However, limitations include the operational difficulty of inserting n-properties and weighting these properties to quantify the hazard degree.

1.1.3. Assessment of the Natural Groundwater Quality

The Piper diagram is one of the methods proposed to assess the groundwater quality. The Piper diagram is used for the initial hydro-chemical analysis. Diagrams such as the U.S. Salinity Laboratory diagram and the Schoeller–Berkalof diagram classify water quality for use in agriculture or for human consumption, respectively. However, these diagrams do not synthesize information about the water quality and do not find organics or special contaminations as indices.

Thus, the works of [23,24,25] developed methods capable of classifying the groundwater quality with a Q-index using subjective, conceptual, and ambiguous formulations (the same is justified about limitations of the Y-index). An FIS may be used to overcome these limitations. Indices are used to classify water quality standards, the parameter weighting of which can be associated with the final intended use (e.g., potable or non-potable water). Thus, adopting a multi-criteria model that manages, in addition to the mentioned limitations, a judgment on the importance between parameters for multiple water use purposes is an asset during this process.

1.1.4. Estimation of the Groundwater Reserves

Although recent work has developed a stochastic assessment of nonpoint source contamination [26], a similar method was not applied to estimate the probability of use of groundwater reserves. In turn, the various forms of deterministic calculations of the groundwater reserves (G), the sum of permanent reserves (G_P), and regulating reserves (G_R) [27,28] consist of variables with large variability. The authors of [29] adopted average values for the saturated thickness of aquifers (b), effective porosity (η_e), area of occurrence (A), and water level variation (Δh) (during the dry and rainy season), which denote uncertainties that can be overcome using a probabilistic approach.

The probability that the G_P reserve will be used can be estimated with the safe use of the G_R reserve (and vice-versa). The objective of this study is thus to propose an fuzzy-probabilistic modelling approach for groundwater contamination risk assessment (FPM-risk) for groundwater contamination assessment, and its methodology is described and prescribed below. Its application was conditioned by strategies for the acquisition, identification, and cataloging of geological and hydrogeological data, source characteristics, water physicochemical parameters, and the physical variables of the shallow aquifer.

2. Proposed Modeling

2.1. FIS for Assessing the Intrinsic Vulnerability

The vulnerability assessment relied on the proposed implementation of a type of FIS, which uses the Fuzzy Logic Toolbox from Matlab® (version R2020a, development by MathWorks, Inc., Natick, MA, USA). For this purpose, the GOD-based method was used to obtain Y_GOD [30]—an index that refers to the acronym of the GOD method; that is, groundwater hydraulic confinement, overlaying strata, and depth to groundwater table—Y_GOD = (x_G) × (x_O) × (x_D), which classifies intrinsic vulnerability as insignificant (I), low (L), moderate (M), high (H), or extreme (E). The GOD method was utilized because it does not demand a significant parameter number and could used in a sedimentary aquifer.

The proposal to use an FIS to build the GOD-fuzzy model considered the following two principles: the hydraulic inaccessibility and contaminant attenuation capacity of the strata overlaying the vadose zone of the aquifer, which are associated with the aquifer lithology. Therefore, the parameters x_G, x_O, and x_D and Y_GOD were fuzzily. Thus, through the interaction among the degree of confinement (x_G) and the distance of the water level from the surface (x_D), the hydraulic accessibility index (HAI) was obtained and qualified in degrees as L, M, or H (Figure 1).

The fuzzification of vulnerability, x_G, x_O, x_D, HAI, and Y_GOD-Fuzzy* is illustrated in Figure 1 and Figure 2. The inference is processed by the Mamdani model, which is driven by the interaction among the x_i parameters and activated by rules. The number of rules is equal to the number of linguistic variables raised to the number of X-parameters [12]. Thus, there are nine rules that express the premise among x_G and x_D and the consequence of the HAI (e.g., Rule: “IF” the degree of confined is “Confined”, and the water level from the surface is “L”, “THAN” the hydraulic accessibility is “M”).

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:

$$RE=\frac{1}{\mathrm{N}}{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}\frac{\left|Y_{\mathrm{GOD}}-Y_{\mathrm{GOD}-\mathrm{Fuzzy}}{}^{*}\right|}{\mathrm{m}\mathrm{á}\mathrm{x}\left(Y_{\mathrm{GOD}},Y_{\mathrm{GOD}-\mathrm{Fuzzy}}{}^{*}\right)}\times 100,$$

(1)

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.

2.2. Fuzzy-Analytic Hierarchy Process (AHP) Model to Assess Hazard and Classify the Groundwater Quality

2.2.1. 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₁), mobility (p₂), degradability (p₃), solubility (p₄), and contaminant volume (p₅).

• p₁: 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₂: 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₃: 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₄: 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₄ 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₅: 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³/fuel station (FS).

Thus, to define the p₁, p₂, p₃, and p₄ 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. Contaminant properties (p_n).

Contaminants	Properties—pn
	p1 (mg/L)	p2	p3	p4 (mg/L)
Benzo(a)pyrene (BP)	10−5	5.95	1.84	10−3–10−1
Benzene (B)	10−2	1.75	2.44	104
Toluene (T)	0.70	2.07	2.94	104
Xylene (X)	0.50	2.25	2.91	104
Ethylbenzene (E)	0.30	2.23	2.81	104
Naphthalene (NT)	0.30	5.95	1.84	104

. Using Table 1, an apparent hazard pattern matrix was developed (

Table 2. Hazard pattern matrix (Z_p). VL—very low; L—low; M—medium; H—high; VH—very high.

Properties—pn	Classes of Apparent Hazard Standards—h
	1 (VL)	2 (L)	3 (M)	4 (H)	5 (VH)
p1 (mg/L)	103	50	5	0.5	10−5
p2	5.95	4.50	3.50	2.50	1.50
p3	4.75	3.75	2.75	1.75	0.50
p4 (mg/L)	10−3	10−1	10	103	104
p5 (m³)	15	30	45	100	150

).

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).

$$p_{n_{\left(m\right)}}={\displaystyle \sum }_{i=1}^m{w}_c{p}_n.$$

(2)

Using this procedure, it is possible, for example, to define that compounds such as BP and B are most important for the value p₁, because they are more toxic than the other compounds.

2.2.2. 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₁), chloride (q₂), total residuals (q₃), hardness (q₄), nitrate (q₅), and fluorine (q₆) 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₁: 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₂: 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₃: to classify the natural groundwater quality, a weight is conferred that is also referenced in the proposal in [35].
• q₄: when ingesting hard and very hard waters, renal and cardiovascular diseases can develop [34]. CaCO₃ 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₅: although NO₃^- 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₂⁻), causing the emergence of illnesses. In addition to methemoglobinemia, the authors of [34] noted that NO₂⁻ in the human body can react with amine substances, resulting in nitrosamine (carcinogenic agent). The maximum permitted concentration of nitrogen-nitrate (N-NO₃) in groundwater for human consumption was set by the WHO at ≤50 mg of N-NO₃/L.
• q₆: 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₅, q₂, q₆, q₁, q₄, and q₃ 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. Natural groundwater quality pattern matrix (Z_q). Un—unfit; Ac—acceptable; Go—good; Op—optimal.

Parameters—qn	Natural Groundwater Quality Degree—Qn
	1 “Un”	2 “Ac”	3 “Go”	4 “Op”
q1	<6–>9	6–9	7–8	7.5
q2 (mg/L)	>350	250–350	200–500	150
q3 (mg/L)	>1000	100–1000	200–500	350
q4 (mg/L)	>3500	350	300	200
q5 (mg/L)	50	45	20	5
q6 (mg/L)	>1.5	1.5	1.2	0.6–1

).

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”.

2.2.3. 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:

$${{\displaystyle \sum }}_{h=1}^h{u}_{h,j}=1 .$$

(3)

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:

$$d_{h,j}={\left\{{\displaystyle \sum }_{\left(p,q\right)=1}^{\left(p,q\right)}{\left[w_{\left(p,q\right)}(x_{\left(p,q\right),j}{}^{\text{'}}-z_{\left(p,q\right),h}{}^{\text{'}})\right]}^2\right\}}^{1/2} ,$$

(4)

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:

$$Min\left\{F\left(u_{h,j}\right)={\displaystyle \sum }_{h=1}^h{D}_{h,j}^2\right\}.$$

(5)

To obtain a Lagrangian function and optimize the classification process by restriction,

$$L\left(u_{h,j},{\lambda }_j\right)={\displaystyle \sum }_{h=1}^h{D}_{h,j}^2-{\lambda }_j\left({\displaystyle \sum }_{h=1}^5{u}_{h,j}-1\right),$$

(6)

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, (i) ∂L(u_h,j,λ_j)/ ∂L(u_h,j) = 0 and (ii) L(u_h,j,λ_j)/ ∂L(λ_j) = 0. Using the partial derivatives of the Lagrange function, the pertinence function of u_h,j is obtained for class h by using the following:

$$u_{h,j}={\left(d_{h,j}^2{\displaystyle \sum }_{h=1}^h\frac{1}{d_{h,j}^2}\right)}^{-1}.$$

(7)

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).

$$H_{j\left(P,Q\right)}=\left(1,2,3,4,5\right)u_j{}^{*}={\displaystyle \sum }_{h=1}^5{u}_{h,j}^{*}h ,$$

(8)

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.

2.2.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:

f(G_R) = A × f(Δh) × f(η_e).

(9)

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.

3. 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³, 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³/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₃, 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 demand for groundwater was characterized in a study conducted by [38], which indicated a water quantity of slightly greater than 430,000 m³/day, of which approximately 20% comes from aquifers (around 86 thousand m³/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³ (54.75 million m³/year).

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³/day (3.65 million m³/year). For the island of Outeiro, with a population of close to 90,000 inhabitants, the daily demand is close to 18,000 m³ (6.57 million m³/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.

4. Results and Discussions

4.1. 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₀” 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.

4.2. 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).

4.3. 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₀” 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).

4.4. 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², 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:

• 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³ (54.75 million m³/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³/day (4.5 million m³/year). On the island of Outeiro, which also depends on groundwater, the population of close to 90,000 inhabitants demands approximately 22,500 m³ per day (8.2 million m³/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³, 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).

4.5. 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).

5. 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.