Vulnerability methods in hard rock formation as a basis for groundwater risk assessment – from resource to source

Groundwater in a hard rock formation is most endangered at places where a potential source can discharge contaminants that can reach the saturated zone of an aquifer. In these circum -stances, an essential tool for groundwater protection is the contamination risk map. This map is based on the integration of two maps: a hazards map, i.e., map of potential sources of contami­ nation and a vulnerability map. The selection of a proper vulnerability method is an important task since the resulting vulnerability map can significantly impact the final risk of contamination map. The most appropriate method for groundwater vulnerability assessment was considered in the case study of Tara National park, in Western Serbia. The four commonly used methods were applied to assess the intrinsic vulnerability maps: DRASTIC, EPIK, PI and COP. All the applied methods resulted in different vulnerability maps in assessing the degree of vulnerability, conse quently influencing the groundwater contamination risk maps. The applied research presents an example of how contamination risk should be assessed in a specific area. Comparison of the results obtained for the area of Tara National Park indicates the preference of the PI method as a well­balanced method, taking into account all the specifics of the study area. A detailed analysis of the assessed risks in the catchments of the existing sources was also conducted to indicate probable sources of contamination and confirm the de gree of accuracy of the created vulnerability and risk maps. The conducted research empha-sizes the necessity to adopt a clear conceptual hydrogeological model and to apply several methods simultaneously to determine the optimal one for each individual area.


INTRODUCTION
Groundwater in hard rocks, especially carbonates, can be very vulnerable to contamination, so protection is important. Although its water quality is generally good, the existing Serbian regulations on groundwater protection do not sufficiently consider the specificity of groundwater circulation in these rocks, and the current state may cause possible contamination (OFFI-CIAL GAZETTE OF RS, 2008). This shortcoming of groundwater management is particularly pronounced in national parks, which are recognized as groundwater-dependent ecosystems with rare and protected species.
For groundwater protection of carbonate aquifers, a general framework for vulnerability and risk mapping was developed by the Working Groups of COST Action 620 (ZWAHLEN, 2004). This European Approach is based on an origin-pathway-target model, which applies to both resource and source protection. The origin is the location of the source of contamination release. The pathway represents the direction of contaminant transport, from the source-origin of contamination toward the target area, which has to be protected , whether it is the protection of resource or source (GOLDSCHEIDER, 2004).
The vulnerability map is essential for groundwater protection and land use planning for protected areas such as national parks. Many different methods of vulnerability assessment have been proposed and the overview of the most current methods and their specificities can be found in the works of several authors (VRBA & ZAPOROZEC, 1994;GOGU & DISSAGRUES, 2000; Vulnerability methods in hard rock formation as a basis for groundwater risk assessment -from resource to source Vladimir Živanović 1 , Igor Jemcov 1 and Veselin Dragišić 1 1 University of Belgrade, Faculty of Mining and Geology, Department of Hydrogeology, Đušina 7, 11000 Belgrade, Serbia;(vladimir.zivanovic@rgf.bg.ac.rs) doi: 10.4154/gc.2022.23 Abstract Groundwater in a hard rock formation is most endangered at places where a potential source can discharge contaminants that can reach the saturated zone of an aquifer. In these circumstances, an essential tool for groundwater protection is the contamination risk map. This map is based on the integration of two maps: a hazards map, i.e., map of potential sources of contami nation and a vulnerability map. The selection of a proper vulnerability method is an important task since the resulting vulnerability map can significantly impact the final risk of contamination map. The most appropriate method for groundwater vulnerability assessment was considered in the case study of Tara National park, in Western Serbia. The four commonly used methods were applied to assess the intrinsic vulnerability maps: DRASTIC, EPIK, PI and COP. All the applied methods resulted in different vulnerability maps in assessing the degree of vulnerability, consequently influencing the groundwater contamination risk maps. The applied research presents an example of how contamination risk should be assessed in a specific area. Comparison of the results obtained for the area of Tara National Park indicates the preference of the PI method as a wellbalanced method, taking into account all the specifics of the study area. A detailed analysis of the assessed risks in the catchments of the existing sources was also conducted to indicate probable sources of contamination and confirm the degree of accuracy of the created vulnerability and risk maps. The conducted research emphasizes the necessity to adopt a clear conceptual hydrogeological model and to apply several methods simultaneously to determine the optimal one for each individual area. DREW & DUNE, 2004;FOSTER et al., 2013;IVAN & MADL-SZONYI, 2017). Intrinsic vulnerability assessments are based on the advection, attenuation, and relative quantity of conservative contaminants (GOLDSCHEIDER & POPESCU, 2004). Thus, the vulnerability accounts only for the hydrogeological characteristics of the system and may contain numerous uncertainties (NA-TIONAL RESEARCH COUNCIL, 1993) so that different vulnerability maps can be obtained for the same area by different authors or the same authors when applying different methods. A step forward was made by the European Approach, although it is non-prescriptive and could be adapted to a method that respects the specifics of the hydrogeological characteristics of the location. Finally, this approach is not solely formed for karst aquifers but can also be applied to other aquifers (ZWAHLEN, 2004).
Hazards are potential sources of contamination that usually originate from anthropogenic activities. They are assessed according to their harmfulness associated with their toxicity and quantity (KETELAERE et al., 2004). Hazard classification has been done based on land use characteristics and the obtained hazard map for nature protection areas can be an effective tool for land use planning.
Finally, Vulnerability and Hazard assessments alone are not sufficient for sustainable groundwater management since groundwater protection requires a comprehensive Risk analysis. Risk evaluation further produces a quality basis for groundwater protection. The risk assessment process depends on the following three elements: the hazard -equivalent to the origin, the intrinsic vulnerability -equivalent to the pathway, and potential conse-quences of groundwater contaminations -equivalent to the target. Following these principles, groundwater contamination risk maps can be obtained by analyzing two factors: hazard index and groundwater vulnerability index (HÖTZL, 2004). Also, it is necessary to emphasize that these maps present a valuable tool for planning a groundwater monitoring network and ensuring longterm groundwater quality (FOSTER et al., 2002).
The European approach to the risk analysis was successfully applied to many examples, mainly in Europe, and the majority of them comprise the methodology of the Risk assessment, developed in the framework of COST Action 620. Even though various methods have been applied to determine the intrinsic vulnerabili ty for interpretation of the European approach including PI ( GOLDSCHEIDER et al., 2000), VULK (CORNATION et al., 2004), LEA (DUNN, 2004), COP (VIAS et al., 2006), The Time-Input (KRALIK & KEIMEL, 2003), but none of these considered represents the only possible interpretation of the European approach (GOLDSCHEIDER, 2004). On the Sierra de Libar, two vulnerability methods (PI and COP) were tested for Risk analysis and the results obtained are consistent . NGUYET & GOLDSCHEIDER (2006) proposed a simplified vulnerability and risk mapping methodology for application in data-poor environments. Alternatively, RAVBAR & GOLDS-CHEIDER (2007) developed a more detailed methodology of vulnerability and contamination risk mapping -the Slovene Approach, which includes the possibility of integrating surface and groundwater protection by applying an extension of the COP method for source vulnerability mapping. The COP method has been applied to evaluate the intrinsic vulnerability (JIMÉNEZ-MADRID et al., 2010) in Spanish transboundary basins to obtaining contamination risk taking into account chemical analyses of carbonate aquifers.
Considering the specifics of the hydrogeological characteristics of each location, the choice of the intrinsic vulnerability methods remains an open question that should be considered through a comparative analysis of several well-established and most used methods. In the case study of the Tara National Park with a large extent of carbonates, several methods of intrinsic vulnerability were comparatively analyzed: DRASTIC -as a widely used and non-specific method for karst aquifer (ALLER et al., 1985); EPIK -as the first and widely used method for karst areas (DOERFLIGER & ZWAHLEN, 1997) and the two most commonly used methods in the European approach PI (GOLD-SCHEIDER et al., 2000) and COP , to obtain a reasonable basis for further risk mapping. Tara National Park represents a good example, mainly because of its importance as needing the highest form of environmental protection. Potentially, any activities could deteriorate the present quality of groundwater, and preventative measures must be applied, particularly at the locations rated with the highest risk of groundwater contamination (ŽIVANOVIĆ et al., 2014).

APPLIED METHODS
Based on the framework of COST Action 620 and the proposed model, the groundwater contamination risk is evaluated through the potential sources of contamination as an origin, the vulnerability of the aquifer as a pathway, and groundwater which has to be protected as a target (DALY et al., 2002). Three essential steps need to be taken to produce a groundwater contamination risk map: 1. Vulnerability mapping; 2. Hazard mapping; 3. Risk mapping -risk assessment. All thematic maps were elaborated using a Geographical Information System -Esri ArcGIS (ESRI, 2014).

The methods of the Intrinsic vulnerability mapping
The European Approach uses three factors in assessing intrinsic vulnerability for resource protection mapping: overlying layers (O), the concentration of flow (C), precipitation regime (P), while karst network development (K) may be used for source protection mapping (DREW & DUNE, 2004). If we take into account the fact that vulnerability is often considered to be a qualitative, nonmeasurable concept rather than a quantitative property, the resulting maps can be different or even contradictory. Thus, there is a need to examine vulnerability concepts from a quantitative point of view (BROUYÈRE, 2004). Since the approach is nonprescriptive for this study, the four most frequently applied methods were used: DRASTIC, EPIK, PI, and COP. Some of the selected methods were purely developed for carbonate karstified aquifers, such as the EPIK method. In contrast, the DRASTIC, PI, and COP methods can be applied to various types of aquifers as well as to fractured hard rocks. The PI method gave the most significant input to the European Approach, and the COP method similarly provides methodological tools for karst (DREW & DUNE, 2004).
The DRASTIC method (ALLER et al., 1985) represents one of the first and frequently applied methods for assessing groundwater vulnerability for all aquifer types. A vulnerability index calculation is based on seven different parameters: depth of groundwater (D), net recharge (R), aquifer media (A), soil media (S), topography (T), the impact of the vadose zone media (I), and the hydraulic conductivity of the aquifer (C). For each parameter, specific maps are created. In order to emphasize the influence of specific parameters, weight factors are included in the equation:

DRASTIC Index = 5·Dr+4·Rr+3·Ar+2·Sr+Tr+5·Ir+3·Cr
where Dr, Rr, Ar, Sr, Tr, Ir, and Cr are values of parameters for a given cell. Depending on the value of this vulnerability index, 6 vulnerability classes are distinguished (from very low to very high).
The EPIK method is a specific method for karst groundwater vulnerability assessment (DOERFLIGER & ZWAHLEN, 1997) created in the early stages of development of the European Approach. A protection index F is calculated based on four parameters which form an acronym: E -Epikarst, P -Protective cover, I -Infiltration Condition, and K -Karst network development. Karst groundwater vulnerability assessment is calculated according to the proposed equation: F = αE i + βP j + γI k + δK l where α, β, γ, and δ are weight coefficients (3, 1, 3 and 2). The protection factor ranges between 9 and 24 and is divided into four groundwater vulnerability classes.
The PI method was developed within the COST 620 project (ZWAHLEN, 2004). This method is applicable to all types of aquifers with a particular focus on karst aquifers. The acronym PI represents two parameters that are taken and thoroughly processed (GOLDSCHEIDER et al., 2000): a protective factor P (protective cover) and infiltration conditions (I). The P factor is based on the GLA method (HOELTING et al., 1995), while the I factor describes the degree to which the protective cover is bypassed regarding the infiltration conditions. The vulnerability index, π, represents the multiplication of these two factors and the final map is created by categorizing the vulnerabilities into five classes.
The COP method  represents the approach for determining intrinsic resource vulnerability mapping (DALY et al., 2002). Vulnerability is determined based on three factors: Concentration of flow (C), Overlying layers (O), and Precipitation (P). The C factor represents the sum of factors that impact infiltration conditions. The O factor depends on the hydrogeological properties of the rocks and sediments above the saturated zone. The P factor represents the amount and intensity of precipitation. Factors O and C are quantified in a similar but slightly simplified way compared to the PI method (GOLDSC-HEIDER et al., 2000). A COP index is calculated by multiplying the results of these three factors. The final vulnerability value varies from 0 to 15, and the final map presents the results of COP index classification divided into five classes.

Hazard mapping
Based on land use characteristics and activities that threaten groundwater quality, an inventory of identified potential contaminant sources was carried out, according to the guidance of COST Action 620 (KETELAERE et al., 2004). Map creation was based on the character of contamination sources (punctual, linear, or diffuse). Following this recommendation, a weight coefficient was introduced in assessing harmfulness and assigned to each hazard.
The hazard index (HI) that describes the degree of harmfulness concerning the type of contaminant that could be released and their likelihood of release, has been calculated for each indicated hazard. The hazard index (HI) obtained is based on three factors, as shown in the following equation (KETELAERE et al., 2004): HI -Hazard index (possible ranges from 0 as a minimum to 120 as maximum scores), H -Weighting factor (defined between 10 and 100), Q n -Ranking factor (ranges between 0.8 and 1.2), R f -Reduction factor (ranges from 0 to 1). The Hazard map is presented based on the potential harmfulness for each assigned hazard, where the appropriate colour has symbolized each hazard.

Risk mapping
A groundwater contamination risk map represents the likelihood that groundwater in a specific aquifer will become contaminated to an unacceptable level (MORRIS & FOSTER, 2000). This concept considers the interaction between contaminant infiltration and the vulnerability of the aquifer.
Groundwater contamination risk maps can be obtained by analyzing two factors: the hazard index and the groundwater vulnerability index (HÖTZL, 2004). The first one relates to hazards, the existing potential sources of contamination. Hazards are usually mapped in a GIS environment, and different hazard indices are assigned depending on the hazard type and the possibility for a contaminant to be discharged on the ground surface (KETE-LAERE et al., 2004).
The risk of contamination is calculated by combining the vulnerability and hazard assessment. A risk intensity index is calculated using the vulnerability index obtained by different methods (π) and previously defined hazard indices (HI) according to the following formula (HÖTZL, 2004):

Risk intensity maps for groundwater source protection
The hazards and their potential contamination, particularly when it comes to the point type of hazards, usually affects only a part of the groundwater body related to the catchment area. Contamination depends mostly on the horizontal components of groundwater flow when it reaches the saturated zone. Thus, the risk is related only to the potentially contaminated part and not to the entire groundwater resource . Risk assessments should only refer to the parts of groundwater which are endangered by a potential impact. The study area needs to be delineated into sub-catchments based on topography, structural-geological settings, and hydrogeological properties to assess the flow direction.

Geostatistical analysis
All maps have been rasterized to have the same extent and discretization (cell size 25 x 25 m) to compare vulnerability and groundwater risk contamination maps. Also, it was required to reclass vulnerability maps to enable comparative analysis of groundwater vulnerability (NEUKUM & HÖTZL, 2005). The maps prepared in this manner are further converted to a table, where each row is filled with central coordinates of each raster cell and with vulnerability, hazard, and risk values, obtained by every applied method. Statistical procedures were applied to the obtained table for comparative observation of the results of the vulnerability methods and, accordingly, the risk methods. In order to compare the risk of groundwater contamination, all three columns in the table that take into account the nature of the contamination sources: punctual, linear or diffuse (NONER, 2004) are sublimated (fused) into one column in the table, taking into account the highest hazard indexes (worst-case scenario principle).

An additional tool for groundwater contamination risk analysis
In order to support contamination risk assessment, a sampling campaign for physicochemical properties of groundwater was carried out. Sampling campaigns were performed after high precipitation events that affected the main springs' discharge. Under these water conditions, the worst springs water quality was expected. During this campaign, 21 samples from springs were collected. Some parameters were determined in situ -pH, electric conductivity (EC), temperature T, redox potential (ORP), total dissolved solids (TDS), and dissolved oxygen (DO) using a portable Hanna HI981984 Multiparameter Meter. All samples were filtered, stored, and transported for further analysis in the certified laboratory (Institute for chemistry, technology and metallurgy IHTM, Belgrade) for Electrometry for Atomic adsorption spectrometry (Ca 2+ , Mg 2+ , Na + , K + ), UV-VIS spectrophotometry (NH 4 -, Fe tot , Mn tot , P), ion chromatography (NO 2 -, NO 3 -, Cl -, SO 4 ), volumetric titration within 24 h of sampling (CO 3 2-, HCO 3 -, Cl -, KMnO 4 ) and microbiology (for Perućac, the largest karst spring). Results were presented on a Piper diagram.
Groundwater chemistry can be applied as an indicator of the potential risk of groundwater contamination, but it is necessary to change the target from the groundwater surface to the source vulnerability, which includes the horizontal pathway in the aquifer. Therefore, catchments of specific sources needed to be defined to analyze the pressure of detected hazards that may indicate contamination risk to a particular source -spring. For this purpose, the impact of the presented risks in the catchments was analyzed on several characteristic examples -springs (karstic and non-karstic).

CHARACTERISTICS OF THE STUDY AREA
Tara National Park represents a natural habitat for unique species of flora and fauna, and it covers an area of 192 km 2 . It is situated in western Serbia, on the border with Bosnia and Hercegovina with the Drina river, which forms the country's natural western border dividing both countries (Figure 1).
The research area belongs to the Dinaric orogenic system and is situated in the zone of direct contact of two large tectonic units -the continental Drina-Ivanjica block in the East and the Western Vardar ophiolite unit in the West. The separations between them appeared due to sequence trusting, significantly modifying the original geometry (SCHMID et al., 2008). The main structural element between those two zones is the Konjska river strike-slip fault zone, oriented NW-SE. There is an extensive thick-bedded to massive organic limestone of Triassic age in the eastern part, representing a well-developed karst aquifer. The basement of this aquifer is claystone, fine-grained sandstone, and schist of Lower Triassic and Palaeozoic age. In the western part, Cretaceous limestone and marl (moderately karstified), transgressively overlie the Jurassic harzburgite, gabbro, and diabase (fissured or non-aquifer).
The geological settings have conditioned the landscape so that most of the area is in the form of plateaus with steep slopes in the north, where the main drainage points of the karst aquifer are located (Perućac, Lađevac, Rača -springs S-1, S-2, and S-3 on Fig.1, respectively). Next to the dolines, dry valleys are the characteristic landforms on the plateau, predetermined by faults (e.g., Mitrovac fault) and a few uvalas ending with a ponor or ponor zone (Mitrovac ponor -P-4 and Ljuto polje ponor zone -P-2). On the north slope, caves formed at three different levels can be observed (approx. 900, 650, and 500 m asl) as a consequence of the intense incision of the Drina River, which led to the development of the karst process to base level (250 m asl -Perućac spring S-1). During construction of the hydrotechnical tunnel (8 km long) through the Tara massive, from the Drina river 610 m upward to another reservoir -the Beli Rzav river catchment, a few different zones were registered. The upper parts of the carbonate sequence (below the epikarstic zone) were intensely karstified but were often filled with secondary material (clay bauxite and sandstone). An intensively karstified zone, with open conduits, was detected at a depth of 150-200 m below ground surface. In the deeper parts, karstification is exclusively related to the fault zones and reaches the base level.
The mean annual air temperature is 8.3 ºC. Mean annual precipitation is about 1000 mm (with about 40 % being snow), and the effective infiltration as a percentage of total precipitation is estimated to be 60 % (based on groundwater budget modeling), estimated for the main drainage point -Perućac spring S-1, and slightly less (55%) for the Rača S-2 and Lađevac S-3 springs. According to the time series analysis, this spring has a large storage capacity and is characterized by a well-structured system, with significant attenuation of the impulse response with base flow index above 0.8 for Perućac spring -S-1 (JEMCOV & PETRIČ, 2009) and 0.56 for the Rača spring -S-2. This is a consequence of a significant lateral flow in the epikarst and lower permeability Figure 1. A geological map of Tara National Park (OLUJUĆ & KAROVIĆ, 1985, modified). al -alluvial deposit -intergranular aquifer; d -deluvial deposit -intergranular aquifer; K 2 1,2 -limestone (micrite, marly micrite, and marls) -karst aquifer; K 2 1 -finegrained sandstone (arenite) and claystone, marls, marly limestone (micrite) -fissured aquifer; J 2,3 -diabase and chert (fissured aquifer); vββ, ββ, σ -gabbro-diabase, diabase, harzburgite -fissured aquifer, T 3 -limestone (sparites) -karst aquifer; T 2 -massive limestone (microsparite and micrite) -karst aquifer; 2 T 1 -sandstone and marlstone -fissured aquifer; T 1 -shale -fissured aquifer; C -Sandstone, schists -non-aquifer. due to infilling in the upper karstified zone. Most of the research area is covered by forest (over 85%), although carbonate rocks are covered with a thin layer of skeletal soil with a relatively high content of humus (JEMCOV et al., 2007).
Special attention was paid to the hydrogeological function of the uvalas. According to geophysical research and the explored open pits, the layer of Upper Cretaceous fine-grained sandstone (arenite), marl, claystone, and marly limestone (micrite) transgressively superimposed on the mostly karstified Triassic rock has a thickness of about 5-30 m in the uvalas. This layer is partially covered by sediments of diluvium (3-5 m thick) with a varia ble hydraulic conductivity of about 10 -5 -10 -7 m·s -1 (KREŠIĆ, 1984, JEMCOV, 2014. Two tracer dye tests conducted at different locations show the full complexity of this hydrological system (JEMCOV et al., 2011). The first location of ponor P-1, predetermined by the fault zone at the contact of two carboniferous Triassic rocks. Injected uranine appeared in the Perućac spring -S-1. The maximum linear flow velocity, calculated from the moment of first tracer detection, was about 107 m/h, classifying it in a group of "very high" vulnerability. The second test was conducted at the largest ponor (about 20 l·s -1 ) in the ponor zone Ljuto Polje -P-2, at the contact with the less permeable Cretaceous arenite, marl and claystone with Upper Triassic carbonates, during the period of snow melt. Though all the important springs were taken into account as observation points (Perućac, Rača, Lađevac, Solotuša) the tracer was not detected in any of the observed springs where the monitoring duration was one month with a sampling frequency of 1-6 hours. Therefore, this could be because the tracer was absorbed by soil or sediments, and thus this zone cannot be characterized as "highly" vulnerable.
Even though the karst aquifer is the most widespread, there are also numerous springs with lower discharge rates (0.1-5.0 l/s) developed in fissured hard-rock formations (harzburgite, gabbrodiabase, sandstones). The general characteristics of the discharge points for the Tara National Park are shown in Table 1.

Creation of groundwater vulnerability maps
For groundwater vulnerability mapping, the following set of maps has been prepared as a basis for vulnerability assessment: geological map, hydrogeological map, terrain elevation model, terrain slope map, vegetation map, pedological map, a map of existing karst features, precipitation map, and a few others. These maps were used for assessing the parameters needed for applying the DRASTIC, EPIK, PI, and COP methods. Vulnerability classes for the entire study area were obtained by combining the selected parameters and calculating the vulnerability index for each selected method.
As a result of applying the four selected vulnerability methods, groundwater vulnerability maps were created and are presented in Fig. 2. These maps show similarities in some parts of the Tara national park, but also specific differences (see also Table 2).
The intrinsic vulnerability of Tara National Park evaluated by the DRASTIC method shows a predominance of low, medium, and high vulnerability types. The very high vulnerability category occupies only 0.5% of the study area and refers to the zones where the groundwater level is close to the surface, near watercourses. High vulnerability (21 %) is related to karst terrains in the northern and western parts of the park. Most of the area (44 %) is characterized by the medium vulnerability type (the karst plateau around Mitrovac in the central part of Tara National Park). A significant part of the study area (35 %) is recognized for its low vulnerability and is mainly associated with areas chara cterized by fissured aquifers or non-aquifers. According to the results of the application of the EPIK method, a significant part of the study area is classified as a very low vulnerability type (39 %), and represented mainly by nonkarstic rocks (schists, peridotites, and harzburgite). The low vulnerability class (12 %) refers to areas with low permeable rocks or sediments covering carbonate rocks. Medium and high vulnerability types are characteristic for the karstified Triassic limestones, while the very high vulnerability type is present in the immediate areas of the ponors (4 %).
Medium groundwater vulnerability is the predominant class (61 %) on the PI vulnerability map, occupying a large area of both karstic and non-karstic terrains. The low vulnerability class (11 %) is mainly present in the area with non-karstic rocks, but also in karstic terrains such as the Mitrovac plateau, outside the ponor zone, covered with low permeable sediments (clays) with a deep groundwater table. The high vulnerability type is present in the Northwest of the study area, while the very high vulnerability class is limited (2.7%) and refers to ponors, ponor zones, and areas with a large extent of dolines near Mitrovac.
The results of the vulnerability map obtained by the COP method show significant differences compared to other vulnerability maps. Nearly 28 % of the study area is characterized by a  very high groundwater vulnerability class (e.g., ponor catchments and outcrops of the carbonate rocks). Most of the area is characterized by low vulnerability (40 %), while small parts of karst areas outside ponor catchments belong to medium to high vulnerability classes (6 % and 7 %, respectively).

Hazard map creation
The analysis of available data and conducted field research in the area of Tara National Park distinguished potential contamination sources. These hazards are divided into three groups: the first hazard group is classified as contaminators that are the result of infrastructure development (wastewater, fuels, and traffic); the second group consists of hazards related to industrial activity (mining operations and industrial facilities); the third group is related to livestock and agricultural activities. General conditions and safety measures were evaluated for all the identified groundwater pressures. Polygonal hazards are mostly related to villages (urbanization without sewer systems) and intense agricultural activities in the northeastern study area. Local and regional roads represent the main linear type of hazard in the area. The point type of hazards refers to cemeteries, gas stations, quarries, old mines, wood processing facilities etc. The highest point hazard indices are associated with existing lime pits and abandoned mines.
Based on the Hazard index values, a hazard map was presented in accordance with the principles of potentially harmful effects. Although the hazard index is classified into 5 classes, all hazards in the study area are classified into three classes: Very low, low, and moderate degrees of harmfulness (Tab. 3, Fig. 3).

Creation of a risk to groundwater contamination map
The risk of groundwater contamination was calculated for all locations where hazards were detected. Based on the groundwater vulnerability level and hazard index, different risks were calculated for locations where contaminants could be discharged.
The risk maps produced indicate that groundwater vulnerability maps can significantly influence the calculation of the potential contamination risk (Fig. 4). The difference in risk classes can be clearly seen for linear types of hazards that are mainly related to existing roads. The Risk of contamination from these hazards is generally low when the DRASTIC method is used. When the EPIK method is used as a basis, several risk classes (low, medium, and high) are distinguished for karstic areas. On the other hand, the risk of contamination for the linear type of hazard is predominantly low for the non-karstic areas because the EPIK method was developed exclusively for karst areas. In contrast to the previous methods, the COP and PI methods are intended to be applicable for all types of aquifers with particular reference to karst aquifers, and therefore several risk classes are distinguished for both karstic and non-karstic terrains. It is essential to note the presence of the very high vulnerability, which appeared as a result of the COP method application, and which further indicated the presence of the high risk for a linear type of hazards.
The polygon type of hazard is present mainly in the eastern parts of the study area (settlements without sewer systems on nonkarst terrain). Risk maps based on applying the four intrinsic vulnerability methods showed primarily low and medium risk levels, while high risks were obtained for only a few settlements (based on EPIK, PI, and COP vulnerability methods). Agricultural areas are located only in the northeast part of the Tara National Park (on the Drina alluvial sediments) and the risk for these polygons is estimated to be of low and moderate level (DRASTIC, PI, and COP method). By applying the EPIK intrinsic vulnerability method, a very low risk of groundwater contamination was obtained.
The risk of groundwater contamination related to the point type of hazard is generally medium, low, and very low for all the risk maps obtained. Medium risk exists for several limekilns located in the East (according to DRASTIC, PI, and COP map) and a few abandoned mines situated in the southeast of the study area (according to DRASTIC and PI). According to all the vulnerabili ty maps, the same risk category is determined for the existing gas station near Mitrovac. The high risk level of contamination is present in the touristic urbanized area at the Mitrovac (vulnerability based on EPIK, PI, and COP methods) and the wood processing factory in Sedaljka village (according to the COP method). Point type hazards with a very high risk of groundwater contamination were not identified.

Comparative analysis of the results obtained
Generally speaking, all four vulnerability maps are significantly influenced by the geological conditions and the distribution of karstic terrains. This impact is particularly pronounced with the EPIK vulnerability map, where non-karstic terrains are of very low vulnerability. This method does not consider the thickness and structure of the vadose zone and, therefore, should not be used for areas outside the karst. Even for karstic areas, soil protection is not considered in detail. Consequently, low vulnerabili ty zones are not present on maps obtained by the application of other methods.
The precipitation regime is directly analyzed only by the COP method, while the DRASTIC and PI analyses consider the infiltration rate. Most of the Tara National Park represents a flat plateau with similar ground elevations and so significant spatial changes in precipitation rate do not exist. Infiltration conditions play an important role for this area since a part of the karst system is recharged by infiltration of surface water through ponor systems. In such circumstances, the EPIK, PI and COP methods are envisioned since groundwater vulnerability is particularly assessed in the ponor catchments area. However, there are significant differences in groundwater vulnerability even in ponor catchments. All three mentioned methods produce very high vulnerability status for the immediate area around a ponor or sinking stream. But the very high vulnerability is also obtained for large parts of ponor catchments when the EPIK and COP me thods are applied. This is particularly pronounced when the COP method is applied, which causes a significant presence of a very high vulnerability in the entire area (28 %). This vulnerability class is not expected in the entire ponor catchments because of small terrain slopes, lush vegetation, and small sinking streams that occasionally occur. Therefore, the PI method produces a vulnerability map that is most reliable for these areas.
Finally, the soil cover could be of crucial importance in vulnerability assessment. Karst aquifers are usually characterized by high or very high vulnerability, but the presence of less permeable soil could significantly increase the protection role of the unsaturated zone. This is the case with the Tara National Park, where carbonate rocks in the plateau area are covered with clay sediments. The DRASTIC and EPIK methods consider the soil characteristics but with small weighting coefficients, and therefore the impact on the vulnerability index is small. On the other hand, the PI and COP methods calculate the soil protection in more detail, and as a result, low vulnerability is present even in the area of karst terrains.
All of the 424.138 data for each variable were included in comparative analysis. Ordinal descriptive statistics show that even though those methods are called "parametric", the analysis of the distribution (e.g., coefficient of variations, kurtosis, etc.) clearly indicates that all (ordinal) parameters required a non-parametric statistical approach. Significant differences between the variables, the vulnerability maps can be seen on the Box-plot diagram (Fig. 5. -left). On the other hand, there are considerably fewer differences between the Risk maps ( Fig. 5 -right), but these differences can also be important for the final selection of the most appropriate single map for the treated area.
According to the nature of the obtained data, Spearman correlations were applied as a non-parametric alternative to Pearson's correlation. The correlation results of the vulnerability methods show high correlation coefficients between the DRAS-TIC, COP, and EPIK methods. The correlation of the results obtained by the PI method is substantially lower when compared to other methods, particularly to the EPIK method (Tab. 4).
The results show a strong influence of the Hazard Index (Tab. 5) on the vulnerability Risk map obtained by the application of the PI method (0.78), as well as by the DRASTIC method (0.72), while its influence on risk maps produced by COP and EPIK methods is more balanced (0.51 to 0.63). The lowest correlation coefficients of vulnerability and risk map were obtained by applying the DRASTIC and PI methods, at 0.33 and 0.45 respectively, while for the COP and EPIK methods, this coefficient is 0.61 for both methods.
The DRASTIC method does not sufficiently consider the specifics of karstic terrains, while on the other hand, the EPIK method has been developed specifically for karst terrains. However, the results indicate that none of the above are suitable for assessing the groundwater contamination risk in this area because of the different types of aquifers present. The impact of the Hazard Index on the risk intensity index is more pronounced when the PI method is used for intrinsic vulnerability assessment compared to the risk map based on the COP vulnerability index. This correlation does not imply which method is better for vulnerability assessment but points out that hazards have a more significant impact on risk assessment, and therefore, special attention should be paid when sources of contamination are evaluated.

The results of hydrochemical analysis
Sampling points for cations show some differences in the Ca/Mg ratio mainly related to the dolomitic component and Na+K as a result of marl or clay components in rock composition, while anions are less dispersed (Fig. 6). Electrical conductivity for all sampling points ranges from 220 to 500 μS/cm, pH values range from 6.65 to 8.37 (for subthermal springs S-3), and temperature for most observation springs corresponds to average annual air temperature (8-10 o C). There are slightly elevated concentrations for some components originating from clay constituents in/overlying karstified rocks or harzburgites at several springs. Overall results of the chemical analysis show that there is no evidence of chemical contamination in the study area. Therefore, the selection of the intrinsic vulnerability map according to the concept of worst-scenario cannot be accepted as it is too restrictive. Instead, applying the vulnerability method with balanced results should be a more realistic option. Furthermore, the occasional presence of bacteria at the largest karst spring S-1 is a consequence of the reaction to the flood phenomenon and the confirmed connection between the ponor and the spring.

Results of the Risk intensity maps for groundwater source protection
Different risk categories were analyzed for the catchment areas of three springs (Sedaljka, Jarevac, and Perućac). These springs were selected as characteristic examples of various outflows of different aquifers, different catchment sizes, and the presence of hazards.
Spring Sedaljka (S-4, Fig. 1.) is a karstic spring near Sedaljka village in the North-West of Tara National Park. It is a small spring with a capacity of 0.15 l/s used for the local water supply. The catchment of the spring covers an area of 0.5 km² with several hazards, mostly related to the existence of houses without a sewerage system and the presence of a cemetery. Prepared risk maps ( fig. 7) show mostly medium and low risk classes except for the one made by the application of the COP method where a high contamination risk has also been obtained. The remaining three methods (EPIK, PI, and COP) produced good risk results bearing in mind that the spring water quality is good, although water sampling was performed after high precipitation events when lower groundwater quality was expected.
The Jarevac spring (S-18, Fig. 1.) drains the harzburgite rock formation situated in the eastern part of the Tara National Park. The spring is untapped and is not being used for a water supply. The spring's catchment covers an area of 1.1 km². The hazards in this area are mostly related to urbanization without a sewerage system and the existence of local roads. There are also several point hazards, mainly related to tourist centres and one gas station. Good water quality without detected contaminants implies that the groundwater is well protected from the surface. The risk identified by different vulnerability methods is generally low class, so all four maps are representative when it comes to the results of spring water analysis (Fig. 8).
The Perućac spring is the largest karst spring in the area located in the northern part of the Tara National Park (S-1, Fig. 1). A part of the spring is tapped and used for the water supply of a local settlement. The catchment of the spring covers an area of about 70 km², mainly composed of karstic rocks. The spring catchment is slightly populated, and detected hazards are asphalt and unpaved roads, resting houses at the South of the catchment (Konjska river), and touristic centre at Mitrovac located in the middle of the catchment area. The water quality at the spring is controlled. Apart from the occasional presence of pathogenic bacteria, which is very common for large karst springs, the water quality is excellent without the presence of any contaminants. The restrictions related to the National Park are of great importance for the absence of significant hazards, but the main reason for good water quality is that the recharge zone of the spring is the karst plateau which is covered with dense vegetation and thick soil cover. In these circumstances, the vulnerability of karst groundwater does not have to be high or very high, and these types of results were obtained by the implementation of the PI method ( Fig. 9). Other methods resulted in higher risk classes that do not correspond to the actual situation of the preserved water quality.
Overall, the results show a good correlation of the risk maps produced with the spring water quality. The absence of chemical contaminants implies that the risk of contamination should be low to medium, and these risk classes were mainly obtained when the PI vulnerability map was applied. The usage of DRAS-TIC and EPIK vulnerability maps also results in low to medium risk for most of the area. But taking into account that the DRAS-TIC method doesn't consider specific karst features and that the EPIK method gives poor results outside karst areas, it is clear that PI vulnerability and risk map is most reliable for the Tara National Park.   Fig. 4).

CONCLUSION
The presented case study shows that it is necessary to apply several different vulnerability assessment methods in order to be able to select the most suitable one. As a result of the comparison of risk maps by weighting hazard and vulnerability indexes conducted for each applied method, several guidelines can be emphasized: -Risk maps strongly depend on vulnerability maps, and therefore, the appropriate vulnerability method should be chosen very carefully. -The most adaptive vulnerability methods designed to assess any type of aquifer (PI and COP) are more likely to be used as a basis for risk mapping. But these methods require a large data set and detailed field observations. For wider areas, their implementation could be complex. -Karstic terrains require specially designed methods for assessing karst groundwater vulnerability (EPIK, PI, COP). Special attention should also be paid when assessing nonkarst areas. -Statistical procedures can be used as a tool for selecting favorable methods, along with the hydrogeological assess-ment of the treated area. Therefore, the application of several vulnerability methods is recommended to realize risk intensity ranges and not underestimate potential groundwater contamination. -Groundwater quality and the presence of contaminants could be an essential add-on in evaluating risk maps. Water quality should be analyzed at the groundwater source level and compared to the existing contaminants in the catchment of a specific spring. The research conducted in the Tara National Park pointed out the differences in results obtained by applying different vulnerability assessment methods. The presence of large karst aquifers and fractured aquifers formed in magmatic and metamorphic rocks imposes the usage of vulnerability methods designed for different types of aquifers, including the COP and PI methods. However, when ponors are present in the area, the COP method results in large areas with very high groundwater vulnerability, and these vulnerability classes cause the further calculation of high contamination risk. This map does not comply with the groundwater quality in the National park area. Therefore, it is evident that the PI method produced optimal results in this area. It should be mentioned that good water quality is also related to the presence of only small-potential sources of contamination.
The resulting map of groundwater contamination risk produced by applying the PI method has significant importance for the Tara National Park. Intense tourism development in this area potentially leads to the risk of incrementally increasing groundwater contamination. The resulting map represents an essential document for spatial planning and further development of this area and a tool for planning the groundwater monitoring program.