Karstic aquifer vulnerability assessment methods and results at a test site ( Apulia , southern Italy )

Abstract. Karstic aquifers are well known for their vulnerability to groundwater contamination. This is due to characteristics such as thin soils and point recharge in dolines, shafts, and swallow holes. In karstic areas, groundwater is often the only freshwater source. This is the case of the Apulia region (south-eastern Italy), where a large and deep carbonate aquifer, affected by karstic and fracturing phenomena, is located. Several methods (GOD, DRASTIC, SINTACS, EPIK, PI, and COP) for the assessment of the intrinsic vulnerability (Iv) were selected and applied to an Apulian test site, for which a complete data set was set up. The intrinsic vulnerability maps, produced using a GIS approach, show vulnerability from low to very high. The maximum vulnerability is always due to karstic features. A comparison approach of the maps is proposed. The advantages and disadvantages of each method are discussed. In general terms, three groups can be distinguished. The GOD method is useful for mapping large areas with high vulnerability contrasts. DRASTIC and SINTACS are "any-type aquifer" methods that have some limitations in applications to karstic aquifers, especially in the case of DRASTIC. EPIK, PI, and COP, which were designed to be applied to carbonate or karstic aquifers, supply affordable results, highly coherent with karstic and hydrogeological features, and reliable procedures, especially in the case of PI and COP. The latter appears simpler to apply and more flexible in considering the role of climatic parameters. If Iv of each method is considered, the highest variability is observed in cells in the neighbourhood of karstic features. In these spatial domains, additional efforts to define more reliable and global methods are required.


Introduction
Karstic aquifers and environments are highly vulnerable to contamination and to anthropogenic modifications.The vulnerability of karstic aquifers to contamination is due to particular characteristics such as thin soils and point recharge in dolines, shafts, and swallow holes (COST, 2003).The residence time of karstic groundwater is generally shorter, and contamination tends to be faster and simpler, than for nonkarstic groundwater (Kac ¸aroglu, 1999).Specific monitoring criteria should be used for karstic groundwater; sampling in karsts should be more frequent, especially in the wake of storms, rainy periods, or snowmelts, rather than being carried out on a regular monthly or seasonal basis (Vrba, 1988;Polemio, 2005;Polemio et al., 2009).
Anthropogenic modifications to karst aquifers are generally caused by population increases and the associated demands for land, as is the case of karstic areas of the Apulia region in south-eastern Italy (Fig. 1) (Ford and William, 2007;Polemio, 2009).The unplanned enlargement of urban and rural areas is the main example of this.In addition to the effects of pollution from agricultural activities, farming improvements can provoke other negative impacts.Fires and the clearing of forests, associated with the advent of machinery, have been widely employed to provide land suitable for farming, and have also reduced grazing areas.Quarrying activities have disrupted the karstic landscape, destroying the epikarst zone and often causing the destruction of caves.Other modifications can be caused by the mismanagement of water resources or flood risks; the use of sinkholes to accommodate the discharge of flood waters can trigger the formation of dolines, sudden subsidence, cave collapse, and groundwater contamination (Parise and Pascali, 2003;Polemio and Limoni, 2006).All these changes cause severe degradations of the typical landscape and of the karstic environment, the effects of which are often soil erosion, bare karst on slopes, rocky desertification, and changes of the infiltration rates.Groundwater vulnerability generally increases as the protective topsoil is removed and, in quarry areas, as water depth is reduced.The effects on groundwater availability are complex and generally negative.
The socioeconomic development of a region depends also on the availability of good quality water.For this reason, in recent decades there has been a worldwide increase in the demand for water.In southern Europe, water consumption has climbed from 7.1 km 3 /year in 1900 to values 15 times higher in 1995 (Shiklomanov, 1999).Further increases are expected in future years, as the increase can be assessed about 17 times at 2002 on the basis of web databases (World Bank).In many European countries, 50% of the drinking water supply comes from karstic groundwater, and in many areas it is the only available source (COST, 2003), as in the case of the Apulia region (Fig. 1).
Around the world, runoff yield is up to ten times greater than groundwater flow (UNESCO, 2004) except in karstic areas, where the latter prevails.This is the case in Murgia and Salento (Fig. 1), where the groundwater discharge is more than the double the surface runoff, notwithstanding the current overexploitation by wells (Polemio and Limoni, 2006).The modifications to the quality and quantity of current coastal karstic groundwater discharge can have severe effects on the hydrological and ecological equilibrium of the sea and of wetlands located near the coast (UNESCO, 2004).
For all these reasons, it is important to preserve groundwater quality, and particularly to prevent contamination of karstic aquifers.Complex but cost effective procedures that are also simple enough to be used by water and land managers and stakeholders, without specific hydrogeological knowledge, need to be defined.Regional mapping of intrinsic groundwater vulnerability is a tool suitable for these purposes, and specifically for groundwater resources management and protection zoning (Civita, 1994;COST, 2003).
The assessment of groundwater vulnerability quantifies the sensitivity of an aquifer system to groundwater degradation due to human activities carried out at the ground surface.Groundwater vulnerability can not be directly measured in situ; it is a complex function of hydrogeological characteristics of the natural environment that acts as a protection system from contamination (Gogu and Dassargues, 2000) (Table 1).
The assessment of groundwater vulnerability should account for the peculiarities and strong heterogeneity of karst systems, including point or diffuse recharge, rapid flow through high permeability conduits and fractures, and slow flow in low permeability volumes (Doerfliger et al., 1999).General research on vulnerability assessment began in the 1960s (Margat, 1968).In recent years it has focused on the optimal adaptation or new definition of methods for karst aquifers (COST, 2003;Ducci, 2007).A relevant contribution was the definition of the so-called European approach (COST, 2003).The European approach to vulnerability, hazard and risk mapping, based on an origin-pathway-target schematisation, was defined to be general, flexible and nonprescriptive.It is a basis that can be adapted to define methods which are appropriate for specific European karst areas, from Alps to lowlands, and from Mediterranean to continental areas.This paper offers a new contribution to these efforts, comparing the characteristics, procedures, and affordability of several methods for intrinsic vulnerability assessment (Civita, 1994), including GOD (Foster, 1987), DRAS-TIC (Aller et al., 1987), SINTACS (Civita and De Maio, 1997), EPIK (Doerfliger and Zwahlen, 1997), PI (Goldscheider et al., 2000), and COP, based on the European approach (COST, 2003).We focused on methods that have been widely applied and adapted for use in any kind of aquifer (GOD, DRASTIC, and SINTACS), and that have been defined for karstic aquifers (EPIK, PI, and COP).
The study was performed in a test area, located in the karstic Murgia aquifer (Fig. 1), hydrogeologically characterised within the framework of the project WATER-MAP, which was supported by European funds through an INTER-REG III B initiative (Polemio et al., 2007;Corbelli et al., 2008).

Methods of assessment and comparison
Parametric methods of different types were selected (Gogu and Dassargues, 2000).GOD, PI, and COP can almost be considered rating system methods in which ratings, which are linked to fixed ranges of each selected parameter, are summed or multiplied to define the vulnerability index.DRASTIC, SINTACS, and EPIK are parameter weighting and rating methods, as they introduce parameter weights that express the contribution of each parameter to the vulnerability.All these methods assess a vulnerability index Iv or a protection factor π. The range of Iv generally changes for each method, but in any case Iv increases as the vulnerability increases, while π decreases as the vulnerability increases.
The variability of Iv is classified using different attributes (as an example low, moderate, high, and very high) on the basis of a defined sub-range.All methods were defined or can be easily adapted to be applied in a GIS environment.
The GOD method (Foster, 1987) considers three parameters: the groundwater occurrence, the lithology of the layers overlying the aquifer, and the depth to groundwater (distinguishing unconfined or confined conditions).The rating of each parameter is defined from 0 to 1.The overlying lithology contributes to the vulnerability index only in the case of unconfined aquifers (in other words, is equal to 1 for other types of aquifers).At each point, Iv is obtained by multiplying the ratings of the three parameters; Iv ranges from 0 to 1.
The US Environmental Protection Agency developed the DRASTIC method (Aller et al., 1987).This method considers seven parameters: depth to water, net recharge, aquifer media, soil media, topography, impact of the vadose zone, and hydraulic conductivity.Each parameter is rated in the range 1 to 10.At each point, parameters with a variable continuous value (such as depth to water and net recharge) assume a ranking value on the basis of plotted relationships.The relationships consider different types of aquifers and soil media in which the rating choices are based on qualitative attributes described by the authors but assigned by the operator.
Parameter weight factors are defined to balance and enhance the contribution of each parameter to the whole vulnerability.The final Iv value is a weighted sum of the seven parameters.DRASTIC provides two weight classifications: one for normal conditions, used in this article, for which Iv ranges from 23 to 230, and the other one for conditions with intense agricultural activity (Gogu and Dassargues, 2000).The former, called the pesticide DRASTIC index, represents a specific vulnerability assessment approach.In a specific area, only one set of weights should be used.
The SINTACS method was defined in an attempt to improve the DRASTIC method, and uses almost the same parameters (Civita, 1994;Civita and De Maio, 1997).SIN-TACS was defined by Italian hydrogeologists on the basis of numerous test sites located in different hydrogeological and climatic conditions.The main purpose was to define a method versatile enough to be homogeneously applied to each Italian environment (Civita and De Maio, 1997).Some test sites were also located in Apulian karstic areas (Polemio and Ricchetti, 2001;Cotecchia et al., 2002).The method peculiarity is the use of five set of weights or strings (normal impact, relevant impact, drainage, karst, fissured aquifer) to be used simultaneously in large areas, and also if dominated by different prevalent conditions as in the case of zones that are deeply modified by anthropogenic activities, as in urban areas or due to intensive use of chemicals in agriculture, and of different types of aquifers such as porous, fissured, or karstic aquifers.Weights are a very effective tool used to adapt the model to different scenarios, increasing the importance of some parameter and minimizing others (Cucchi et al., 2008).In each cell in which the area is discretised, the set of weights is selected on the basis of conditions that mainly contribute to the local vulnerability.Iv ranges from 26 to 260 and it is classified in six classes.Some authors proposed a modification of the SINTACS score criteria in the case of holokarst areas; nevertheless, up to now this modification has only been tested in Alpine environment (Cucchi et al., 2008).
The EPIK method was defined in Switzerland to be applied only to the vulnerability assessment of karst aquifers (Doerfliger and Zwahlen, 1997;Doerfliger et al., 1999).Four main parameters are considered and mapped: epikarst (E), protective cover (P ), infiltration conditions (I ), and karst network development (K).The E parameter considers the effects in terms of water storage (during rainfall or snow melt) and of the concentration of flow toward vertical conduits; it is assessed on the basis of geomorphological maps considering three sub-parameters.The P parameter describes the protective function of the layers between the ground surface and the groundwater table, mainly soil, subsoil, non-karst rock, and unsaturated karst rock.The I parameter is assessed by distinguishing concentrated infiltration areas and areas in which diffuse infiltration prevails, where the slope and land use are the key sub-factors.K represents the degree of karst network development in the aquifer.The contribution of each parameter to π is assessed by multiplying each parameter by a weighting factor and adding the four contributions.π ranges from 9 to 34 and it is ranked in four classes.The classes can be labelled on the basis of the logical low protection-high vulnerability criterion: a protection factor of 34 indicates the highest protection and also a very low or minimum vulnerability.The vulnerability index is not explicitly defined, but we propose to define Iv=34−π.Practical experiences highlight E and I are the prevalent parameters in the EPIK Iv assessment (COST, 2003).
The PI method was defined in a project founded and realised by German Institutions considering the experience of EPIK method (Goldscheider et al., 2000).The P factor considers the protective function of the layers between the ground surface and groundwater.The I factor describes the infiltration conditions.The P factor ranges from 1 to 5 and can be calculated using a logical schema based on tabled values of considered local conditions, such as topsoil, recharge, subsoil, lithology, and fracturing (Hölting et al., 1995).The I factor ranges from 0 to 1 and describes the infiltration conditions, particularly the degree to which the protective cover is bypassed as a result of lateral surface and subsurface flow in the catchment of swallow holes and sinking streams.The protection factor π is the product of P and I , and ranges from 0 to 5. It is subdivided into five classes of protection/vulnerability.The vulnerability index is not defined, but we propose to define Iv=5−π.The method is also used as basis for hazard and risk mapping (COST, 2003;Mimi and Assi, 2009).The PI method inspired the European approach which considers four parameters: overlying layers (O), concentration of flow (C), precipitation regime (P ), and K (COST, 2003).The parameters O and C are respectively similar to parameters P and I of the PI method.The COP method represents one possible way in which the European approach may be applied (Daly et al., 2002;COST, 2003).The COP method considers the special hydrogeological properties of karsts, and can be applied by considering variable climatic conditions and different types of carbonate aquifers, both diffuse and conduit flow systems (Vías et al., 2006).The COP acronym is justified by the use of three parameters of the European approach: C, O, and P .The conceptual basis of this method is to assess the role of the natural protection of groundwater determined by the properties of the overlying soils and the unsaturated zone (the O factor), of protection weakening due to diffuse or concentrated infiltration processes (the C factor), and of the variable climatic conditions (the P factor).The P factor is innovative as it takes into account not only the mean rainfall but also factors that influence the rate of infiltration, i.e., the frequency, temporal distribution, duration, and intensity of extreme rainfall events.The author's final vulnerability index Ia is determined by multiplying the score of each parameter.Ia is a protection index (if Ia increases, the vulnerability decreases).Ia ranges from 0 to 15 and is grouped into five vulnerability classes: 0-0.5, 0.5-1, 1-2, 2-4, and 4-15 (Daly et al., 2002;COST, 2003).We propose the use of Ib values, defined from 0 to 100, on the basis of five classes, by stretching Ia values of each class with a linear transformation for Ia<1, and with a logarithmic transformation for Ia>1.
In order to compare the results obtained with the different methods, it was necessary to define a standard classification, in which the assessed vulnerability can range from the minimum to the maximum predisposition of groundwater to suffer contamination originating from surface or shallow anthropogenic activities.In, the normalised vulnerability index, was calculated by the normalisation of Iv for each model.The normalisation was realised for each method by attributing In=0 to the minimum Iv score and In=100 to the highest Iv score.The normalization in the I range was performed following a linear law for DRASTIC, SINTACS, EPIK, and PI.In the case of the COP method, In was calculated as In=100−Ib.
In was analysed for the entire set of methods using three steps.The first two steps considered the results of each method.The first step was based on the use of simple statistical tools, such as the calculation of minimum, mean, maximum, standard deviation, range, and the spatial occurrence of the In values.The second step was the qualitative analysis, which utilised In maps and the maps of the main hydrogeological parameters.In this case a visual approach, which can be defined as expert-eye-scanning, was used.The main purpose was to recognise the role of the considered hydrogeological parameters on the spatial variability of In.The third step was based on correlation analysis between pairs of methods, and allowed the quantification of the correlation between the In assessments in each cell with the considered methods.

The study area
The Apulian karstic areas are made up of Mesozoic rocks of the Apulian foreland (Fig. 1).The Apulian platform emerged at the end of the Cretaceous and became part of the foreland of the southern Apennine chain.It was composed of three structural domains: Gargano, Murgia, which includes the study area, and Salento.These structural highs were transformed into islands by subsidence during the middle Pliocene.Transgression led to the widespread deposition of Tertiary-Quaternary formations partially covering the carbonatic platform of Murgia and Salento with thin strata of sand, conglomerates, calcarenites, limestones, and clays.From the middle Pleistocene onwards, the whole region began to be uplifted (Cotecchia et al., 2005).
Apulia hosts large coastal karstic aquifers that form the main regional water source, which is very often affected by degradation in quality due to seawater intrusion (Polemio et al., 2009).The Apulian karstic aquifers are highly permeable due to fracturing and dissolution well below the current sea level, where groundwater flow mainly happens.
The Murgia (maximum height 680 m a.s.l.) is a large asymmetric horst affected by two neotectonic fault systems (NW-SE and NE-SW).Because of these faults, the morphological structure slopes down towards the Adriatic Sea and towards the adjoining regions, by means of a succession of step-shaped ledges bounded by small fault throws.The karst environment consists of platform Cretaceous limestone and dolostone covered with thin layers of Pliocene-Quaternary rocks and soils.The carbonate rock is bedded, jointed, and subject to karst phenomena (Polemio, 2005).The test site (Fig. 2), as large portions of Murgia, is characterised by developed karst landforms that form a network due to the chemical dissolution of limestone (Grassi, 1983).As a result, an underground network of cavities, caves, and conduits has developed, some of which are very large (Regione Puglia a ).The landscape is typically low-relief karst, characterised by many dolines and by very flat-bottom valleys (locally called lame) filled with thin alluvial deposits; mainly residual clayey deposits from karst processes (the so-called terre rosse) and detritus.Outside the test site in the southwestern Murgia sector, typical deep canyon valleys, called gravine, can be observed (Parise and Pascali, 2003).Valleys and water lines are the remnants of the original hydrographic network, which is discontinuous at present.The high Murgia plateau (higher than 300 m a.s.l.) is dominated by clear karstic features and by the existence of discontinuous and thin topsoil layers (less than 50 cm, mainly residual clayey soils and limestone pebbles).Moving from the inland to the Adriatic coast and decreasing in altitude, the shallow karstic features are less evident due to the progressively more continuous and often thicker topsoil cover (Provincia di Bari; Regione Puglia, 1999).

Test site and data
The selected test site (Fig. 2) is 78.2 km 2 wide, and is a portion of the Murgia aquifer.All kinds of geological, hydrogeological, and climatic factors were investigated.The main sources of data are briefly described.
Geological, hydrogeological, and chemical-physical groundwater data were collected from 250 boreholes.The information was extracted from a vast array of sources (from public sources or institutions and secondly from private companies) and dates back to various periods of time.Hydrogeological data, well pumping tests, and piezometric head measurements were checked and upgraded with author's surveys.The quality of the historical data was carefully evaluated considering the location, diameter, depth and screen position of wells, date, type of parameter, and method of measurement (Jousma et al., 2006;Polemio et al., 2009).After the validation process was completed, low quality data were deleted from the data set.
All the data collected for this study were converted into digital format to be implemented in a GIS for the vulnerability assessment.A relational geo-database was designed to permit the simultaneous analysis of all types of data.
The considered parameters yielded some additional data that helped in characterising the aquifer, such as the depth to groundwater, the thickness of the unsaturated zone and of the saturated aquifer, the hydraulic conductivity, and the type of groundwater flow.In order to produce maps of all the different variables measured in the boreholes and useful in assessing vulnerability, the well data were interpolated using a geostatistical approach.A semi-variogram for each variable was made from point data to define the interpolation function.Several raster maps, with 10 m resolution, were then computed for the different variables using kriging techniques or the inverse square distance method (Deutsch and Juornel, 1992) and software (Surfer, Arcview, Arcmap).
The topographical data were mainly obtained from official 1:25 000-1:50 000 hard copies of maps and from digital data provided by the Cartographic Office (Regione Puglia b ).A 10 m grid was used for the DEM calculation, using 5 m spaced contour lines and bench marks.Each gridded calculation (as an example, the distance from a sinking stream or the catchments of karstic features) was realised using a 10 m resolution.Geological information and karstic features (mainly dolines, swallow holes, caves, and sinking streams) were defined using official 1:100 000 geological maps, topographical maps, aerial photo analysis, published speleological and morphological data, and author's field investigation (Regione Puglia a,b ).Land use and cover were defined considering data defined on the basis of Corine classification (Regione Puglia b ).Soil types and thicknesses were obtained from pedological maps and published data derived from field investigations (Provincia di Bari, Regione Puglia 1999).The actual and net yearly rainfall, the number of rainy days, the rainfall rate, and other hydrological calculations were performed by interpolating data from 11 gauge stations, using inverse squared distance interpolation techniques.The used rainfall and temperature monthly data are from 1921 to 2005.The wet years, defined as years in which the precipitation is 15% above annual mean value (COST, 2003), were identified for each rain gauge; using the same interpolation algorithm, the rainfall map of the wet years was obtained (Table 2).
The elevation in the test site ranges from 90 to 330 m a.s.l., and the maximum distance from the sea (Adriatic Sea) is 15 km.The climate is typically Mediterranean; the mean annual rainfall and temperature range from 600 to 657 mm and from 14.8 to 16.2 • C, respectively.The mean annual rainfall of the wet years ranges from 782 to 856 mm.The mean annual value of real evapotranspiration ranges from 494 to 513 mm, and the net rainfall ranges from 105 to 156 mm.The slope is generally low (in 96% of area is less than 6 • ) except where morphological scarps are observed, where the slope is greater than 12 • .A loamy clayey soil less than 1 m thick covers the limestone bedrock; the soil cover is generally thin or negligible (<0.50 m) with the exception of a narrow area located in the north-east (1% of total area).Thick layers of clay soil are present along the drainage network and in the dolines.The mean spatial value of annual infiltration can be assessed equal to 105 mm on the basis of SINTACS method (Table 2).
A wide and thick aquifer is located in the carbonate rocks of Murgia.The aquifer is divided into more permeable strata due to the variable distribution of fractured and karstified strata, confined between less permeable levels of various extents and thicknesses, mainly due to tectonic events that fractured the carbonate rocks in a discontinuous manner and to the variation of base level of flow.The groundwater flow is generally confined except along a narrow coastal strip; faults govern the major preferential flow paths and seawater intrusion.
The hydraulic conductivity ranges about from 10 −4 to 10 −3 m/s, and the depth to groundwater ranges from 158 to 444 m below the ground surface.

Results
Figure 3 shows the final maps of In.The In GOD map is not shown, as In was equal to 20 everywhere.This homogeneous value was due to two circumstances: the aquifer is not phreatic (was considered semi-phreatic in the GOD application) and so the effect of unsaturated strata was considered null in terms of protection, and the depth to water was everywhere in the lower vulnerability class (>100 m).
The In values can be classified in five classes of vulnerability: very low (0.0-19.9), low (20.0-39.9),moderate (40.0-59.9),high (60.0-79.9),and very high (80.0-100.0).This classification scheme was applied to each In output; the results in terms of percentage spatial occurrence are summarised in Table 3.
The GOD assessment indicated the lowest In values with the unique exception of the COP method, which indicated In<20 over 2.5% of the test site and 18.1 as minimum In value (Table 3).This method seemed to underestimate the vulnerability and to provide a low sensitivity to the spatial variation of key hydrogeological parameters.These characteristics were also observed by other authors (Civita and De Regibus, 1995;Corniello et al., 1997;Gogu and Dassargues, 2000).GOD should be mainly used to map vulnerability assessment for very large areas, with small scales (1 100 000-200 000) or with high vulnerability contrasts, especially if the data set is poor.For these reasons, GOD was not considered further.
The DRASTIC values were in the range of 23-63 (Fig. 3; Table 3); in practical terms, they defined two classes, low and mainly moderate vulnerability.A peculiarity of the In DRASTIC map, as well as the SINTACS map, was the evidence of faults in terms of a slight increase of vulnerability with respect to neighbouring pixels (Fig. 3).This was due to the assumption that faults could create an increase of fracturing within a strip, defined with a 20 m-wide buffer in this case.The buffer was used only in the case of parameter "A" in both methods, which considers the role of the aquifer media or material.The DRASTIC assessment showed sensitivity to parameters such as some karstic features (caves, dolines, and swallow holes).
The SINTACS assessment was in the range of 39-79.Excluding GOD, SINTACS shows the lowest variability or standard deviation, and is about equal to the DRASTIC values.Two vulnerability classes were again observed, but these were of one level higher than the DRASTIC values, from moderate to the more prevalent high vulnerability.At a preliminary comparison with DRASTIC, SINTACS In values and statistics seemed to be similar, but were higher by about 10-20.A detailed eye-scanning of the SINTACS map highlighted that it was also somewhat sensitive to the drainage network, the effect of which was absent in the DRASTIC assessment.This was due to the effect of the drainage set of weights (Civita and De Maio, 1997).The disadvantage was that the traditional drainage set of weights generally reduced In, with respect to neighbouring pixels in these hydrogeological conditions, in which the karstic set of weight was used, and determined a In decrease which did not seem to be realistic everywhere.The SINTACS output appeared complex and was influenced by relevant parameters as depth to water, hydraulic conductivity, tectonic faults, drainage network, lithology, soil cover, and urban areas.
The EPIK results spanned a large In range, from 24 to 100, and had a low mean value equal to 39.8.High to very high vulnerability was observed in 2.5% of the test site, mostly in narrow areas surrounding karstic features and sinking streams.The largest spatial occurrence was observed in the low vulnerability class, which included 89.2% of the test site.The Low values observed in the north-eastern sector were mainly due to the effect of discharging outside the karst area (Fig. 2).The moderate class assessment covered 8.3% of the area, but the pattern did not appear to be clearly linked to evident or objective features or parameters.The whole assessment, as in the case of the GOD method, appeared underestimated.
The PI method defined a very large range, 20-100, and the minimum of the means of each method.These values were almost equal to the respective EPIK statistics; the distribution of the class occurrence was quite similar to EPIK as well, but there were some relevant differences.The classes from moderate to very high vulnerability showed slightly larger spatial occurrence due to the reduction of the low class occurrence.The main difference appeared after the expert-eye-scanning: the In variability over the whole observed range was determined in large areas surrounding karstic features, and was clearly due to the effect of the role of catchments of sinking streams or to areas discharging into the karstic subsurface as it was reasonable.Low variability was observed only in the north-eastern sector, the main part of which discharges outside the karst region.
The COP method defined the largest range, from 18 to 100, and so was the only method that spans all the vulnerability classes.The mean In value, equal to 74, was the maximum of all the methods.The class occurrence was statistically bimodal (as was EPIK, while the remaining cases were unimodal) due to the secondary peak (2.5%) of very low vulnerability, and to the absolute maximum (84%) of high vulnerability.The very low vulnerability class was due to the effect of alluvial deposits outcropping in the north-western sector, while the low vulnerability class (0.4%, the minimum spatial occurrence) was due to the effect of the thicker soils in the north-eastern sector.In any case, the COP output was mostly from high to very high vulnerability, classes that cover 96.5% of the area.
The third step of the analysis highlighted that the highest correlation was between the DRASTIC and SINTACS outputs, which appeared similar, though the In values appear to be translated or shifted to greater values from the former method to the latter (Table 4).The correlation analysis confirmed the good agreement between the EPIK and PI results, which was already observed qualitatively.It also seemed relevant that the correlation analysis determined two groups of methods, SINTACS and DRASTIC in one group and PI and EPIK in the other, that show very low correlations.The COP method was the only method whose output was almost correlated with any other method output (Table 4).

Conclusions
The results of this study highlighted that the karstic test site shows aquifer vulnerability that ranges from low to very high.The highest contribution to vulnerability was due to karstic features such as dolines, swallow holes, caves, and sinking streams.
If the In value of each method except the GOD method was considered for each cell, the highest variability between methods was defined in the neighbourhood of karstic features.In these spatial domains, additional efforts to define more reliable and global methods are required.The minimum vulnerability in each cell was generally generated by the PI, EPIK, and DRASTIC methods, while the maximum was due mainly to the COP method and secondly to the SIN-TACS method.
The GOD method should be used mainly to map vulnerability assessment for very wide areas, with small scales or with high vulnerability contrasts.
The DRASTIC and SINTACS methods are "any-aquifer" methods that consider the greatest number of parameters.For this reason, they produce results that can be used at very high resolution or at large scales (if the dataset is robust enough, as in this case).They also define a fine variation of In over a narrower range.On the contrary, the methods defined for karstic or carbonate aquifers, EPIK, PI, and COP, seem to be black-white or low-high methods; in other words, they are able to discriminate in the karstic environment the potentially most dangerous locations for pollution sources.
The SINTACS assessment seems to be more accurate and flexible than the DRASTIC method.It takes into account the role of different prevalent conditions, such as drainage and high anthropogenic modifications, but some of the weights probably require further optimisation.
Between the karstic aquifer methods, PI and COP show more realistic results on the basis of the expert-eye-scanning approach and reliable procedures.The latter appears simpler to apply and more flexible in considering the role of climatic parameters; the final maps seem to reflect the author's knowledge of the vulnerability characteristics of the specific aquifer and test site.
Further efforts and research should be performed.The role of preferential flow paths and of epikarst or shallow subsoil need to be better schematised with tools that permit the af-fordable quantification of their contributions to the vulnerability.Additional efforts are necessary to improve the vulnerability assessment methods, moving from intrinsic vulnerability to a complete vulnerability assessment that includes hazard mapping, the characterisation of groundwater quality degradation due to existing potential sources of pollution, and the critical analysis of the assessment procedure.

Fig. 3 .
Fig. 3. Maps of the vulnerability normalised index In of the test site.

Table 1 .
Methods and used parameters for the groundwater vulnerability assessment.

Table 2 .
Statistics of spatial values of morphological, climate and hydrogeological parameters in the test site.(1)Infiltrationassessed on the basis of the SINTACS method; (2) wet years, defined as years in which the precipitation is 15% above annual mean value, are considered in the COP method(COST, 2003).

Table 3 .
Statistical values of the normalised vulnerability index (In) and spatial percentage occurrence in the test site of In classes.

Table 4 .
Correlation coefficients between In maps.