Study of groundwater vulnerability to contamination using the DRANTHVP method in wates groundwater Basin, Indonesia

Wates Groundwater Basin in Kulon Progo Regency, the Special Region of Yogyakarta, covers an area of 152.67 sq.m. With human occupations continuously growing in the regency, the groundwater in this basin is subject to change. This research was intended to identify the spatial distribution of groundwater vulnerability of Wates Groundwater Basin using the DRANTHVP method, which factors in Depth to Water Table (D), Recharge (R), Aquifer Thickness (A), Nitrate Attenuation Intensity (N), Topography (T), Groundwater Velocity (V), and Pollutant Input Intensity (P). The data were analyzed quantitatively using a weighted tiered approach―the rating was based on to what degree each of these factors contributes to the calculated vulnerability, then overlaid to produce a groundwater vulnerability map. From the analysis, the index values were clustered into three classes of groundwater vulnerability: low 4.48-6.35, medium 6.35-7.39, and high 7.39-7.97. A large proportion of the basin had medium groundwater vulnerability (109.16 sq.m), and less than half of it had low (39.63 sq.m) and high groundwater vulnerability (3.88 sq.m); these severity levels vary depending on the hydrogeological characteristics and land use of a study area.


Introduction
Groundwater vulnerability is an intrinsic or natural characteristic of a groundwater system that depends on the sensitivity of the system to natural and human impacts [1]. It can be differentiated into two: intrinsic groundwater vulnerability and specific groundwater vulnerability. Intrinsic groundwater vulnerability concentrates on physical factors (i.e., rocks, soil, and hydrogeological features) that can naturally protect groundwater from contamination regardless of the pollutant source [2], whereas specific vulnerability factors in potential sources of pollutants, such as the product of human activities like land use and population density [3]. Nowadays, the former is perceived to have not as much significance as the latter because some of its parameters, namely soil medium, the impact of the vadose zone, and water conductivity, are frequently altered as human activities grow [4].
Wates Groundwater Basin (WGB) is situated in Kulon Progo, a regency in the Special Region of Yogyakarta, Indonesia, that currently has the potential to physically develop, especially into industrial estates, office buildings, and hotels, leading to increased population size and rise in groundwater withdrawal to support human activities. To fulfill their daily needs, the local communities extract groundwater in this basin by installing boreholes and dug wells [5]. If improperly controlled, this exploitation will create unstable groundwater conditions and increase threats of pollution coming from, among others, domestic, agricultural, industrial, and development-related activities.
Geomorphologically, WGB is mainly composed of a fluviomarine plain, one precursor of abundant groundwater quantity. On this highly arable morphology, agricultural fields inevitably dominate the area and contaminate the environment due to excessive use of fertilizers. Besides, coupled with the new airport project on the coast of Kulon Progo, the ongoing physical development is highly likely to accelerate and, consequently, alter the groundwater. If all conditions are met, reduced groundwater recharge and ubiquitous groundwater over-extraction in a coastal area can lead to seawater intrusion [6].
Based on the above problems, a groundwater vulnerability study intended to assess the ease with which contaminants are introduced to WGB becomes necessary. One of the most widely used methods to assess groundwater vulnerability is DRASTIC. However, after discovering some significant weaknesses, many prior scholars have modified it into, one of which, the DRANTHVP method [7]. This method was selected as it accommodates a more optimal assessment of groundwater vulnerability to contamination, particularly nitrates.
The research set out to identify the spatial distribution of groundwater vulnerability of WGB using the DRANTHVP method. This assessment lays out groundwater quality indicators that can be useful as a basis for making decisions on short, mid, and long-term groundwater protection, contamination risk, and prevention against pollution.

Research area
The study took place in the Wates Groundwater Basin (WGB), which has a total area of 152,67 sq.m and administratively covers eight districts of Kulon Progo Regency, Indonesia (i.e., Galur, Kokap, Lendah, Pejantan, Pengasih, Sentolo, Temon, and Wates). Astronomically, it is situated at the UTM coordinates 393082-417138 mE and 9117555-9134403 mN. WGB directly borders Purworejo Regency (Bogowonto River) to the west, Bantul Regency (Progo River) to the east, Kokap, Pengasih, and Sentolo Districts (Sentolo Hills) to the north, and the Indian Ocean to the south (Figure 1). These borders follow the research of [12], which delineates this groundwater using both horizontal and vertical boundaries. Horizontal boundaries are extracted from geological and hydrogeological maps, resulting in hydrogeological units, and topographic maps to determine the boundaries of surface water divides. Vertical boundaries are defined from geoelectric and bore data, which provide information on subsurface lithology and aquifer/non-aquifer distribution and dimension.
Based on the 1:100,000 Geologic Map of Yogyakarta [8], WGB geologically comprises two formations, namely the dominant Alluvial Deposits (Qa) and the Young Merapi Volcano Deposits (Qmi) in the east of the basin or the floodplain of Progo River. The regional stratigraphic column of Kulon Progo [8] shows that Qa is the youngest formation formed in the Quaternary Holocene Epoch. Meanwhile, Qmi was formed before Qa in the Quaternary Pleistocene Epoch. There is one group of rocks that do not belong to either formation, namely a complex of beach ridges/beaches and dunes composed of Holocene-aged sands [9]. Geomorphologically, WGB is composed of two landscapes: fluvial-marine and marine-aeolian. The fluvialmarine landscape has Fluvio-marine Plains (Fm) and Natural Levees and Floodplains (F), while the marineaeolian landscape includes Old Beach Ridges (M2), Sand Dunes and Swale Complex (E), and Young Beach Ridges (M1) [10].

DRANTHVP method for groundwater vulnerability assessment
The DRANTHVP method [7] was developed and first implemented in the west of the Jilin Province, China, where the majority of the area belongs to the Songnen Basin. It is a development of DRASTIC that optimizes several parameters while maintaining the basic structure of the original method and combines two main factors of groundwater vulnerability, namely hydrogeological features and human activities. For these reasons, DRANTHVP has become a method of assessing specific groundwater vulnerability.
Some of the modifications applied to the original DRASTIC aim for optimization by replacing qualitative parameters into quantitative ones and adding another parameter, better selection of the assessment scale of every parameter, and the application of Projection Pursuit Dynamic Clustering (PPDC) model for parameter weighting. This method consists of research variable determination and data acquisition and processing techniques. The research variables were eight groundwater vulnerability parameters factored in the DRANTHVP method, namely Depth to Water

Data acquisition and processing techniques
The primary data were acquired from the measurements of depth to the water table and groundwater velocity and land use validation. The secondary data included recharge, aquifer thickness, nitrate attenuation intensity, slope gradient, hydraulic resistance, and pollutant input intensity.
For the depth to water table, this parameter was measured at the water wells selected by systematic random sampling. First, the entire area of WGB was divided into 1000 ×1000 m 2 grids, and then the water well samples were selected randomly from each of the resulting 201 grids. Eighteen years of rainfall data, from 2002 until 2019, at the nearest station were collected from the Main Station of Serayu-Opak River Basin (BBWS-SO). Aquifer thickness was obtained from the results of geoelectrical measurements conducted by previous researchers [10][11][12]. As for nitrate attenuation intensity, it was extracted from soil type and texture data in the study area, as presented in the Semi-Detail Soil Map of Kulon Progo Regency Sheet No. 1408-21 Scale 1:50,000 produced by the Main Station of Research and Development on Agricultural Land Resources (BBSDLP). Slope gradients were derived from the 1:25,000 Indonesia Topographic Map (RBI). The test water wells were selected by purposive sampling technique based on landform as the sampling unit. This sampling nature complements data from previous research [12]. Hydraulic resistance was extracted from depth to water table as the thickness of unsaturated zones, i.e., the distance from the ground surface to the water table or saturated zone, and from hydraulic conductivity of unsaturated zone [13]. A pumping test was conducted to determine the hydraulic conductivity. As for the pollutant input intensity, this data was extracted from the land use map (the 1:25,000 Indonesia Topographic Map) and validated indirectly using the interpretation results of the Google Earth satellite image and directly in the field.
Groundwater recharge was obtained by calculating data on annual rainfall, annual transportation, crop evapotranspiration, soil texture, and land use (the type of vegetation) using the equation below: where P is rainfall, ET is crop evapotranspiration, and Ro is surface runoff [7]. The evapotranspiration was computed using the Turc Lungbein formula, then multiplied by crop coefficients (Kc), while the surface runoff was estimated by the Soil Conservation Service (SCS)-Curve Number (CN) method as a function of land use and soil type [14]. Geoelectric data of aquifer thickness were processed in IP2WIN to produce resistivity values of detected underground materials and validate drill data, and the properties of these materials tell the type, depth, and thickness of the aquifer. The range of nitrate attenuation intensity represents differences in nitrate content in groundwater, porosity, and percent organic matter between existing soil texture. Hydraulic resistance is the ratio of the thickness of unsaturated zones to hydraulic conductivity, which was calculated using the equation below [15]: where c is hydraulic resistance, d is the thickness of unsaturated zone (m), and k is hydraulic conductivity (m/day). Darcy's law was adopted to calculate groundwater velocity. The equation used was as follows: where V is groundwater velocity (m/day), K is permeability (m/day), dh/dl is hydraulic gradient (m/m), and ne effective porosity. The hydraulic gradient was defined by first creating a groundwater contour map and cross-section for every landform unit, while the effective porosity was determined following [16]. Land use data were classified into a certain range of pollutant input intensity respective to the amounts of pollutants entering the groundwater described in [17]. Every parameter was processed in the ArcGIS 10.3 program by geostatistical interpolation and then classified according to the requirements of the DRANTHVP method, as summarized in Table 1.

Data analysis technique
Groundwater vulnerability assessment in WGB ultimately aims at implementing the DRANTHVP method. The data were analyzed quantitatively using a weighted tiered approach, and the rating was based on to what extent the parameters contribute to the resulting vulnerability. Because every parameter was given a numerical rating, then a weighting factor based on its contribution, it had different weight values.
Eight maps, each representing the DRANTHVP parameter, were overlaid in the ArcGIS 10.3 program to produce and map groundwater vulnerability index using the equation below [7]:

Groundwater vulnerability parameters
Depth to water table (D) data were measured at 156 bore wells. Data processing showed that WGB had four out of the ten predetermined classes of depth, each of which has a different response to the movement of pollutants. A deeper water table means that pollutants require a longer time to reach aquifer; hence, lower groundwater vulnerability. In WGB, depths to water table were mostly distributed into the range of 2 m to 4 m, followed by 0 m to 2 m, 4 m to 6 m, and 6 m to 9 m. Spatially, depths to water table at the range of 2-4 m were distributed almost evenly, from the middle to the west of WGB. The second-largest class, 0-2 m, was found from the middle to the east of the basin. The other class, 4-6 m, were distributed in part of the southern coast and the north of the basin, including the groundwater recharge area in Wates and Pengasih Districts. Meanwhile, water tables located 6-9 m from the ground surface were only detected in several locations in Sentolo District directly adjacent to Sentolo Hills. In this case, since Sentolo Hills is a topographical factor that affects this parameter, it was categorized as a groundwater recharge area.
Groundwater recharge is the difference between total rainfall and cumulative loss due to surface runoff and effective evapotranspiration [14]. The larger the groundwater recharge of an aquifer, the more pollutants the groundwater carry and the higher the groundwater vulnerability. WGB received a regional rainfall of 1793 mm/year and experienced water losses by evaporation (shaped by the type of land use and vegetation cover) and runoff (depending on the type of land use and soil). Out of 10 classes of groundwater recharge, data analysis revealed that the research area had four. Groundwater recharge of 130-250 and 115-130 mm/year was distributed in most parts of the south coast. Physically, this coast was composed of sand-textured regosols and was dominantly used for agricultural practices, namely, plantation and dry cultivated area. Groundwater recharges within the ranges of 75-85 and 85-95 mm/year were found in most parts of areas composed of sandy loam and sandy clay-textured alluvial soils and in varying land uses, including agricultural and built-up areas. Coarse soil textures have a large porosity and, consequently, high drainage value, meaning that it creates an environment that produces a large quantity of water input to the soils. Land use determines the amount of runoff, e.g., built-up land tends to have a larger runoff volume than agricultural land, reducing the portion of rainfall that seeps into the soil (groundwater recharge).
In DRANTHVP, aquifer thickness replaces the qualitative parameter of the aquifer media in DRASTIC that, according to many previous scholars, fails to reflect the spatial dimension of groundwater storage space. Besides, it plays a crucial role in determining the ability of groundwater to dilute nitrate [7]. A thicker aquifer can dilute more nitrate in the groundwater, resulting in low vulnerability [18].
Interpretation of subsurface materials from the geoelectrical sounding results and bore data showed that the aquifers in WGB had varying thicknesses. The aquifer thickness was mainly in the range of 10-45 m, which was This aquifer was included in the Fluviomarine Plains aquifer system composed of two layers: (1) Alluvium deposits alternating with marine clay and (2) marine-clay alternating with marine sand [10]. The thickest aquifer (45-60 m) was found on the south coast-belonging to the aquifer system of the Sand Dunes Complex-and composed of a bowlshaped layer of marine sand, which facilitates rainwater collection [10]. In the east, a 30m-thick local aquifer was formed of Young Merapi Volcano Deposits and due to the fluvial process of the Progo River and.
Nitrate attenuation intensity or the ability to reduce nitrate content is a quantitative parameter that replaces the soil texture in DRASTIC [7]. The consideration is that soil media is considered to significantly affect the ability of pollutants to transport vertically to the vadose zone. That is because there are differences in groundwater content, porosity, and organic matter content in soil media that can affect the attenuation of nitrates that enter the groundwater [19]. This condition modifies the dilution process of nitrate by groundwater; the coarser the soil texture, the smaller the nitrate attenuation intensity.
Nitrate attenuation intensity at a relative range of 10-30 was sand-textured regosols distributed on the south coast. Nitrate attenuation intensity at a relative range of 19-20 was sandy clay-textured soils distributed in most parts of the basin, i.e., in the middle. In the study area, soils with sandy clay texture are alluvial soils, gleisol, and cambisol. Claytextured gleisols were found in the north of the basin, with a nitrate attenuation intensity of 20-25.
Slope gradients determine runoff formation and rainwater infiltration process, which potentially carries pollutants into the aquifer. The higher the slope gradient, the larger the runoff volume and the smaller the quantity of water seeping through the soils; hence, lower groundwater vulnerability. A large proportion of the study area had 0-1% slopes, meaning that most of the basin lies in the middle of the fluviomarine plains and has flat topography. Slopes with a gradient of 1.5-2.5% were detected on the south coast with nearly flat topography and dunes. Slopes with gradients of 2.5-5%, 5-7.5%, and 7.5-10% were found in the north and directly adjacent to the Sentolo Hills, creating a flat to undulating topography.
Hydraulic resistance is a substitute parameter for the impact of unsaturated or vadose zones in the DRASTIC method because these zones are considered providing a small discretization in a porous aquifer (mainly due to dominant clay-sand-gravel soils) and, hence, a qualitative parameter. Hydraulic resistance of the vadose zone can be used to measure the relative time of transit of pollutants from the surface through the vadose zone and then into groundwater [15]. It is a ratio of the thickness of unsaturated zones to hydraulic conductivity. The higher the hydraulic resistance, the longer the transit time of pollutants before reaching the unsaturated zone; hence, lower groundwater vulnerability.
Like the depth to water table (D) data, the thickness of unsaturated zones in the hydraulic resistance calculation was measured at 156 dug wells during the dry season. Unsaturated zones with a thickness of 6-9 m were distributed on the south coast, with marine sand deposits as the constituent rocks and hydraulic conductivity of 10 m/d. Meanwhile, thinner unsaturated zones, 0-6 m, were found in the majority of the basin are (in the middle), while thicker ones, 9-13 m, were in the north. Both classes were composed of alluvium materials interspersed with clay and fine sand, resulting in hydraulic conductivity of 10 -3 m/d.
Hydraulic resistance is directly proportional to the thickness of saturated zones but inversely proportional to hydraulic conductivity. The results of hydraulic resistance at CAT Wates The majority of WGB was in the range of 2500-4000 years, followed by 0-1000 years and 1000-2500 years that stretched from the middle to the south (the closer the aquifer to the south, the smaller the hydraulic resistance). Low hydraulic resistance in the middle of the basin was attributable to the relatively shallow unsaturated zone, 0-4 m. On the contrary, high hydraulic resistance in the south was associated with the constituent materials of the unsaturated zone, which were sand. As for the other ranges of hydraulic resistance, i.e., 4000-5000, 5000-6000, and 6000-7000 years, they covered increasingly smaller areas in the northern part of the basin. Such high hydraulic resistance was caused by deeper water table and clay materials in the unsaturated zone, creating a layer that has a low ability to transmit water (low hydraulic conductivity).
Groundwater velocity is a substitute parameter for hydraulic conductivity because previous scholars [20] found a negative correlation between hydraulic conductivity and nitrate but a positive correlation between nitrate and groundwater velocity. A higher groundwater velocity indicates a faster distribution of pollutants into the aquifer, or in other words, a higher groundwater vulnerability. This parameter was computed using Darcy's law [16] by factoring in hydraulic conductivity, hydraulic gradient, and effective porosity of every landform unit.
The calculation results showed that the highest hydraulic conductivity, 50.62 m/d, was in the Natural Levees and Fluviomarine Plains, followed by 5.80 m/d in Beach Ridges, 4.7 m/d in Fluviomarine Plains, and 1.32 m/d in Sand Dunes. Hydraulic conductivity largely depends on geological characteristics. The highest hydraulic gradient was found in Fluviomarine Plains in the north of the Wates District, which functions as a discharge area. Generally, WGB had very small hydraulic gradients because most of the areas overlying it were flat. The effective porosity of sand is 0.43, while that of clay is 4.2 [16].
The calculation results showed that, on average, the groundwater velocity in Beach Ridges was 0.077 m/d, Sand Dunes Complex 0.0096 m/d, Fluviomarine Plain 0.0513 m/d, and Natural Levees and Floodplains 0.2997 m/d. Based on the class range proposed in the DRANTHVP method, all of these values were included in the same class, 0-1 m/d (the smallest score), meaning that from the aspect of groundwater velocity, WGB has very low groundwater vulnerability.
Land use is a specific parameter in groundwater vulnerability and is sensitive to the human factor. Since land use is a qualitative parameter, every type of land use was assigned with values that were relatively different from pollutant input intensity and according to the variety of sources of pollutants in the soil [17]. Pollutant input intensity had the highest weight value in the groundwater vulnerability assessment, implying its strong effect on groundwater vulnerability. The higher its range, the higher the groundwater vulnerability score. In WGB, there were seven (out of 10) classes of this parameter. Accordingly, the ascending order of land-use types based on pollutant input intensity was as follows: industrial/commercial areas, urban areas, dry cultivated land, paddy field, water body, grassland, and vacant land.
The most widely distributed pollutant input intensity was 250-300 (dry cultivated land), followed by 100-250 (paddy) in nearly all of the basin area, particularly in the middle. In urban areas, as marked by vast settlement areas and many places of activities, it varied between 250 and 500, which also was the most widely distributed range in Wates District, the center of Kulon Progo Regency. The highest pollutant input intensity was 500-800, which covered industrial/commercial estates like buildings/factories (industrial) and the project of the new Yogyakarta International Airport (commercial). Meanwhile, the least distributed one was the lowest class, 1-50, found in grassland, vacant land, and water body.

Groundwater vulnerability to pollution
Groundwater vulnerability to pollution is the overlay of eight parameters computed in the DRANTHVP method [7]. The results showed that WGB had varying index values of groundwater vulnerability, from 4.953-7.704, which were categorized into three classes (out of 5-scale classification), namely low, medium, and high ( Figure 3). The low groundwater vulnerability covered an area of 39.63 sq.m or 25.9 6% of the total area and was identified in most areas of Fluviomarine Plains and Beach Ridges, or administratively in Temon, Lendah, and Kokap Districts.
Meanwhile, a large proportion of WGB (109.16 sq.m, 71.50 %) fell into the category of medium groundwater vulnerability (index= 6.35-7.39). This class was distributed in almost the entire landforms in the study area, namely Fluviomarine Plains, Natural Levees and Floodplains, Sand Dunes Complex, and Beach ridges. Except for Sentolo, all districts within the WGB area had medium groundwater vulnerability. Only the smallest share of it (3.88 sq.m, 2.54%) had a groundwater vulnerability index of 7.39-7.97, which was categorically high. This class of vulnerability was distributed in the west, including the Sand Dunes Complex, Beach Ridges, and Fluviomarine Plains in Temon and Wates Districts.

Fig. 3 Groundwater Vulnerability Map of Wates Groundwater Basin
Groundwater vulnerability levels in the study area are dependent on existing hydrogeological conditions, as reflected by the seven parameters computed in the DRANTHVP method. In this case, hydrogeological variation is the product of different rock formations and landform units. Unlike the passive hydrogeologic characteristics, land use-a representation of sources of pollutants-is a rather dynamic parameter. Groundwater vulnerability is likely to worsen if the underlying hydrogeological conditions become more sensitive to pollution and land-use change as the highest contributors of pollutant loads.