Abstract

Increased anthropogenic activity in urban areas has exacerbated the vulnerability of groundwater resources. The AVI, GOD, SINTACS, and DRASTIC methods were used to analyze groundwater vulnerability in Pangkalpinang City. Schlumberger vertical electrical sounding was used to determine the lithology and aquifer configuration in the study area. There are three vulnerability index areas in the city of Pangkalpinang. Low levels of aquifer vulnerability were generally found in the southeastern and northwestern parts of the study area, whereas high levels of aquifer vulnerability were discovered in the northern and southern parts of the study area. Areas with low aquifer vulnerability levels generally have low hydraulic conductivity values on the protective layer. In these areas, groundwater extraction is possible with a reasonable extraction pattern. Industrial areas can also be built by considering environmental aspects. In an area with high-level aquifer vulnerability, groundwater pollution must be considerably managed. The areas should not be designated for industrial areas and excess groundwater extraction.

1. Introduction

1.1. Background

Pangkalpinang is the capital of Bangka Belitung Islands Province (Figure 1). It has a population of 212,727 people, with a growth rate of 0.85% and an area of 118.41 km2 (Pangkalpinang Central Bureau of Statistics, 2020). It is the center of the economy in Bangka Belitung Islands Province with rising economic growth rates from year to year. As the economy of a city develops, the demand for water resources will continue to grow. The amount of water needs in the city of Pangkalpinang reached at least 2 million m3 annually (Pangkalpinang Central Bureau of Statistics, 2020). With the decline in surface water quality, massive exploration of groundwater is unavoidable [1].

Increased anthropogenic activity in urban areas will increase vulnerability of existing groundwater resources [2]. Anthropogenic wastes will threaten groundwater quality and quantity [3]. Proper groundwater identification and management are needed to support the sustainability of groundwater resources. Indexing areas that have potential groundwater resource vulnerability will facilitate better groundwater management [4].

Aquifers are highly threatened by urban waste [5]. Contaminated groundwater will cause severe health problems for humans, such as cancer [6]. This chain of problems certainly will greatly impact the economy of a city if not taken seriously [7].

Many researchers used different approaches to identifying groundwater vulnerability [8]. Some researchers used the aquifer vulnerability index (AVI) method to identify and analyze groundwater vulnerability [9]. GOD, SINTACS, and DRASTIC methods were commonly used by other researchers [10].

The aquifer vulnerability index (AVI) is an approach introduced by Van Stempvoort et al. [11]. The aquifer vulnerability index (AVI) method attempts to calculate the aquifer susceptibility index based on the thickness value and the hydraulic conductivity of the aquifer cover layer. The thicker and lower the hydraulic conductivity value of the aquifer cover layer, the lower the aquifer susceptibility index of the region [11]. In hydrogeological systems, this layer is commonly referred to as aquitard, aquiclude, or aquifuge layer [12]. Obiora and Ibuot [2] and Van Stempvoort et al. [11] refer to this layer as the protective layer of the aquifer. This research used the geoelectric method to get the value of aquifer thickness. This threat will greatly increase if the system does not have a good protective layer of the aquifer [12]. In urban areas, vertical percolation is a threat to groundwater [2].

The GOD method utilizes three parameters to determine groundwater vulnerability, i.e., groundwater confinement, overlying strata, and depth of groundwater. This method is introduced by Foster [13]. The SINTACS method used seven parameters to identify groundwater vulnerability. The SINTACS method used in this study was developed by Civita and De Maio [14] to evaluate relative groundwater pollution vulnerability. The seven parameters were water table depth (S), effective infiltration (I), unsaturated zone (N), soil media (T), aquifer media (A), hydraulic conductivity zone (C), and topographic slope (soggiacenza (depth of water), infiltrazione efficace (effective infiltration), non saturo (vadose zone), tipologia della copertura (soil cover), acquifero (aquifer), conducibilità idraulica (hydraulic conductivity), and superficie topografica (slope of topographic surface)). DRASTIC methods were introduced by Aller et al. [15]. DRASTIC methods also used seven parameters to determine groundwater vulnerability. Depth to groundwater, net recharge, aquifer media, soil media, topography, impact of vadose zone, and hydraulic conductivity were the parameters that are needed to analyze groundwater vulnerability.

According to Maria [16], the GOD method is suitable for designing large areas such as land management while DRASTIC has good accuracy and more real use in geoenvironmental detailed studies. The SINTACS method generates very high vulnerability zones in the areas concerned with surface waters and aquifer interactions [16]. AVI methods are suitable for regional groundwater basin analysis [11].

The geoelectric method is one of the geophysical methods for mapping subsurface conditions based on the electrical parameters [17]. The resistivity method maps subsurface conditions based on rock resistivity parameters [18]. The resistivity method has been widely used in hydrogeological research [19]. It is also widely used in research on aquifer susceptibility [20]. Vertical electrical sounding is a data collection technique widely used in geoelectrical research [21]. Schlumberger arrays are electrode configuration techniques in vertical electrical sounding that are widely used in hydrogeological investigations [22]. Schlumberger arrays have a good vertical resolution in providing a clear view of subsurface conditions [23].

This study’s objective was to analyze groundwater vulnerability in Pangkalpinang City and proposed a vulnerability index by integrating commonly used methods (AVI, GOD, SINTACS, and DRASTIC). Possible groundwater management patterns are based on the new proposed vulnerability index considerably.

2. Geology and Hydrogeology

2.1. Geology

According to Mangga and Djamal (1994), there are three geological formations in the study area, i.e., alluvium, Tanjunggenting Formation, and Pemali Complex (Figure 2). The alluvium is a Holocene surficial deposit. The alluvium in the study area consists of boulders, cobbles, pebbles, sand, clay, and peat. The Tanjunggenting Formation is a Triassic sedimentary rock, consisting of an alternation of metasandstone, sandstone, clayey sandstone, and claystone with lenses of limestone. The formation is well-bedded, strongly folded, jointed, and faulted. The Pemali Complex is an upper Paleocene to lower Permian formation consisting of metamorphic rocks. It consists of phyllite and schist with intercalation of quartzite and limestone lenses; the formation is jointed, folded, and faulted.

2.2. Hydrogeology

According to Sukrisna and Sudadi (2002), the occurrence of groundwater and productivity of aquifers in the study area consist of aquifers in which flow is intergranular and aquifers (fissured or porous) of poor productivity and regions without exploitable groundwater (Figure 3). For an aquifer in which flow is intergranular, it is a locally moderately productive aquifer. It is mostly an incoherent aquifer of low thickness and transmissivity. The groundwater table is generally less than 3 m below the ground surface; the well field is less than five l/sec—alluvium and swamp deposits are composed of pebble, sand, clay, and peat with generally moderate to high permeability. For aquifers (fissured or porous) of poor productivity and regions without exploitable groundwater, those are almost poorly productive aquifers of local importance. With generally low transmissivity, locally in favorable sites, a small discharge of springs can be expected. Limited shallow groundwater can be obtained in the valleys and weathered or fractured zone of solid rocks. There are alternations of metasandstone, sandstone, clayey sandstone, claystone, slate, mudstone, shale, and cherts with lenses of limestone, with low permeability and locally moderate permeability in the weathered zone. There are regions without exploitable groundwater in the study area. It mostly consists of phyllite and schist with intercalation of quartzite and limestone lenses and is generally low permeabile to impermeable.

According to Sukrisna and Sudadi (2002), there are no seawater intrusion phenomena in the study area. Groundwater flow is towards the eastern and northeastern parts of the study area (Sukrisna and Sudadi, 2002).

The annual precipitation in the study area was 2099.6 mm (Pangkalpinang Central Bureau of Statistics, 2020). Figure 4 shows the daily rainfall in Pangkalpinang City. The average air temperature of Pangkalpinang in 2019 was 27°C with the highest one at 32.40°C and the lowest at 23.50°C (Pangkalpinang Central Bureau of Statistics, 2020). Figure 5 shows the daily temperature in the study area.

Table 1 shows groundwater quality data in the study area in 2018.

3. Methods

3.1. Vertical Electrical Sounding

Electrical Resistivity Tomography (ERT) and vertical electrical sounding (VES) are measurement techniques widely used in groundwater exploration [23]. Vertical electrical sounding (VES) is effective in characterizing aquifer conditions [25]. The Schlumberger array electrode is an effective configuration in mapping subsurface conditions because it has an excellent vertical resolution [23]. The output value of the field measurement is apparent resistivity calculated by the equation where :

AB and MN are current electrodes, and potential electrode (meter) and are voltage and current at a particular measurement, respectively. The field data calculations resulted in apparent resistivity values, meaning the values are still affected by the resistivity values above the layer. Materials used for the VES survey included a set of GL-4200 geoelectric devices, Garmin 64 GPS, laptops, and IPI2Win software. IPI2Win software was developed by Moscow University. This software was used to get the value of true resistivity through the inversion process. The output of the IPI2Win software is the value of the resistivity, thickness, and depth of a layer. Apart from that, the output from IPI2Win can also be in the form of inversion curves.

In this study, eight measurement points of vertical electrical sounding were scattered around the city of Pangkalpinang (Figure 2). This study used the Schlumberger array electrode for electrode configurations. The maximum stretch of meters.

Qualitative interpretation was made to get subsurface lithology (drawn from the research of Mangga and Djamal (1994)). Resistivity layer values were also interpreted according to the geological information provided in Mangga and Djamal (1994). The results of this interpretation were subsequently used in the calculation of the aquifer vulnerability index.

3.2. Aquifer Vulnerability Index

The aquifer vulnerability index method is introduced by Van Stempvoort et al. [11] based on the thickness and hydraulic conductivity values of an aquifer protective layer. Lithology layer thickness values were obtained based on the interpretation of the results of vertical electrical sounding inversion, while hydraulic conductivity values were  m/day [26],  m/day [27], and 1 m/day [28] for schist, sandstone, and silty sand, respectively [2, 9]. Once these two values are determined, the value of hydraulic resistance can be obtained using the following equation [2, 11]: where is the thickness of the protective aquifer layer and is the hydraulic conductivity of the protective layer.

Table 2 shows the relationship of the aquifer vulnerability index to hydraulic resistance. A contingency was used to visualize the aquifer vulnerability index in Pangkalpinang City, with data gridding using inverse distance weighted [2931]. Contingency uses ArcGIS Pro 2.6 software to process and visualize the data (ESRI 2016).

3.3. GOD Method

According to [8], the DRASTIC and GOD models were chosen as the most appropriate for groundwater vulnerability study because (1)they are simple and take into account several parameters of the natural environment which control the contamination process of the aquifers(2)they include a high number of layers of input data that limit the impacts of errors of the individual parameters on the final result(3)these parameters are in general available for their evaluation

Many authors have evaluated the vulnerability rate of aquifers using these two methods and established good correlations. According to Polemio et al. [32], the GOD method is useful for mapping large areas with high vulnerability contrasts, whereas DRASTIC is useful for any aquifer [33].

The GOD method is a vulnerability assessment method developed in the United Kingdom by Foster [13]. The GOD index can be divided into six categories (Table 3), from 0 to 1 [13]. The GOD parameter indexes (groundwater confinement, overlying strata, and depth to groundwater) were created as described for the DRASTIC model. The index of vulnerability is given according to the following formula [13]:

3.3.1. Groundwater Confinement

Groundwater confinement is one of the parameters used in the analysis of groundwater vulnerability [34]. Groundwater confinement can be interpreted based on subsurface conditions using vertical electrical sounding (VES). The groundwater confinement type in the study area consists of unconfined and confined aquifers. Estimating the groundwater confinement value can be obtained from GOD classification, as seen in Methods. The unconfined aquifer in the GOD method category has a value of 1, while the confined aquifer has a value of 0.2 (Table 4).

Estimating values for unconfined aquifers based on the characteristics of unconfined aquifers tends very close to the ground level and has no aquitard. The confined aquifer assessment is based on the features of the confined aquifer, where the aquifer lies between the confined layer of the aquifer.

3.3.2. Overlying Strata

It is a lithological character and degree of consolidation of the vadose zone or confining beds—estimation of vadose zones and confining beds based on vertical electrical sounding (VES). There are three types of vadose zone and confining bed in the study area. The three types of overlying strata are igneous/metamorphic, sandstone, and alluvial. Igneous/metamorphic and alluvial have a value of 0.6, while alluvial and sandstone have a value of 0.7 (Table 5).

3.3.3. Depth of Groundwater

Groundwater level depth is another important parameter in determining groundwater vulnerability [34]. Table 6 shows the depth of groundwater value/rating in the study area. The groundwater level condition is known through the results of monitoring wells of 10 wells scattered around Pangkalpinang. Groundwater depth data is also used in the analysis of groundwater vulnerability using DRASTIC and SINTACS methods. Groundwater level depth data in 10 monitoring wells were then interpolated using the IDW interpolation technique in ArcGIS Pro 2.6 software. The IDW interpolated gridding results were then extracted to 10 VES points in the study area to obtain the groundwater level depth at each VES point. The shallower the groundwater level, the higher the groundwater vulnerability, while the deeper the groundwater level in a place, the less vulnerability to contamination of groundwater. The determination of the depth of groundwater value refers to the GOD classification. Figure 6(c) shows the spread of depth of groundwater in the study area.

3.4. SINTACS Method

The SINTACS method is a PCMS developed by Civita and De Maio [14] to assess groundwater’s intrinsic vulnerability, and it is partially derived from the worldwide used DRASTIC [15]. The acronym SINTACS originates from the Italian words of the seven parameters employed in the method: (i)Soggiacenza (depth of water)(ii)Infiltrazione efficace (effective infiltration)(iii)Non-saturo (vadose zone)(iv)Tipologia della copertura (soil cover)(v)Acquifero (aquifer)(vi)Conducibilità idraulica (hydraulic conductivity)(vii)Superficie topografica (slope of a topographic surface)

The parameters are multiplied by a weight related to a given hydrogeological setting. The form of the equation is the following [35, 36]: where is the rating of each parameter and is the corresponding weight suggested from [14].

3.5. DRASTIC Method

US EPA developed the DRASTIC model and extensively used it to assess pollution potentials in the United States and other countries around the world [15]. This model uses seven physical and hydrodynamic parameters of the natural environment to assess the aquifer’s vulnerability, such as depth to the water table, recharge, aquifer media, soil media, topography, vadose zone, and hydraulic conductivity. This method’s index is often used to standardize the evaluation of groundwater pollution potential within various hydrogeological settings. Its calculation assumes that (1) the contaminant is introduced at the ground surface, (2) the contaminant is flushed into the groundwater by precipitation, (3) the contaminant has the mobility of water, and (4) the area evaluated is 0.4 km2 or larger [15]. Each of the seven parameters is classified in a range of values or according to the dominating medium. It reflects its impact on aquifer vulnerability. Each parameter has been assigned a rating between 1 and 10, corresponding to the low and high potential of contamination. Each parameter’s rating is combined with its weight for all the study area sectors. The DRASTIC index is finally computed by implying linear combinations of the products of rating and weights:

The index corresponds to the rating and the index to the weight. The higher DRASTIC index corresponds to the high potential contamination and vice versa. This index is a relative value without dimension, allowing a comparison between the sites. The result can vary between 23 and 230. These extreme values are rarely reached. The computed values generally vary between 50 and 200 [8].

4. Result

4.1. Aquifer Vulnerability Index

In general, the resistivity values in the study area were in the range of 96-1988 Ωm (Table 7). The silty sand layer, which is part of the alluvium formation, was in the range 213-1043 Ωm with an average value of 628 m. The sandstone layer, which is part of the Tanjunggenting Formation, was in the range of 96-1254 Ωm with an average value of 675 Ωm. The schist layer, which is part of the Pemali Complex formation, was in the range of 682-1988 Ωm with an average value of 1335 Ωm. Schists had the highest resistivity value of all the other rock layers (sandstone and silty sand). According to Telford et al. [37] and Reynolds [23], metamorphic rocks have relatively higher resistivity values compared to sedimentary rock and loose sediment. The high resistivity value of schist rocks is influenced by age and degree of compaction and lithification of these rocks [38]. The depth and thickness of each rock layer vary greatly (Figure 7).

The aquifer layer had a range of resistivity values that vary from 682 to 1347 Ωm, where the aquifer was found in the fissure and fracture of schist and sandstone rocks as revealed by Sukrisna (2004). High aquifer values are common in aged rocks such as schist and sandstone [21, 39]. Apart from being a limited aquifer, schist can also act as an aquifuge layer, that is, a layer that cannot hold or drain water [28, 40].

The values of resistivity, depth, and thickness of the protective aquifer layer at the study area were varied (Figure 8). This layer is located above the aquifer layer. This layer is also known as the vadose zone [2, 41]. Generally, in the study area, this layer was lithologically filled with silty sand and upper sandstone layers, except in VES01 where the protective layer consisted of a schist layer. The inversion results show that the resistivity of this layer tended to be smaller than the layer below it based on the interpretation of each geoelectric measurement point (Table 7).

The value of aquifer resistivity tended to be lower in the southeast part of the study area. This is likely due to intensive secondary porosity in the form of rock fractures resulting from tectonic activity on the island of Bangka (Mangga and Djamal, 1994). Secondary porosity in fractures can reduce rock resistivity [37].

The values of the sum of the hydraulic conductivity protective layer varied in the study area, ranging from 0.00021 to 1.00911 m/day. Hydraulic conductivity values tended to be low in the southeastern part of the study area, but high in the southern part of the study area. Different lithologies cause this discrepancy of conductivity values at each VES measurement point. High hydraulic conductivity is closely related to rock properties that tend to be more permeable. The southeastern part, especially the measurement point of VES01, was composed of schist. Schist has a smaller hydraulic conductivity index due to this rock’s ability to store and transmit groundwater than sedimentary rock (sandstone) and loose sediment (silty sand) [28]. In the southwestern and northern parts of the study area, high hydraulic conductivity values in the protective layer can be attributed to the loose sediment in the form of silty sand. The silty sand layer is relatively easier to store and loose water compared to schist and sandstone [42].

The sum of the thickness layer values ranges from 63.695 to 95 m. The southern part of the study area tended to have a thick protective layer, while the northern part tended to have a thinner protective aquifer. This difference might be caused by the pattern of depositional facies existing in the study area (Mangga and Djamal, 1994). VES03 had the thickest protective layer (95 m), while VES04 had the thinnest protective layer (63.695-95 m).

There are four regional categories based on the aquifer vulnerability index in Pangkalpinang City (Figure 9). The regions are divided into areas with high, moderate, low, and extremely low levels of aquifer vulnerability (Table 8). The categorization was based on the aquifer vulnerability index proposed by Stemvoort et al. [11]. There were no areas with very high levels of aquifer vulnerability in the study area. Very low aquifer vulnerability levels were generally found in the southeastern and northwestern parts of the study area.

In contrast, high aquifer vulnerability levels were found in the northern and southern parts of the study area. Areas with very low aquifer vulnerability levels generally have low hydraulic conductivity values on the protective layer. These areas are protected from vertical pollutant percolation [2].

VES01 had a small vulnerability index caused by the low hydraulic conductivity of the protective layer. The protective layer on VES01 was in the form of schists with the smallest hydraulic conductivity value among the rocks at the research location, namely, sandstone and silty sand. Even though the thickness of the protective layer in VES01 was only around 73 meters, the small hydraulic conductivity value has impeded the vertical water movement from entering the aquifer [43].

VES04, VES07, and VES08 had high vulnerability because the protective layers tended to have a high hydraulic conductivity value. These three points consisted of three types of rock lithology, namely, silty sand, sandstone, and schist. The silty sand and the sandstone layers act as the protective layers. These two layers are relatively more permeable than schists, thus allowing water to move vertically towards the aquifer layer [28].

VES02, VES03, and VES06 had a low vulnerability index because the protective layers at this point consisted of sandstone lithology. Although sandstone can hold and drain water, this layer tends to have a low hydraulic conductivity compared to the silty sand layer. The sandstone layer tends to be more massive and solid (compacted and verified) [40, 44]. VES05 had a moderate vulnerability index due to the thickness of the protective layer, although it had the same sum of hydraulic conductivity as VES04.

4.2. GOD Method

VES01, VES04, VES05, VES07, and VES08 are unconfined aquifers (Table 9). VES04, VES05, and VES08 are alluvial-type aquifers (Table 9) with groundwater flowing through pore media. VES01 and VES07 are metamorphic aquifers with groundwater flowing through fracture and fissure. Figure 6(a) shows groundwater confinement in the study area.

There are two-three types of overlying strata in the study area, i.e., weathered metamorphic, weathered sandstone, and alluvial. VES01 and VES07 were types of weathered metamorphic overlying strata (Table 9). VES02, VES03, and VES06 were weathered sandstone overlying strata (Table 9). VES04, VES05, and VES08 were overlying alluvial strata (Table 9). Weathered metamorphic has a value of 0.6; alluvial and sandstone have a value/rating of 0.7. Figure 6(b) shows the overlying strata in the study area.

VES07 and VES08 have groundwater levels of <5 meters, while the other six VES points have groundwater depths between 5 and 20 meters (Table 9). Rating/value for groundwater level meters is 0.9 (Table 6), while rating/value for groundwater level depth 5–20 meters is 0.8 (Table 6).

The groundwater vulnerability index in the study area has a GOD index range of 0.1–0.6 (Table 9). The groundwater vulnerability level at the study area can be divided into three categories: low vulnerability level, moderate vulnerability level, and high vulnerability level. The low vulnerability level has index values of 0.1-0.3, visualized with the green area in Figure 6. Low vulnerability values are in the southeastern and northwestern study areas. Moderate vulnerability values are located almost throughout the study area, generally in the middle and west of the study area. The moderate vulnerability area has a vulnerability index value of 0.3-0.5. The moderate vulnerability area is visualized in yellow on the image. Areas with high vulnerability indexes are in the southwestern, northern, and northeastern parts of the study area. Areas with high vulnerability indexes have vulnerability value indexes of 0.5-0.7. Areas with a low vulnerability index have characteristic types of confined aquifers, vadose zones consisting of weathered rock, and groundwater level depths of 5–20 meters. Areas with a high vulnerability index have characteristics of an unconfined aquifer type, vadose/overlying strata consisting of loose sediment, and shallow groundwater level depth.

4.3. SINTACS Method
4.3.1. Soggiacenza (Depth of Water)

Groundwater level data were obtained from extrapolation of groundwater level monitoring wells from gridding results. Depth of water is defined as the piezometric level’s depth (both for confined or unconfined aquifers) regarding the ground surface. It was a great impact on the vulnerability because its absolute value, together with the unsaturated zone characteristics, determines the time of travel (TOT) of a hydrovectored or fluid contaminant and the duration of the attenuation process of the unsaturated thickness, in particular, the oxidation process due to atmospheric O2. Therefore, the SINTACS method rating of depth-to-groundwater decreases with an increase of the depth, i.e., an increase in the unsaturated zone’s thickness within the range 10-1.

The groundwater level value is in the range of 4.3–7.7 m. The rating value for groundwater level depth is 6.2–7.2 (Table 10). The weighting value for this parameter is 5. The total weighting value and rating in the study area are in the range of 30–36. VES07 and VES08 have the shallowest depths. In contrast, VES01 and VES02 have the deepest depth of the groundwater level.

4.3.2. Infiltrazione Efficace (Effective Infiltration)

The second SINTACS parameter is effective infiltration/groundwater recharge. The role that the effective infiltration plays in aquifer vulnerability assessment is very significant because of the dragging down of the pollutant’s surface and their dilution, first during the travel through the unsaturated zone and then within the saturated zone. Direct infiltration is the only or widely prevalent component of the net recharge in all the areas where there is no interflow linking aquifers or surficial water bodies or no irrigation practices using large water volumes. The groundwater recharge value for each VES measurement point is obtained from the gridding recharge value’s extrapolation based on calculations using ESPERE v2 software.

Recharge groundwater at the study area is in the range of 215–334.9 mm/year. To get the rating value, the recharge value that has been obtained is matched to the curve. The rating value for the effective infiltration/recharge parameter ranges from 8.4 to 9.2 (Table 10). The effective infiltration/recharge parameter has a weighting value of 5.

4.3.3. Non Saturo (Vadose Zones)

The third parameter in the SINTACS is non saturo (vadose zone). The unsaturated zone is the “second defense line” of the hydrogeologic system against fluids or hydrovectored contaminants. A four-dimensional process takes place inside the unsaturated thickness in which physical and chemical factors synergically work to promote the contaminant attenuation. The unsaturated zone attenuation capacity is assessed starting from the hydrolithologic features (texture, mineral composition, grain size, fracturing, karst development, etc.).

The study area’s vadose zone consists of weathered igneous/metamorphic, weathered sandstone, and alluvial (Table 10). The vadose zone with weathered metamorphic type is generally located in the southeastern part of the study area where in this area, there are complex Pemali formations. Vadose zones with weathered sandstones are typically scattered in parts of the study area where the parent rock is the sandstone of the Tanjunggenting Formation. Vadose zones with alluvial types typically fill coastal regions and near rivers.

The results of matching values with curves indicate that in the study area, the rating value is in the range of 4–7.5. This parameter has a weighting value of 5. The total result between the rating value and the weighting value indicates that the vadose zone has a range of 20–37.5.

4.3.4. Type of Coverage (Soil Cover)

The fourth parameter in SINTACS is soil cover. This is the “first defense line” of the hydrogeologic system: several important processes take place inside the soil that built up the attenuation capacity of a contaminant traveling inside a hydrogeologic system and, therefore, in aquifer vulnerability assessment and mapping. Soil is identified as an open, three-phase accumulator and transformer of matter and an energy subsystem which develops through the physical, chemical, and biological alterations of the bottom lithotypes and of the organic matter that it is made up of. Soil types in the study area based on FAO DSMW consist of Dystic Histosols (Od) and Orthic Podzols (Po). More specifically, there are two types of soil according to the criteria of Civita and De Maio [14], namely, silty clay loam and sandy loam. The soil types’ ratings are obtained by matching the soil type to the curves in the figure.

The study area’s soil type rating is in the range of 3.5–6.5 (Table 10). The weighting value for this parameter is 3. The total weighting value multiplied by the rating is in the range of 10.5–19.5.

4.3.5. Acquifero (Aquifer)

The aquifers in the study area are fracture, fissure, and porous aquifer. Based on the category of aquifers based on Civita and De Maio [14], aquifers in the study area consist of metamorphic fissure, sandstone, and coarse alluvial deposits. The rating value is obtained by using a curve match [14].

The rating is in the range of 3.5–8.5 (Table 10). The weighting value for this parameter is 3. The total result times between weighting values and ratings are in the range of 10.5–25.5.

4.3.6. Hydraulic Conductivity

Hydraulic conductivity represents the groundwater’s capacity to move inside the saturated media, thus the mobility potential of a hydrovectored contaminant, which in density and viscosity is almost the same as the groundwater. In the SINTACS assessment context, the hydraulic gradient and the flux cross-section are equal. This parameter determines the aquifer unit yield and flow velocity towards the tapping work effluences indicating the risk targets.

The value of hydraulic conductivity in the study area varies greatly. The range of hydraulic conductivity values in the study area ranges from 2.4-09 m/s to 1.1-05 (Table 10). The rating for each hydraulic conductivity is obtained by matching the curves [14]. The rating is in the range of 0.2–3.5. The weighting size for this parameter is 3. The total yield times between rating and weighting are in the range of 0.6–10.

4.3.7. Superficie Topografica (Slope of Topographic Surface)

The last parameter in SINTACS is topography. The slope for each VES point is in the range of 2–12%. The rating is in the range of 6-8 (Table 10). The weighting value for this parameter is 3. The total result of times between rating and dissing is in the range of 18–30.

4.3.8. SINTACS Vulnerability Index

Vulnerability index values are in the range of 90–189.8. Regions with a low vulnerability index is in the southeastern section (Figure 10(g)). The area is visualized in green. Regions with intermediate vulnerability indexes are located in almost all parts of the research area. This area is visualized in yellow. Areas with high vulnerability indexes are in the southwestern, north/northwestern, and northeastern parts of the study area. The area is visualized in red. Regions with high vulnerability indexes are associated with shallow mat depth, high groundwater recharge, water-flowing zone type, water-relative soil type, high permeability aquifer type, high hydraulic conductivity value, and ramped topographic index. A low vulnerability index is associated with deep mat depth, low groundwater recharge, water-flowing vadose zone, soil type that does not pass water, aquifer type with low permeability, low hydraulic conductivity value, and steep topographic index.

4.12. DRASTIC Method

`

4.12.1. Depth of Groundwater

The groundwater surface depth condition is known through the extrapolation of 10 wells from monitoring well measurement data to 8 VES points scattered around the city of Pangkalpinang. The groundwater level’s depth is in the range of 4.3-7.7 m below the local ground level. The measurement results are used to determine the parameter class, and each is multiplied by the weight. Based on the depth of the groundwater surface, this area is a confined aquifer and fractured aquifer type. The groundwater surface depth in the study area is the groundwater surface in groundwater dominance, where deep groundwater will be difficult to pass by contaminants.

4.12.2. Net Recharge

The data used to obtain groundwater recharge values is rainfall data derived from the Tropical Rainfall Measuring Mission (TRMM) data which can be obtained through the website https://giovanni.gsfc.nasa.gov/giovanni/, and potential evaporation data from Climate Forecast System Reanalysis (CFSR) can be obtained free of charge via the https://clim-engine-development.appspot.com/fewsNet; this data is then processed using the ESPERE v.2 software which can be obtained via https://www.brgm.eu/scientific-output/scientific-software/espere-estimating-effective-rainfall-aquifer-recharge.

Based on the distribution map analysis results and the value of the net recharge parameter, the net recharge at the study area is in the range 215-334.9 mm/year, including in the high net recharge category (Table 11). A high net recharge level affects groundwater pollution because incoming rainwater dilution causes the contaminants to dissolve easily and move towards groundwater.

4.12.3. Aquifer Media

Based on the aquifer lithology, groundwater availability at the study area can be grouped into three aquifer systems with different rock passage rates, namely, the aquifer system with flow through the space between grains, fractures and spaces between grains, and fractures or fissures (Sukrisna and Sudadi, 2002).

An aquifer system with flow through the space between grains is found in loose sediments, namely, alluvium and coastal deposits, composed of loose material from clay to gravel size with varying degrees of passing. This aquifer distribution occupies the coastal plain area, and there are localities in river valleys in hilly areas. This aquifer system is classified as locally productive and moderately productive.

An aquifer system with flow through the space between grains and fractures is composed of solid to less cohesive rocks. Rocks included in this aquifer system are a group of rock formations of the Tanjunggenting Formation consisting of malih sandstones, sandstones, clay sandstones, claystone, shale mudstone, shale, and chert with limestone lenses (Table 11). The distribution is quite evenly distributed in the study area. This formation is generally low graduated, locally graded moderately in the rock weathering zone.

Aquifer systems with flow through fractures/fissures are found in igneous and metamorphic rocks. The digger is complex, consisting of metamorphic rocks such as filites and schists with quartzite inserts and marble lenses. Transmisivty is generally low; local transmissivity is moderate in the weathering and fracture zones.

4.12.4. Soil Media

Soil texture parameters are related to the soil type. The basis for determining the soil type is obtained from the standards set by the FAO Digital Soil Map of the World (DSMW). The results of weighting the soil texture parameters are shown in the table. Soil types in the study area based on FAO DSMW consist of Dystic Histosols (Od) and Orthic Podzols (Po). Histosol is a soil composed mainly of organic materials. They are defined as having 40 centimeters (16 in) or more of organic soil material in the upper 80 centimeters (31 in). Mystic histosols were other histosols with a pH H2O (1 : 5) of less than 5.5, at least in some soil between 20 and 50 cm from the surface. Podzols are the typical soils of coniferous or boreal forests. Podzols can occur on almost any parent material but generally derive from either quartz-rich sands and sandstone or sedimentary debris from magmatic rocks, provided there is high precipitation. Most podzols are poor soils for agriculture due to the sandy portion, resulting in a low level of moisture and nutrients. There are two types of soil in the study area, namely, silty clay loam and sandy loam (Table 11).

4.12.5. Topography

Slope parameters were obtained from the SRTM Digital Elevation Model (DEM) and weighting results. The study area, which is dominated by a flat slope (Table 11), tends to collect water and increase infiltration, thereby accelerating the movement of contaminants. Meanwhile, the slope is steep-very steep because it increases the run-off so that groundwater is not easily contaminated. Thus, the principle of gravity will accelerate the movement of contaminants.

4.12.6. Impact of Vadose Zone

The unsaturated zone parameter is obtained from the soil texture map based on the soil’s grain size, assumed to be the surface layer, and the weighting result (Table 8). This parameter is very influential on pollution: a larger soil grain size and porous soil conditions will help move contaminants to the aquifer.

4.12.7. Hydraulic Conductivity

Soil hydraulic conductivity () is closely related to soil grain size distribution and porosity. This aquifer parameter greatly determines the sustainability of groundwater in an area (Hutasoit, 2009), while hydraulic conductivity values were  m/day [26],  m/day [27], and 1 m/day [28] for schist, sandstone, and silty sand, respectively [2, 9]. High values will accelerate contaminants, moving towards aquifers.

4.12.8. DRASTIC Vulnerability Index

Low vulnerability index areas can be found in the southeastern and northwestern parts of the study area (Figure 11(h)). The moderate vulnerability index area can be found in almost all parts of the study area. Moderate vulnerability index areas dominate the study area. High vulnerability index areas can be found in southwestern, northern, and northeastern parts of the study area.

5. Discussion

The classification of the vulnerability index is based on the GOD method [13], DRASTIC (Civita and De regibus, 1995, Corniello drr., 1997), SINTACS (Civita and De Maio, 1997; Al Kuisi et al., 2006), and AVI [11] presented in Table 12. There are three classifications for each method, namely, high, moderate, and low.

Figure 12 visualizes the vulnerability index at the study area based on the AVI, GOD, SINTACS, and DRASTIC methods, where each method requires specific parameters from one method to another.

5.1. Comparison Test

Based on Table 13, the study area is dominated by regions with a moderate vulnerability index. Areas with a high vulnerability index have a smaller area than those with a low vulnerability index, except for the AVI method. The area of the high vulnerability index has a larger area than areas with low vulnerability. The DRASTIC method has a total area of the highest moderate vulnerability (66.06558%). The AVI method has a high vulnerability index area, which is the highest among other methods (21.56044). The difference in input parameters is one of the main factors for the difference in the vulnerability index’s value.

To test how significant the difference is, a statistical test in the form of a -test is carried out. To test the correlation between methods, a correlation test was carried out. The correlation test results are presented in Table 14, while the results for the -test are shown in Table 15.

5.1.1. Correlation Test

The strongest correlation method was SINTACS and DRASTIC (0.997), while the method with the weakest correlation was AVI-GOD (0.742). The AVI method strongly correlates with the SINTACS method, while the AVI method has the lowest correlation with the GOD method. In general, the correlation index value between methods ranges from 0.742 to 0.997. This correlation value can be categorized as moderate-strong correlation.

5.1.2. -Test (Two-Tailed)

If the value is less than the significance level, the difference between means is statistically significant. Otherwise, the value is greater than the significance level; the difference between means is not statistically significant. In this case, we will use two-tailed, which is the value for the two-tailed form of the -test. Because all value is greater than the standard significance level of 0.05, we can conclude that the vulnerability index area is not different.

The -test result shows that the difference in the four methods’ vulnerability area’s value is not significant. The correlation results show that the area of the four methods has a moderate-strong correlation. This shows the similarity in results between the four methods. Figure 13 shows the comparison of the area of vulnerability based on the four methods. The histogram shows that only AVI has a different value for the three other methods for the low vulnerability area. The four methods confirm the area value for the moderate vulnerability area, which has similarities between the four methods. For the high vulnerability area, only AVI shows a higher area than the other three methods.

To obtain a reliable vulnerability area and because the four methods have similarities in area to one another, a vulnerability index map is created by combining the four methods overlay (Figure 14).

Areas with a moderate vulnerability index still dominate the research area, where this area is generally located in the middle and spreads to the northwestern part. Most of the areas with a low vulnerability index are located in the study location’s southeastern part. Low vulnerability areas are found in the northwestern and northeastern parts of the study area. Areas with a high vulnerability index are located in the southern and northern parts of the study area. The study area’s hydrological and hydrogeological conditions greatly influence groundwater vulnerability, where the hydrogeological parameters in the form of hydraulic conductivity, aquifer type, and groundwater level depth have a significant effect on the vulnerability value index. Areas dominate Pangkalpinang City with a medium vulnerability index (Table 16).

5.2. Groundwater Management and Protection

The northern and southern parts of the study area had a high index of aquifer vulnerability. The protective layer of the aquifer in these regions generally consisted of silty sand and sandstones. This situation has caused the hydraulic conductivity in this region to be high and resulted in high aquifer vulnerability index values. These regions need serious attention. The lithological conditions in these regions allow water to percolate through the aquifer layer more easily and quickly [43, 45]. Proper groundwater management is needed for these conditions [7]. Industrial activities in these two regions must be limited, and the areas should be used as green zones or urban forests [46]. In such regions, excessive groundwater extraction must be prohibited to maintain groundwater sustainability. Regulations are also needed [47]. In areas with a low aquifer vulnerability index (based on AVI, GOD, SINTACS, and DRASTIC methods), groundwater extraction can be permitted, and industrial activities can be carried out while still observing the industrial area’s environmental conditions.

However, every single method above incorporates a limitation. For the aquifer vulnerability index (AVI) method, certain parameters are ignored, including climate, hydraulic gradient, porosity and water content of the porous media, and sorptive or reactive properties of the layers, which can be contaminant-specific [11]. For the DRASTIC method, the main disadvantage of DRASTIC is the difficulty related to identifying appropriate rating and weight assignments for every parameter [48]. SINTACS methods are “any-aquifer” methods that consider the best number of parameters. For this reason, they produce results that may be used at very high resolution or at large scales (if the dataset is strong enough, as during this case) [32]. The GOD method does not give enough importance to recharge parameters and aquifer permeability. Putting these methods together can improve the vulnerability map.

6. Conclusions

The data analysis results show three categories of regions based on the vulnerability index in Pangkalpinang City (Figure 14). The regions were divided into those with high, medium, and low levels of aquifer vulnerability. Each area requires a different aquifer management and protection system. This differentiation is very important because groundwater resources in the research location are very limited. Areas with high levels of aquifer vulnerability are associated with a high hydraulic protective layer, allowing for vertical fluid percolation to very easily penetrate the aquifer layer. In areas like this, factors that can cause groundwater pollution must be considered. These areas should not be designated for industrial areas and excess groundwater extraction. The upper protective layers in areas with high levels of aquifer vulnerability generally have silty sand lithology. This lithology has a higher hydraulic conductivity value compared to sandstone and schist lithology. Fluid will flow very quickly through this layer to further contaminate aquifers composed of sandstone lithology. Areas with low and very low aquifer vulnerability levels generally have a small hydraulic conductivity protective layer. In these areas, the protective layer is usually composed of sandstone and schist. Groundwater extraction is thus possible with a reasonable extraction pattern. Industrial areas can be built by still paying attention to environmental factors. Areas with medium level aquifer vulnerability generally have a hydraulic conductivity protective layer similar to areas with low and very low aquifer vulnerability levels but with thinner protective layers. In these areas, groundwater extraction can be done with a reasonable extraction pattern. The areas can be used for residential or office purposes with due regard to environmental aspects. This management is important so that anthropogenic waste does not pollute the aquifers in the areas. Management and protection of groundwater will not run properly without proper regulations. Special regulations on proper urban planning to support the sustainability of groundwater are needed.

Data Availability

The data used to support the findings of this study were supplied by Gumilar Utamas Nugraha under license and so cannot be made freely available. Requests for access to these data should be made to Gumilar Utamas Nugraha ([email protected]).

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

All authors contributed to the study of Aquifer Vulnerability: Its Protection and Management—A Case Study in Pangkalpinang City, Indonesia. Gumilar Utamas Nugraha provided VES and AVI analysis and participated in manuscript drafting and interpretation and presentation of the result. Karit Lumban Gaol contributed to software analysis. Priyo Hartanto took part in the evaluation of VES and AVI analysis and also in interpretation and presentation of the result. Hendra Bakti contributed with deep insight into the hydrogeological conceptual model of the study area.

Acknowledgments

The authors would like to acknowledge the Deputy of Earth Sciences, Director of Research Center for Geotechnology Indonesian Institute of Sciences (LIPI). We also would like to thank Dr.Sci. Rachmat Fajar Lubis for insightful discussion of the groundwater system in the study area, Dr. Taat Setiawan for providing the hydrogeological map, and the Department of Energy and Mineral Resources (Dinas ESDM) Bangka Belitung for providing hydrochemistry data.