Karst groundwater vulnerability assessment: application of an integrative index-based approach to main catchments of middle Valseriana springs (Northern Italy)

The study aim was to assess the intrinsic vulnerability of the middle Valseriana (northern Italy) through a methodology that could fully consider the peculiar characteristics of the Pre-Alps region. After applying classical methods for karst environments (COP and EPIK), the integrative COPA+K method was developed for the Nossana (dominant drainage system) and Ponte del Costone (dispersive circulation system) spring catchments. The COPA+K approach exploits the potential of the COP method with the addition of two factors: the influence of the main discontinuities on the water system according to their distance from the spring (A factor), and the development of the karst network (K factor). The COP and COPA+K maps were validated based on the characteristic correlation, for the study area, between the values of δ18O and the elevation of recharge areas. Compared to COP, the COPA+K results appeared more stringent in the identification of areas characterized by high vulnerability, which passed from 35.6% to 23.6% of the entire study area. COPA+K also made possible to better differentiate the behavior of the two catchments. In terms of high susceptible areas, their percentage difference increased by 12.3%, correctly emphasizing the greater susceptibility of the Nossana drainage system in comparison with the Ponte del Costone dispersive one. For the COPA+K map, in the Nossana catchment, the δ18O values indicated a good agreement between the mean elevation of high vulnerable and water recharge areas (± 106 m). For Ponte del Costone, the isotopic data indicated a possible mixing with the waters of the Serio River.


Introduction
Karst environments and karst aquifers are widely known for their vulnerability to pollution. The causes of this high predisposition to pollution can be found in the nature and characteristics of these environments: They are mainly characterized by significant lateral extension, shallow soils, and highly effective recharge points such as sinkholes, sinking streams, and polje (Zwahlen et al. 2003). The carbonate deposits, the main actors of karst phenomena, cover 15.2% of the Earth (21.8% of the European territory) (Goldscheider et al. 2020). The strategic importance of these environments is also given by the large amount of drinking water they provide (Stevanović 2019); about 20-25% of the world's population depends partly or entirely on them (Ford and Williams 2013). Their contribution in terms of drinking water supply is of great value, but considering the predisposition of these environments to human pollution, more attention should be paid to the preservation of this resource and the maintenance of the balance of the entire karst system (Bakalowicz 2005;Zhang et al. 2017;Padilla and Vesper 2018). Selfpurification processes, which are normally occurring in large surface aquifers on the plains, are less efficient due to fast infiltration processes, low filtration capacity, and high flow rates (Kačaroğlu 1999). The transfer time of the water from the surface to a spring depends on the intensity of the karstification: Through some channels the water flows quickly, while in others it can be retained for a longer period in the rock matrix. Therefore, contaminants can take from a few hours to some days to reach a spring, or they can be stored underground for weeks, months, and, in some cases, years (Vesper et al. 2000;Green et al. 2006).
The vulnerability of aquifers identifies the propensity of groundwater systems to be impacted by anthropogenic contaminants from the surface (or subsurface). This concept implies that the natural environment and geologic context may yield some degree of protection to groundwater with respect to human impacts (Vrba and Zaporozec 1994). If not better specified, the term vulnerability can be interpreted with three different meanings (Civita 1987(Civita , 1994. The first is intrinsic vulnerability of an aquifer to pollution, which represents the specific susceptibility of groundwater systems to handle and diffuse a fluid or waterborne pollutant such as to produce groundwater quality deterioration, in space and time. Intrinsic vulnerability is a function only of the geological, hydrological, and hydrogeological characteristics of the area. The second is integrated vulnerability, which refers to the intrinsic vulnerability in association with certain potential sources of pollution existing in the considered area. The third is specific vulnerability, which is instead referred to a particular contaminant (or group of contaminants) and is a function of the properties of the contaminant that characterize its behavior underground.
The concept of intrinsic vulnerability, proposed in the European Action COST 620 (Zwahlen et al. 2003), is based on the conceptual origin-flow-target model for aquifers in carbonate rocks (Daly et al. 2002). Using this model, it is possible to distinguish between the protection of resources and the protection of water supply (Marín et al. 2012). Considering the protection of resources, it is assumed that the element to be protected is the water stored in the aquifer below the piezometric surface. The flow of the contaminant from the point of origin is considered vertical, passing through the soil and the unsaturated zone. In the case of the protection of supply, the goal of protection is the point of extraction of the groundwater, in a well or a spring. Conceptually, it is assumed that the contaminant is transported from the groundwater table and from there to a spring or a well. This displacement has a vertical component from the origin to the piezometric surface and a horizontal component.
The assessment maps of the vulnerability of aquifers to pollution represent a basic tool for the protection of the water resource to be implemented over the entire recharge area of an aquifer. Vulnerability maps were first presented by Albinet and Margat (1970). During the 1970s and 1980s, numerous authors addressed the issue, attempting to deal with the concerning pollution phenomena in highly urbanized areas and suggest specific groundwater protection interventions (Montolio et al. 1979;Foster and Hirata 1988). These vulnerability maps are designed to provide the targets for the implementation of prevention and protection measures as well as management indications for the location of new industrial plants, landfills, intake, or input units. They constitute a system of protection of the groundwater resource at basin scale that is particularly useful because of the direct and immediate connection of the spring with the entire recharge area in karst environments. Over the years, different types of techniques for the assessment of vulnerability have been defined: statistical, process-based simulation, and index-based techniques (Beretta 1992;Focazio et al. 2002). The latter are the most used and preferred because they are not encumbered by computational complexities and data shortage (Kumar et al. 2015). Furthermore, the combination of these techniques with the GIS software management potential has undoubtedly facilitated and expanded the application horizons of these methodologies (Doerfliger et al. 1999;Goldscheider 2005;Nguyet and Goldscheider 2006;Goldscheider 2007, 2009;Mimi and Assi 2009;Polemio et al. 2009;Mimi et al. 2012;Doummar et al. 2012a;Bardai et al. 2015;Kazakis et al. 2015Kazakis et al. , 2018Iván and Mádl-Szőnyi 2017;Nanou and Zagana 2018).
The main goal of this study was to define an integrative methodology that could represent, in the most correct way possible, the conditions of intrinsic vulnerability of the middle Valseriana (northern Italy, Province of Bergamo). In the past, several authors have tried to adapt or create new index-based techniques to assess the vulnerability of the study area as correctly as possible (i.e., Leyland 2009;Jenifer and Jha 2018;Tziritis et al. 2020). The purpose of this study was to define an approach that does not require a lot of data and expensive investigations available in order to make it applicable in mountain contexts, where many hydrogeological and hydrochemical data are not available as in the plain areas. Moreover, it was necessary to formulate a method that respected the constraints imposed by the drinking water exploitation of the springs, which excluded the possibility to perform tracer tests to avoid compromising water quality. For map validation, tracer tests are the most commonly applied methodology too (Goldscheider et al. 2008). Alternative possibilities include: hydraulics and spring hydrographs analyses, with which karst storage and infiltration processes can be defined; analytical and numerical modeling, which are theoretically the best solution but at the same time complex and demanding in terms of data requirements; and natural tracers, which are particularly useful for the identification of recharge zones (Zwahlen et al. 2003). Given an ongoing water sampling campaign for environmental tracer analysis, these data represent a validation opportunity for the study area.
Two of the most important aquifers in the Middle Valseriana valley, the Nossana and Ponte del Costone spring catchments, were studied characterizing the natural conditions, surface, and subsoil, regardless of the actual presence of sources of pollution. The data used were related only to the geological and hydrogeological characteristics of the basin, including climate data, vegetation cover, land use, topography, lithology, structural data, and karst features.

Study area and data
Nossana and Ponte del Costone springs are located in Valseriana valley (Fig. 1), in the Province of Bergamo (Northern Italy). These spring systems are set within an area characterized by a vast outcrop of carbonate platform deposits oriented toward southeast (Calcare di Esino and Dolomia Principale formations), which constitute the aquifers; the Argilliti di Riva di Solto and the Carnian formations (base layer, characterized respectively by successions of argillites and micritic limestone-marl alternations) represent aquicludes, while the Calcare di Zu formation (roof layer) shows an aquitard behavior. In this context dominated by Triassic formations, the permeability of the aquifer is mainly related to the interconnection of fracture systems and karst phenomena within the Calcare di Esino and the Dolomia Principale formations. The catchments of Nossana and Ponte del Costone springs show substantial differences. The Nossana basin is a dominant drainage system spring (Jadoul et al. 1985;Vigna and Banzato 2015) characterized by a high permeability, which is linked to a considerable karstification process. Conversely, within the Ponte del Costone system the circulation is mostly dispersive and the permeability of the aquifer is mainly related to the families of interconnected discontinuities (Spada and Bertuletti 1988;GeoTer 2000;Ecogeo S.r.l. 2007). The Ponte del Costone karst network has a reduced development compared to that of Nossana.
Both springs are exploited by UniAcque S.p.A., the company that ensures the drinking water service to the city of Bergamo and the communities of the valley. The two basins investigated show very different behavior in terms of discharge ( Fig. 2) due to the different areal extension (area of about 80 km 2 for Nossana and about 10 km 2 for Ponte del Costone catchment) and the different degree of karstification of the rocks. The Nossana spring has a discharge that varies between 0.5 m 3 s −1 and 18 m 3 s −1 , while the Ponte del Costone system, composed of 13 aligned springs (Fig. 1), has a cumulative discharge between 0.15 m 3 s −1 and 0.45 m 3 s −1 . The richness from the point of view of water resources is also given by the abundant precipitation that characterizes the whole area of Valseriana: the average precipitation is close to 2000 mm/year and can reach peaks of about 3000 mm/year (Ceriani et al. 2000).
In the study area, the karst phenomena affect the most superficial portions of the rock masses. The Federazione Speleologica Lombarda (FSLo 2020) reports the presence of several smaller caves, with a generally descending development for tens of meter. Moreover, following geomorphological field works, swallow holes, sinkholes, and infiltrations within the detrital deposits have been identified: This study reports these elements because they would represent, other than good recharge zones, also contaminants catchers and therefore areas of high vulnerability. However, it should be noted that in the study area there are few and limited centers of danger.
The geological, geomorphological, and elevation data of the area were obtained from the Geoportale della Regione Lombardia (Regione Lombardia 2020a, b), supported by field evidence and previous studies (GeoTer 2000;Ecogeo S.r.l. 2007). For the development of the karst network and the evaluation of the propensity to karstification of the area, the work of FSLo (2011) was relied on. The meteorological data are available on the open-data section of the Environmental Regional Agency (ARPA Lombardia 2020a, b, Table 1), while the discharge data were provided by Uni-Acque S.p.A.
A water sampling campaign was conducted in the study area between May 2018 and July 2019. Water samples from Nossana (ID-1 in Fig. 1) and some Ponte del Costone springs (Costone 1, Merlo, and Sottoplatea, which are ID-2, ID-3 and ID-4 in Fig. 1, respectively) were collected monthly for the whole sampling period. In addition, water samples related to the Serio River were collected once per year at two locations, upstream (ID-5) and downstream (ID-6) of the confluence of the Riso Stream into the Serio River ( Fig. 1). Measurements of δ 18 O and δ 2 H were taken on the water samples using an isotope ratio mass spectrometer (IRMS) referring to the SMOW (Standard Mean Ocean Water) standard, represented by the mean isotopic content of ocean waters (Craig 1961).

Methodologies applied
Initially, two of the best-known methods to assess the vulnerability of karst environments were applied, EPIK and COP. Following, a new integrative approach was developed that could best simulate the territorial conditions of the study area. The data, available in vector or raster format, were processed in GIS environment using various analysis tools (buffering, clipping, etc.), conversion tools to transform vector data into raster, and Raster calculator functions for mathematical calculations and overlapping of various layers. The area considered in the data processing measures 150 km 2 (× min: 562,000 m, × max: 572,000 m, y min: 5,075,000 m, y max: 5,090,000 m; WGS 84 UTM zone 32 N) and the raster cells have dimensions of 5 m × 5 m.

EPIK: method and application
EPIK is an assessment method specifically developed for karst rock aquifers, regardless of current land use and  socioeconomic considerations. The method leads to the mapping of the vulnerability of the recharge area of a karst spring considering four parameters (Doerfliger 1996;Doerfliger et al. 1999): epikarst development (E), protective cover conditions (P), infiltration conditions (I), and karst network development (K). The sum of these four weighted attributes produces a protection factor (F), which is assigned to each cell of the investigated catchment: where α, β, γ, δ represent the weights of the respective E, P, I, and K factors and they are equal to 3, 1, 3, and 2, respectively. Finally, the values of the reduction of the protection factor (F) are divided into four classes: very high, high, moderate, and low. In this study, the determination of the E factor was based essentially on the detection of the presence of karstification on the surface; this parameter is therefore representative of the water collection process. E factor has been divided into three classes (E1, E2, and E3), outlined according to the distribution of sinkholes, caves with the presence of water, and sinking streams in the study area. Recharge areas located within 50 m around a sinkhole are generally defined as highly vulnerable (Dunne 2004); therefore, they (and equally for sinking streams, swallowholes, and caves) were assigned the E1 class, with the corresponding score of 1 point. Buffer areas of 500-m around sinkholes, caves, swallow-holes, and sinking streams (Doummar et al. 2012b) and 50-m around watercourses were assigned a value E2 (3 points), while the rest of the catchment was assigned the value E3 (4 points).
The surface protection is the thickness of the soil or other no-karst geological formations above the main aquifer. The thicker the protective cover, the more significant is the protection, which means the higher the value of the P factor. This parameter is divided into four classes, which vary from 1 to 4 points. The evaluation of the P factor was performed using the soil thickness data extrapolated from the soil map at the 1:250.000 scale (ERSAF 2010) and the data from the 077 CLUSONE geological map (Jadoul et al. 2012). The soil map divides the present soils into five classes based on depth (very deep, deep, shallow, thin, and very thin), also providing a value in cm of depth. In the study area, the prevailing lithology are made of carbonate, and thus soils mainly belong to group A. Very thin soils were attributed class P1 (1 point), thin and shallow soils class P2 (2 points), and deep and very deep soils class P3 (3 points). Soils with depths greater than 20 cm placed in correspondence with low permeability and no-karst formations, such as Argilliti di Riva di Solto and Carnian formations (GeoTer 2000), were assigned to classes defined for group B (2 or 3 points). In the study area, the P4 class was not assigned.
Following Dounmar et al. (2012b), for the mapping of the I factor, the recharge areas of sinkholes, caves, swallow-holes, and sinking rivers (buffer of 50 m) as well as the buffer areas of 50 m around the watershed were assigned class I1 (1 point). Only the waterflows directly connected to the sinkholes, both perennial and temporary, were considered. For the remaining area, a slope map was derived, thanks to which it has been possible to assign, within the recharge area of the aforementioned watercourses, class I2 (2 points) to the slope areas > 25% and class I3 (3 points) to the slope areas < 25%, being present pastures and fields. Class I3 was also assigned to areas with slopes larger than 25% and located outside of the spring catchments. As for the P factor, no areas attributable to I factor class 4 were recognized.
For the evaluation of the K factor, the study of FSLo (2011) was used. Based on the karstification propensity index (of values 1, 2, 2 + , and 3) of the map of karstification degree of the Lombardy Pre-Alpine belt and the caves surveyed, classes K1 (e.g., Calcare di Esino and Calcare di Zu), K2 (e.g., Dolomia Principale), and K3 (e.g., Argilliti di Riva di Solto and Carnian low permeable formations) were assigned to the different lithology (1, 2, and 3 points, respectively).

COP method
COP is the acronym of the three main factors used to assess the groundwater vulnerability: concentration of water flow (C), overlapping layers considering both soil and rock substrate (O), and precipitation (P). This method evaluates the ability of the overlying layers, in particular, soil and unsaturated zone (O factor), to mitigate the contaminant. In addition, since karst aquifers are characterized by widespread and concentrated infiltration, the C factor defines the importance of infiltration processes, while the P factor emphasizes the role of climatic conditions (Vías et al. 2006). Likewise to EPIK, the method consists of analyzing the three main parameters and producing the corresponding maps. The obtained values are finally multiplied to define the COP index, a protection index to which different vulnerability degrees are attributed.
The final map is rated according to the five classes of vulnerability: very high, high, moderate, low, and very low.
In the COP method, the basin is divided into two main zones: The first zone (Zone 1) includes the karst elements recharge area, and the second zone (Zone 2) consists of the rest of the area.
The C factor for Zone 1 consists of multiplying three main sub-factors: where dh is the distance from the karst elements, sv refers to vegetation and slopes, and ds is the distance from the secondary streams that feed the karst network. The factor C for Zone 2 is instead equal to: where sv has the same meaning as in Eq. 3, while sf refers to the development of the karst network and the permeability of the more superficial layers.
The sv parameter is composed of two maps: one for the slopes and another for the density of vegetation. The slope values, in percentage, were reclassified into four categories (≤ 8%, 8-31%, 31-76%, and > 76%), to which the relative scores were assigned. Based on the data of the DUSAF project (ERSAF 2018; Regional land-use map), it was possible to discriminate three classes: areas with vegetation from areas without vegetation or with low vegetation development. Through a matrix combination among the cells of the two raster, the map related to the sv parameter was generated. For the distance from swallow-holes, sinkholes, and sinking streams (dh parameter) a series of buffer areas located at certain distances from such karst elements were considered. Thus, a series of concentric buffer areas around each karst element was generated. For the calculation of the distance from the karst drainage channels (ds parameter), two buffer areas around the watercourses directly connected to karst features (the same chosen in the EPIK method) were created: one with a distance of 10 m from the watershed and one of 100 m.
The O factor represents the overlying layers, namely the ground cover (O S ) and the bedrock (O L ). The COP method considers through the O factor the hydrogeological characteristics of the rock, including the actual porosity and hydraulic conductivity, the degree of fracturing (Iy), the thickness of the layers (m), and the degree of confinement of the aquifer (cn). The O S parameter evaluates the soil, considering its framework and thickness. The O S parameter map for the study area was developed based on the 1:250,000 scale soil map, assigning the scores suggested in the COP procedure to combinations created according to soil type and depth. The O L parameter is instead representative of the unsaturated zone. It is obtained according to the following equation, in which the summation of the various products of ly and m, calculated for each lithology present, is obtained separately, reclassified, and multiplied by the degree of confinement: The representative values of each type of lithology are given to the various types of outcropping rocks in the catchment area; the higher the value, the lower the vulnerability (for example, fractured rocks are assigned a value of 3, while marls are assigned a value of 1000).
The P factor is the sum of two sub-parameters (P I and P Q ) that define the amount and intensity of annual precipitation, respectively. P I considers the intensity of precipitation, in other words, the ratio between the amount of precipitation and the number of rainy days. It varies between 0.2 and 0.6, and it is the result of interpolating annual precipitation average from ten meteorological stations (ARPAL Regional agency) around the study area from 2013 to 2018. P Q instead represents the amount of precipitation. It is scored according to a table. As the annual precipitation at the considered meteorological stations varies between 1650 and 1100 mm/ year, the scores vary spatially between 0.4 and 0.2.

COPA+K integrative method
In order to decide the best applicable method between EPIK and COP, to perform a more extensive and detailed analysis, the validity of the final maps should be verified by testing with artificial tracers (Doummar et al. 2012b). Their application would allow to verify or not the actual vulnerability of the various areas. However, given the strategic importance of the springs for the drinking water supply of Bergamo and surroundings, in accordance with the service company managing the aqueduct network, it was preferred to avoid the use of tracers in the absence of an appropriate tracing plan. Besides, the wide extension of the catchments (altogether about 90 km 2 ) makes the implementation complicated. An isotopic study of the catchments water is underway to identify the main recharge zones so to allow the drafting of an operative tracing plan. Without tracer tests information, based on previous studies conducted in karst areas that compare the two methods (Doummar et al. 2012b(Doummar et al. , 2012aGuastaldi et al. 2014;Hamdan et al. 2016), COP was considered the most reliable. It evaluates the precipitation of the area and has a better construction for the calculation of the final score without assigning indiscriminately specific weights to the different factors.
To further improve the effectiveness of the tool at the local level, it is important to evaluate some aspects not considered in the resource vulnerability map obtained by applying COP. As suggested by Zwahlen et al. (2003), it is possible to extend the target by also assessing the horizontal flow within the system, shifting the focus to the protection of spring drinking water. In this way we obtain a source vulnerability map from the combination of the new factors and O, C, and P.
The new first element to consider is the influence of tectonic features on the water system. Areas characterized by high discontinuity density often correspond to areas with relatively high groundwater pollution (Sarikhani et al. 2014). Several examples can be found in the literature of integration to classical methods or new methods of vulnerability assessment that consider a discontinuity factor (Mendoza and Barmen 2006;Pochon et al. 2008;Awawdeh and Jaradat 2010). The innovative element of this study is the different scores given according to the distance from the spring. This approach arose from the fact that it was possible to represent two specific aspects of the study area. The first is to highlight the rapid water contribution that discontinuities located closer to the spring provide to the aquifer. Karst areas are strongly conditioned by the tectonic history of the rock mass (Eraso 1986). When the hydraulic gradient coincides with the main drainage direction (as it happens in this study area, being oriented southeast toward the springs), the karst conduit system tends to develop linearly (Eraso 1986;Bakalowicz 2005). The structure of a linear, and therefore simple, karst system increases the likelihood that discontinuities and fractures located near the source will favor a rapid water contribution. The second aspect is to distinguish the different behavior of the studied karst systems. The hydrograms in Fig. 2 show how the Nossana spring discharge can increase of 12 m 3 s −1 in a single day as a response to an intense rainfall event (Citrini et al. 2020), while the Ponte del Costone springs record minimal variations that cannot be associated with rainfall (GeoTer 2000; Ecogeo S.r.l. 2007). Thus, the main discontinuities within the two catchments were selected based on the geological surveys summarized in the CLUSONE geological map (1:25.000, Jadoul et al. 2012) and a 200 m buffer was run around them. Based on the area of the spring hydrogeological catchment, two different radii were calculated for the buffer areas originating from the spring, according to the following formula: where n is equal to 2 for the first buffer area (equal to half of the catchment area) and equal to 1 for the second buffer area (equal to the whole catchment area). As an example, for the Nossana spring catchment (≈ 80 km 2 ), a first buffer equal to 3.5 km and a second buffer equal to 5 km were obtained. The main discontinuities, which are those shown in Fig. 1, were attributed a score of 0.2 if within the first buffer area, a score of 0.5 if within the second buffer area, and a score of 0.8 if anywhere else in the catchment. Such a scoring scale allows giving more importance to a structural discontinuity that lies within the first buffer than to one that lies at a further distance from the spring. This new factor is called factor A, namely "Association between the main discontinuities and their distance to spring," and it aims to (6) r buffer = √ catchment area∕n identify areas distributed around the main discontinuities that can considerably influence the infiltration conditions of the catchment area. The second element to be combined with the COP method comes from the European Approach (Daly et al. 2002;Zwahlen et al. 2003): the evaluation of karst network development (K factor). In this way, the vertical path of recharge water (with or without contaminant) is combined with the mainly horizontal path through the saturated substrate to the considered springs. Examples of the addition of the K factor to the COP method was already proposed in the past in two different ways. The first is the so-called Slovenian Approach of Ravbar & Goldscheider (2007). It involves a limited modification of some values of the original COP parameters to make the method better applicable to karst aquifers in Slovenia. These method variations include the addition of a sub-factor tv (temporal variability) to the C index, two new sub-factors rd (rainy days) and se (storm events) for the calculation of P index in place of the original ones, and the consideration of very thin or completely absent protective cover of soil for the O index. The additional K parameter is calculated based on the residence time in the saturated zone under high-flow conditions (based on tracer tests), the development of the karst network and the degree of channel connection, and the contribution of different parts of the aquifer system to the spring. The second approach to include K in COP was presented by Andreo et al. (2009). The authors proposed the COP + K method, which involves the use of COP without modification added to the K factor, calculated as suggested by Ravbar and Goldscheider (2007).
In this study, the second approach of Andreo et al. (2009) was followed, but using the K factor as derived for EPIK, because of the impossibility to perform tests with tracers. The K factor was reclassified using the score proposed by Andreo et al. (2009): class K1 was given the score 0 points, class K2 the score 1 point, and K3 the score 2 points. In this way, in the final sum of the raster of K to that of COP, in the presence of class K1, the vulnerability remains unchanged, not reducing the score in the possible areas of high vulnerability. Instead, in the presence of classes K2 and K3 the vulnerability changes by 1 class or 2 classes, respectively: In the presence of K2 and K3 the protection is higher, since the karst network is less developed and, thus, the vulnerability decreases.
The integrative COPA+K approach (Fig. 3) is therefore created by combining the results of the COP method and the A and K factors, according to the equation:

Creation of intrinsic vulnerability maps
For all methods, each factor map was calculated individually and then the weighted sum or multiplication that characterizes the specific vulnerability assessment methodology was applied. In the EPIK application, the weighted sum of the four factors produces the F protection index map, where each zone is assigned a value between 9 and 34 Fig. 3 Integrative COPA+K* method. The main steps of the COP method (Vías et al. 2006) are highlighted in purple, the integration of the karst network development according to Andreo et al. (2009) in green, and the integrations of this study in orange points. The values were then divided into four classes: low (> 25 with P4 and I3, or I4), moderate (> 25 points), high (20-25 points), and very high (9-19 points). For the COP approach, the multiplication of the three maps by each factor produced the resulting COP map for the catchment areas. The final map was then reclassified according to the five classes of vulnerability suggested by the COP method (0-0.5 very high, 0.5-1 high, 1-2 moderate, 2-4 low, and 4-15 very low). In the same way as COP, the classification of the COPA+K method divides the study area into five classes from very high to very low. The same classification is maintained but the score range of the lowest class of vulnerability is extended to 17 for the consideration of the development of the karst network. To perform the comparison between EPIK and COP, it was possible to overcome the different number of classes of the two methods by reclassifying the COP map into four classes. This was done by combining the high and moderate vulnerability classes of the COP method in a single high vulnerability class.

Validation process
Early results from the isotopic analyses conducted between 2018 and 2019 were exploited for the vulnerability map validation. Indeed, environmental isotopes represent a tracking method for infiltrating precipitation or surface water (Zwahlen et al. 2003). In particular, the δ 18 O value measured in groundwater samples is an effective tool for estimating recharge areas, as its variations are strongly correlated with elevation and geomorphological features (Tang et al. 1998). In this study, for map validation purposes, it was verified that the mean elevation of the high vulnerability zones (Very High and High classes) obtained by applying the COP and COPA+K methods match the mean spring recharge elevations as derived from δ 18 O data. To perform the validation, four methodological phases were followed.
First, the mean elevation of the recharge areas of the Nossana and Ponte del Costone springs were estimated from the reference δ 18 O and δ 2 H values of Valseriana precipitation and the local isotopic line (Table 2) provided by Longinelli and Selmo (2003). Based on these data, a linear regression was created to derive elevation as a function of δ 18 O values.
Second, mean annual precipitation data from ten ARPA meteorological stations in the area surrounding Valseriana valley (Table 1) were regionalized to obtain a gridded (50 m × 50 m) precipitation distribution. For the regionalization, two methods were applied, consisting in a linear regression (LR) with altitude and a following interpolation of the residuals either by thin plate spline (Duchon 1977, TPS) or inverse distance weighting (Shepard 1968, IDW). Third, the elevation values of the high vulnerability areas and related precipitation amounts were extrapolated from the DTM and the previously calculated precipitation maps (LR + TPS and LR + IDW), respectively. Then, through a precipitation-weighted average, the average elevation of the high vulnerability areas was obtained for both the Nossana and Ponte del Costone catchments. Such values were finally compared with the results of the isotopic correlation.

Vulnerability maps
The final intrinsic vulnerability maps for the EPIK and COP methods are presented in Fig. 4. Focusing on the EPIK approach, the areas characterized by the highest vulnerability are mainly the recharge areas of the sinkholes and watercourses directly connected to them, while the areas near the springs belong to the moderate vulnerability class (Fig. 4, a). Regarding the COP technique (Fig. 4, b), the C factor strongly influences the result, mainly because the areas located at a distance of 50 m from the sinkholes and watershed were attributed 0 points. This significantly influences the multiplication with the vegetation and slope parameters and the other two factors (O and P), reducing their score to zero.
Even though both methods are applicable to karst aquifers and show several similarities, the results are not completely comparable. High vulnerability areas are more or less delineated in the same positions in EPIK and COP. In the map obtained through COP, high vulnerability areas mainly include buffer zones around sinkholes, caves, and swallowholes. In the map obtained through EPIK, these areas are more irregular and restricted since the method itself is more selective. The high and moderate vulnerability areas of the COP method generally correspond to the high vulnerability areas of EPIK, while the low and very low vulnerability areas of COP are usually equivalent to the moderate and low vulnerability areas of EPIK (Fig. 5, a). Similar observations were done by Hamdan et al. (2016). The differences are partly due to the different number of vulnerability classes, which in COP are 5 while in EPIK are 4. Indeed, by combining the high and moderate vulnerability classes of the COP method in a single high vulnerability class, the results of the two methods become more similar (Fig. 5, b). Compared to the previous classification, the reclassified COP has an unchanged class VH, a class H resulting from the combination of classes H and M, a class M corresponding to the previous L, and a class L resulting from the original VL. Assigning a value from 1 to 4 to the classes defined in EPIK and reclassified COP, starting with the minimum values for high vulnerability classes, it is possible to observe the differences between the map generated through COP and the map generated through EPIK (Fig. 5, c). The differences in vulnerability maps are due to the different importance given in the two methods to the vulnerability assessment parameters. In fact, some parameters of EPIK correspond to those of COP, but they are given different weights. Despite this, about 57% of the area considered does not show any difference between the results of the two methods, about 40% differs for one class, and only about 3% of the total area shows a difference greater than two classes.
The fundamental COP factors are those concerning the framework and thickness of the soil, which strongly influence the overlying layers (O factor) and consequently the final vulnerability value, and the thickness of the unsaturated zone (O L ), especially in karst and fractured carbonate rocks. Parameter dh (karst elements, in the calculation of parameter C) is the most relevant, since areas within a distance of 50 m from the sinkhole are given a value equal to 0, which is able to cancel both the factor O and the factor P in the final calculations (Tayer and Velásques 2017;Kovarik et al. 2017). Finally, the P factor plays an important role only in the case of a moderate to low vulnerability, while in cases where the C and O factors are very low, the precipitation factor (P) plays a negligible role.
For the EPIK method, the main parameters that influence the protection factor (F) are the epikarst (E) and the development of the karst network (K), since they are assigned a weight of 3. The protective cover (P) also plays an important role given the weighting coefficient (2), although it has less influence than the factors mentioned above. Infiltration conditions (I), which in other words correspond to slope variations, influence the boundaries of the vulnerability areas. Due to the importance of some parameters, there are uncertainties in both methods due to the quality of the data used for the evaluation of the various factors. In both methods, a regional scale soil map (scale 1:250,000) was used, Fig. 4 Results of the application of (a) EPIK method and (b) COP method consequently with a not high degree of detail; in EPIK, this map was used for the parameter P, which has minor importance in the method. In COP, it was used for the parameter O S (O factor), which has a big influence on the final value. In EPIK, further uncertainty is related to the E factor because literature values (Doerfliger 1996;Doerfliger et al. 1999;Doummar et al. 2012b) were used for the definition of recharge areas around sinkholes and not values specific to the case in question due to lack of data.
Regarding the COPA+K approach, the maps of the five considered factors are shown in Fig. 6. The final map of this integrative method distinguishes five classes of vulnerability (Fig. 7). The consideration of the main discontinuities has produced areas of very high and high vulnerability along the main tectonic lineaments, but without neglecting the epikarst evidence. Focusing on the effects of the K factor, the K2 class is the most widespread and, for this reason, the protection rises in most of the area. The vulnerability classes very high (VH), high (H), moderate (M), low (L), and very low (LW) within the considered area have percentage of 11.2%, 12.4%, 21.0%, 46.4%, and 9.0%, respectively. Considering the most vulnerable classes (VH and H), the COPA+K method allows the identification of more restricted areas than COP. In fact, the values move from 35.6% (COP) to 23.6% (COPA+K) over the whole study area. This is therefore indicative of a better prioritization of the areas on which action needs to be taken to preserve the water system or to implement effective monitoring systems. Furthermore, COPA+K makes it possible to highlight better than COP the different behavior of the two spring systems. Considering only the highest vulnerability classes VH and H, the percentage difference between the Nossana and Ponte del Costone catchments increases from 5.2% to 17.5%, for COP and COPA+K, respectively. These results reflect the evidence of the discharge rates and the initial assumptions, the Nossana spring is more vulnerable than the Ponte del Costone springs due to its geological setting and the dominant drainage behavior (Vigna and Banzato 2015). In addition, all the high vulnerability areas have been identified in areas where groundwater resource prevention measures have already been implemented (spring protection area) or are mountainous areas with very low population density (seasonal silvicultural and breeding zones). This is comforting and ensures greater preservation of aquifer quality.    Given the conservative properties of the stable isotopes in the water (Gat 1996), mixing with the Serio waters can be demonstrated. In fact, the δ 18 O value of a hypothetical recharge water A, originated at an elevation of 1000 m a.s.l., can be derived from the local isotopic correlation (Fig. 8). Such Water A would be consistent with the orography of the Ponte del Costone catchment and the derived very high and high vulnerability areas. Figure 10 shows how the average isotopic characteristics of the Ponte del Costone springs water, in terms of δ 18 O vs δ 2 H, fall between Water A and the Serio River water (Fig. 10). Considering the isotopic characteristics of Water A and the Serio River water (ID-6) as shown in Fig. 10, it is possible to estimate a river water contribution at the Ponte del Costone springs that ranges from about 35% at Merlo (ID-3) to about 75% at Sottoplatea (ID-4). At Sottoplatea, water surfaces at a couple of meters above the Serio River, making the difference between the two springs reasonable.

Maps Validation
The evidence of the possible mixing between Ponte del Costone springs and Serio River does not exclude that the estimate of the average recharge elevation within the catchment (i.e., the estimate of the most vulnerable areas) might be subject to an error. However, its demonstration coupled with the excellent elevation match obtained for the Nossana catchment encourages considering the derived vulnerability map as a reliable product. Furthermore, the comparison of the COP and COPA+K methods (Table 3) shows how the new proposed approach allows to achieve the best validation results. Through an exhaustive isotopic study, it could be possible to ascertain and accurately define the mixing of Ponte del Costone springs water with that of the Serio River. This would help not only to completely validate the COPA+K map for the Ponte del Costone catchment, but also to define the relationships between the springs and the river to prevent any problem, originated outside the catchment, that could threaten the quality of the exploited water.
Many vulnerability assessment methods aim to have the widest possible range of application for different geological environments (i.e., Tziritis et al. 2020), but no vulnerability assessment approach is robust enough to meet the needs of all types of geological settings (Kumar et al. 2015). Each study area has, in fact, its own characteristics. Hence, trying to adapt classical methods to the study area is the best way to achieve a good result and to best represent the real conditions. However, this can lead to one of the most common problems when applying indexbased methods: the subjectivity of the author in assigning values to the different factors involved (Tziritis et al. 2020). Therefore, a balance must be found between the faithful application of the method and the modifications dictated by the study area. Following this modus operandi, an improvement in the classic COP method has been achieved to obtain an integrative and innovative method that can be applied and reproduced throughout the entire karst Pre-Alpine belt.

Conclusions
In this study, the COP, EPIK, and the newly designed COPA+K methods were applied and evaluated for a groundwater vulnerability assessment on the middle Valseriana area. The study area is characterized by two main karst spring systems with different characteristics. The Nossana spring is a drainage system, while the Ponte del Costone springs system is dispersive.
The two classical methods showed similar and comparable results but were not able to discern the different behavior of the two spring systems. For this reason, the idea was born to modify one of the two applied methods to create a new integrative approach with suitable sensitivity. This new method was created from COP since, compared to EPIK, it considers the rainfall contribution in the catchment and does not apply weights in the final calculation.
Two additional factors were proposed for the innovative assessment approach: the A factor, which considers the influence of the main discontinuities concerning their distance from the karst spring, and the K factor, which considers the development of the karst network. The K factor was not determined as in the integrative methods proposed in the literature through tracer tests (Ravbar and Goldscheider 2007;Andreo et al. 2009), but using the one already calculated in EPIK, due to the impossibility to apply these types of tests. Despite this constraint, the estimation of the K factor using literature data proved to be sufficiently accurate. VH-H areas correspond to the extension of the most vulnerable areas according to the COPA+K method The COPA+K map defined five classes of vulnerability (very high, high, moderate, low, and very low) that allowed determining more precise areas to prioritize the implementation of contaminant risk prevention and mitigation measures. Most of the areas classified as highly vulnerable are fortunately already subject to protection measures or are located in mountainous areas and not at the valley floor. These areas are intended for seasonal silvicultural and breeding use and currently do not raise any concern regarding the preservation of water quality. Moreover, the application of the COPA+K method allowed to underline the different responses of the two considered water systems. The percentage difference between the high vulnerability areas (VH and H) of the Nossana and Ponte del Costone catchments is increased by 12.3% compared to the COP method (from 5.2% for COP to 17.5% for COPA+K approach).
The validation of the vulnerability maps carried out based on the relationship between recharge areas elevation and δ 18 O values in spring waters confirmed the reliability of the map obtained with the proposed new method and its higher quality than that derived with COP.
Maps of intrinsic vulnerability such as this one, in association with maps of territorial planning and constant monitoring of water quality, represent a strategic tool for the management and protection of mountain water systems such as Nossana and Ponte del Costone. Finally, the COPA+K method is applicable to the entire karst environment of the Pre-Alpine belt due to its easy applicability (it does not require much data and a considerable computational effort) and because it has been verified on two of the most important spring catchments in Northern Italy in terms of discharge and water supply. This approach could also be useful in karst contexts where it is not possible to use tracer tests to assess the development of the karst network. As this case demonstrates, the availability of environmental isotopic data represents an asset to validate the approach.
Funding Open access funding provided by Università degli Studi di Milano within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.