Challenges and Limitations of Karst Aquifer Vulnerability Mapping Based on the PaPRIKa Method—Application to a Large European Karst Aquifer (Fontaine de Vaucluse, France)

Abstract

Aquifer vulnerability maps can improve groundwater management for sustainable anthropogenic development. The latest update of karst aquifer vulnerability mapping is named: the Protection of Aquifers base on Protection, Rock type, Infiltration and KArstification (PaPRIKa). This multi-criteria assessment method is based on a weighting system whose criteria are selected according to the aquifer under study. In this study, the PaPRIKa method has been applied in the Fontaine de Vaucluse karst aquifer using the novel plugin for Quantum Geographic Information System (QGIS) software. The Fontaine de Vaucluse karst aquifer is the largest European karst hydrosystem with a catchment area that measures approximately 1162 km${}^2$. Four thematic maps were produced according to the criteria of protection, rock type, infiltration, and karst development. The plugin expedites the weighting system test and generates the final vulnerability map. At a large scale the vulnerability map is globally linked with primary geomorphological units and at the local scale is mostly affected by karst features that drive hydrodynamics. In conclusion, the novel QGIS plugin standardizes the application of the PaPRIKa method, saves time and prevents user omissions. The final vulnerability map provides useful contributions that are most relevant to groundwater managers and decision-makers. We highlight the sensibility of the vulnerability map to the weighting system and validation issues of the vulnerability map are raised.

1. Introduction

Independent of the nature of a contaminant and contamination scenario, intrinsic vulnerability is an expression of an aquifer’s geological and hydrogeological properties that define its susceptibility to pollution [1,2]. To determine aquifer vulnerability, most methods use the origin-pathway-target model [3]. Vulnerability is a qualitative indexing of a flow pathway through a geological sequence between a contaminant release and a target [4]. The target might be groundwater (resource protection) or an outlet (spring protection).

Graf [5] reviewed a large panel of vulnerability mapping methods and classified them into three approaches: Physical distributed modelling, index methods, and statistical models. The most common approach used in the groundwater vulnerability assessment of karst aquifers may be index methods. According to this method, primary geological and hydrological characteristics are mapped and clustered from least to most likely to promote vulnerability. Identified aquifer components are combined according to a weighting relationship to compute the vulnerability index. Index methods are based on average long-term aquifer components and are easy to use. But the subjective index weighting system prevents a comparison of vulnerability maps [6].

The European Cooperative in Science and Technology program, action 620 [3] for karst groundwater vulnerability mapping established a European approach using key aquifer infiltration factors [7]. Based on its own groundwater regulations, each country has developed vulnerability mapping methods based on the European approach [8,9,10,11]. To develop a standardized strategy for karst groundwater protection, a French consortium of groundwater decision-makers and researchers developed the PaPRIKa method [4]. In detail, PaPRIKa stands for the Protection of Aquifers incorporating four criteria: P for protection (considering the most protective aspects among parameters related to soil cover, unsaturated zone, and epikarst behaviour), R for rock type, I for infiltration, and Ka for karstification degree [12].

The novel plugin for Quantum Geographic Information System (QGIS) software [13] makes it possible to apply the PaPRIKa method to karst aquifer vulnerability mapping in a standardised way. Developed within an open-source environment, the toolbox provides a clear workflow to establish vulnerability maps. In this study, the QGIS toolbox was applied by conducting vulnerability mapping on the Fontaine de Vaucluse catchment area (France) which is the largest European karst hydrosystem. We highlight both the advantages of the PaPRIKa method and the choices available to hydrogeologists by using the QGIS toolbox, as well as the challenges and the limitations of the method.

2. The Fontaine de Vaucluse Karst Aquifer

2.1. Catchment Area

The Fontaine de Vaucluse spring is the main outlet of a watershed that covers 1162 km${}^2$ in southern France (Figure 1, [14]). The watershed is bounded to the north by the mountain range of Mont-Ventoux (1912 m a.s.l.) to the west, and the Montagne de Lure (1826 m a.s.l.) to the east. It is also bounded to the east by the Durance valley and to the west by the Rhône valley (Figure 1, Plaine du Comtat Venaissin). The karst system is limited to the south by the Apt plain (Figure 1. Forest and bush cover 84% of the surface. Farming occupies 15% of the watershed and 1% is human habitation structures. Hence, soils are mainly lithosol, calcisol, and fersiallitics in pockets or in cracks, and locally very rocky [15]. Soil texture varies from sand to clayey loam.

2.2. Climatological and Hydrological Background

The impluvium of the Fontaine de Vaucluse is under a Mediterranean climate. Winters are mild and humid and summers are hot and dry. The average temperature is 10 °C and the average annual rainfall is 960 mm·y${}^{-1}$. There is no permanent hydrological network (Figure 1). The landscape is marked by numerous dry valleys and dry canyons, which proves that infiltration processes are prevalent over runoff processes. The aquifer is fed by precipitation alone, there is no exchange with adjacent aquifers. The flow of the Fontaine de Vaucluse has been accurately measured by a gauging station since the end of 2003. Daily flows range between 2.8 and 62 m${}^3$·s${}^{-1}$ with an average of 15 m${}^3$·s${}^{-1}$.

2.3. Geological Background

Sedimentary conditions in the Provence area during the Cretaceous resulted in a thick limestone series, about 1500 m [14]. This sedimentary series is composed by a large panel of carbonate types, with marls, cherty limestone, sandy limestone or bioclastic limestone.

Cretaceous limestones crop out on almost the entire catchment area (Figure 2). The limestone is covered by tertiary deposits only at the edge of the catchment area or within large ditch collapses that occurred during the extensional phase that began during the Late Priabonian. Due to lateral limestone facies variation and discontinuous marl facies, the entire thickness may be karstified [14,16].

The catchment area is delimited on the northern side by the mountain chain of Ventoux-Lure. This mountain chain is the southern limb of an anticline (Figure 1). The south side of this mountain chain forms regional highlands (Plateau d’Albion or Albion highlands). Between Mont-Ventoux and the Albion Mountain is the Sault rift with a NNE-SSW orientation. Post Gargasian deposits cover the rift; this material is more impervious than Cretaceous limestones. The Banon rift is located to the east of the catchment area between the Albion highlands and the southern side of the Lure Mountain. The rift is oriented NNE-SSW and is highly fractured. The presence of the two rifts is due to overstretching conditions during the Cretaceous period [17]. Since the Cretaceous, numerous tectonic phases have affected this structure; thus highly variable stresses induced possible fracturing in all directions. NNW-SSE faults are primarily due to compression events, whereas NNE-SSW faulting is mainly due to extensional events [18]. The entire catchment area is affected by both fault families. A few faults have an E-W orientation as a result of an extensional event during the Miocene; these are primarily observed at the outcrop scale. But this is a predominant orientation of cave development, so it seems that this family plays an important role in karst development [17]. Since the upper Albian, the Durancian uplift has eroded post Cretaceous deposits and intensified karstification processes. Since that time compressive and extensional events have also promoted erosion. Bartonian stage in particular are marked by uplift of the Mont-Ventoux mountain chain and formation of canyons and sinkholes. Karstification continues today. As of today, speleologists have explored about 364 cavities (Figure 1).

2.4. Hydrogeological Background

The catchment area has a mean elevation of 880 m a.s.l. The elevation of the aquifer outlet at the Fontaine de Vaucluse spring is 84 m a.s.l. [16]. Thus, mean thickness of the unsaturated zone (UZ) is about 800 m [14]. Due to its thickness and lithology, the UZ plays an important hydrological role as a buffer stock of water [19].

Flows from the unsaturated zone are observed within the Low Noise Underground Laboratory of Rustrel (LSBB, Figure 3). Flows from shallow underground to 500 m below the surface are measured. The observed flows are intermittent or permanent with variable or stable discharges. Based on water chemical quality, Barbel-Perineau et al. [20] classified the organisation of flows with depth within the unsaturated zone.

The flows in the saturated zone were characterized by tracing tests. Six tracing tests were carried out between 1963 and 1974. They provide evidences about the hydraulic connection between the Nesque canyon, the Albion plateau and the area east of the Banon ditch with the Fontaine de Vaucluse spring. The maximum apparent transfer velocities are between 12 and 208 m·h${}^{-1}$, depending on whether the injection occurred in underground runoff under low water conditions or near a major drain under flood conditions [18].

The rainfall-discharge correlation analysis of the Fontaine de Vaucluse hydrosystem shows that the pressure transfer between the surface and the saturated zone is very fast, between 1 and 6 days [21]. Wavelet analysis of the rainfall-discharge relationship shows that flow variations are clearly influenced by monthly and annual rainfall [22]. This analysis also highlights that long-term trends of discharge evolution are mainly due to changes of the measurement method rather than changes in rainfall patterns. Climate change has not significantly affected regional precipitation yet and is therefore not visible in the flows of the Fontaine de Vaucluse spring. However, the very good hydraulic connectivity between the surface and the spring shows that if the precipitation regime changes, it will have an almost immediate effect on the hydrosystem of the Fontaine de Vaucluse.

3. Outline of the PaPRIKa Method

The PaPRIKa method is based on EPIK (E: epikarst, P: protective cover, I: infiltration conditions, K: karst network development) [24] and RISK (R: rock of the aquifer, I: infiltration conditions, S: soil and protective cover, K: karstification degree) methods [25]. It use an analysis of geological and hydrological aquifer properties to map intrinsic vulnerability of the studied karst aquifer. For more details the reader is invited to read the PaPRIKa method guide [4,12,13]. In addition, useful application cases can be found in Marìn et al. [26], Huneau et al. [27], Kavouri et al. [28], Kazakis et al. [29]. The complete workflow of the method also as the QGIS plugin steps are illustrated in Figure 3.

Four thematic criteria are considered to build the related thematic maps: (a) Protection (P) refers to aquifer properties that delay infiltration (soil, epikarst, unsaturated zone), (b) Rock type (R) considers the lithology and the fracturing degree of the saturated zone, (c) Infiltration (I) differentiates diffuse from concentrated infiltration and (d) Karstification degree (Ka) represents the functionality of the karst system, based on karst network organization and development.

Criteria P and R are related to aquifer structure and criteria I and Ka represent hydrosystem functioning. These four criteria are mapped independently and categorized into five classes, from the least to the most likely to promote vulnerability [13]. Weighted criteria are used to compute the vulnerability index (Equation (1)).

$$V=p\ast P+i\ast I+r\ast R+k\ast Ka$$

(1)

where V is the vulnerability index, p, r, i and k are weights of P, R, I and $Ka$, respectively. This weighting relationship is specific to each study site. The exact values of p, r, i, k, are thus defined by the user. However, Dörfliger and Plagnes [4] advise that the sum of i and k should range between 50 and 65% and the sum of p and r between 35 and 50%. The PaPRIKa QGIS toolbox provides an intuitive workflow for producing P, R, I, $Ka$ maps and computing a vulnerability index (Figure 4). It is important to mention that the data layers to be informed within the QGIS plugin are the result of a preliminary data collection and analysis by the user. Each data layer contains information that the user has previously collected and classified regarding the theme of the data layer. After their accurate implementation, the QGIS workflow facilitate calculation of parameters P, I, R and Ka according to the scheme in Figure 4.

4. Results

4.1. Aquifer Structure Criteria (P and R)

The R criterion refers to saturated zone lithology classified from the least to most potentially karstified: Predominantly limestone, marly limestone, and predominantly marl. The expected nature and thickness of saturated zone lithology is established on the basis of geological maps, local geology series, and geological sections [14]. Structural properties such as faults are classified according to their direction and it is assumed that faults directed toward the outlet are less protective than others. The related map is illustrated in Figure 4.

The P criterion at Fontaine de Vaucluse is based on the leakage ability of soils, unsaturated zone (UZ) thickness, and subsurface water reserves (Figure 5). The higher the protection is, the lower the vulnerability index remains. The P criterion related thematic map is illustrated in Figure 4. Due to the Fontaine de Vaucluse structural complexity additional details about the P criterion are given for each parameter in the following parts.

4.1.1. Soil

The soil map of the study area [15] provides information on the proportion of different soil types (rendzine, lithosol, calcosol...) within each morphological unit, such as the south side of Mont-Ventoux. We propose to use the mean leakage index of soils of these morphological units to qualify aquifer protection due to soil as illustrated to

Table 1. Soil protection index according to soil leakage index.

Leakage Index	Descritpion	Equivalent Protection Index
1	Very porous soil with very fast infiltration	4
2	Fast infiltration	3
3	Moderate infiltration	2
4	Infiltration is not total	2
5	Minor infiltration	2
6	Infiltration weak	2
7	Infiltration very weak	1
8	Infiltration close to nil	1
9	Impervious soil	1

.

4.1.2. Epikarst

The presence of 84 small epikarstic springs in the catchment area allows subsurface water transfer to the surface. These subsystems might contain surface pollution or delay water infiltration toward the groundwater table. The catchment area has no perennial hydrographic network, and as a result, all surface water mainly infiltrates. As a result, the catchment area of a spring may be considered to be protective of the Fontaine de Vaucluse aquifer, but depending on flow velocity, the spring may be a point where surface pollution concentrates. Moreover, it is challenging to define the catchment area of each of all the epikarstic springs because their area depends on landscape and subsurface geology. So, we use spring density to delineate epikarst aquifer location at the regional scale. Following the proposed epikarst classification by Kavouri et al. [12], most of the study site has an index of 3, and area with epikarst aquifer has an index of 2.

4.1.3. Unsaturated Zone (UZ)

The criterion called unsaturated zone is the conjugation of two properties of the aquifer structure, the thickness of its unsaturated zone and the tectonic faults. Faults with a buffer zone of 25 m are assumed to be less protective than zones where there are no faults, thus tectonic faults have an index of 4. The UZ thickness is greater than the size established by the PaPRIKa method [4]. The unsaturated zone flows are hierarchic organised depending on depth Barbel-Perineau et al. [20]. As a result, in

Table 2. Unsaturated zone thickness and corresponding index.

Index	UZ Thickness
1	>400 m
2	200 to 400
3	15 to 400
4	<15

, we propose an updated indexing of the unsaturated zone according to the flow ranking.

4.2. Aquifer Funtionning Criteria (I and Ka)

Morphological components such as slope gradient and surface karst features are used to construct the I criterion. The Ka criterion represents the impact of karstification on the hydrological behaviour of the aquifer. Karst features are used for both maps but in different ways.

4.2.1. Karst Features

We assume that karst features observed at the surface are linked to the karst network promote rapid pollutant transfer. Three kinds of surface karst features are considered: Caves, sinkholes and dry valleys.

• Cave descriptions are used to determine their vulnerability index, including maximum cave depth and the horizontal conduit size. Threshold values of cave depth are the same as for the unsaturated zone; the deeper the cave, the higher the index. This index also depends on horizontal cave development depending on a classification: 0 to 20 m, 20 to 100 m and more than 100 m. Moreover, cave density is used as a karstification index, so if numerous caves with small extent are near each other (less than 100 m), their index is higher.
• Analysis of the digital elevation model makes it possible to locate endorheic areas such as sinkholes. Sinkhole impact on surface flows is ambiguous because it can enable rapid infiltration in the case of open holes [30], or delay infiltration and promote evaporation in the case of sinkholes filled with clay [31]. Sinkholes have an index of 3.
• Dry valleys concentrate surface runoff, so their index is set to 2.

Surface karst features, dolines, and caves occupy a surface area of about 7 km${}^2$ of the watershed.

4.2.2. Karst Spring Behavior

The discharge classification proposed by Mangin [32] characterises this system as a karst aquifer with a mature karst network, well-developed conduits, and a large area of drowned karst downstream. According to Mangin’s classification, the Fontaine de Vaucluse aquifer belongs to class 1. So according to the PaPRIKA guide, the entire catchment area has a Ka index of 2.

The Ka map is globally set to an index of 2. This index increases locally due to karst feature density or exceptionally extensive caves. An example of such a cave is the “Trou souffleur” (Blowing Hole), which has been explored by speleologists to a depth of 600 m under the land surface [14]. Both Ka and I thematic maps are illustrated in Figure 6.

4.3. Intrinsinc Vulnerability Map

The intrinsic vulnerability map of the Fontaine de Vaucluse karst aquifer is illustrated in Figure 7. The chosen weighting system places major emphasis on the infiltration factor, with a weight of 50% on the I factor, 20% each for the P factor and R factor, and 10% on the Ka factor. The “very low” class is absent on P, R, Ka, and I maps; therefore the intrinsic vulnerability of the aquifer is always higher than “very low”. The weighting system minimizes the influence of P and R and weights slope highly. The presence of a fault systematically leads to the attribution of a “high” to “very high” vulnerability index, particularly, in the centre of the catchment area where numerous tectonic features, karst features, and caves were identified. Steep slopes are areas of low vulnerability. The northern part of the catchment area is a mountain chain with a high slope gradient. There, the unsaturated thickness is estimated to be greater than elsewhere in the catchment.

The most vulnerable areas are surface karst features (e.g., dolines, sinkholes) and major faults. These karst and tectonic objests are supposed to act as drains that bypass the buffer role of the unsaturated zone. They represent 2% of the watershed’s area. Forty-five percent (45%) of the watershed has a high vulnerability index. These areas correspond mainly to the high plateaus where flat areas and karst shapes promote rapid infiltration. An intermediate vulnerability index is assigned to 48% of the area, corresponding mainly to the steep slopes of the mountains. Finally, the lowest vulnerability index is assigned to 5% of the area which is mainly characterized by very steep slopes where runoff processes are the majority.

5. Discussion

5.1. The Appropriate Weighting System

The QGIS PaPRIKa toolbox facilitates the vulnerability index computation step. In this project we tested six weighting systems and covered a wide range of index changes (

Table 3. The tested weighing system for intrinsic vulnerability mapping of Fontaine de Vaucluse.

Version	p	r	i	k
V0	0.25	0.25	0.25	0.25
V1	0.2	0.2	0.5	0.1
V2	0.2	0.2	0.3	0.3
V3	0.2	0.2	0.4	0.2
V4	0.3	0.1	0.5	0.1
V5	0.3	0.1	0.3	0.3

and Figure 8). The difference between the weighting systems tested never exceeds one point, which means the variation is always between two consecutive vulnerability classes. The reference weighing system $V0$ is an equal weighting of each factor. The comparison of the five other systems to $V0$ shows that induced changes at most 30% of the studied area. Regardless of the system used, two areas always have a vulnerability index higher than $V0$. The first is the airport located in the centre of the area. The second is the north side of the Lure Mountains, located in northeast of the area. This is due to the R factor, these two areas have moderate protection due to their lithology. In fact, the slope classification makes it possible to determine whether the vulnerability map will result in a pessimistic or optimistic assessment of aquifer protection.

5.2. Vulnerability Assesment Validation Issues

The reliability of the vulnerability maps is critically important. Hence, new methods are suggested and modern statistical and simulation techniques are tested so as to validate and improve the accuracy of vulnerability maps [33]. However, before testing new methodological approaches we should answer the following question:

What does reliability of a vulnerability map mean?

Concerning karst aquifers this question has been answered in previous works (e.g., [29]) highlighting the necessity for high frequency and spatial distribution of hydrological and hydrogeological data, to understand the function of the karst system and considering the response of the system under different simulation scenarios. Undeniably, the abundance of such data sets is usually rare for karst hydrosystems and validation methods are used in order to increase the reliability of the vulnerability map. In this study, the validation process includes a trial and error approach in order to choose the best weight combination of the parameters. Each set of weights provided the corresponding vulnerability map and then evaluated according to the empirical knowledge of the experts (in the present case, the authors). The weights of each criterion in index-based methods can be also determined from the expert’s opinion [34]. In this study, the novel plug in allowed the opposite approach due to the fast application of the PaPRIKa method. Nevertheless, the validation of the weights could be obtained using methods such as sensitivity analysis [35,36].

6. Conclusions

6.1. The PaPRIKa Method and the Related QGIS Plugin

The QGIS plugin toolbox that has been developed based on the PaPRIKa method provides a fast and secure method for creating thematic maps and vulnerability maps. The vulnerability map is original because it is based on an individual estimation of input parameters and a weighting system [8,26,27], and the plugin makes it possible to standardize the application of the method. The operator can easily share data and decision rules used for vulnerability mapping and obtain repeatable results. This is the first step toward an inter-comparison of vulnerability maps. The most useful vulnerability map is the one that best describes the hydrogeological system and the operator’s subjectivity. The possibility of easily testing several weighting systems supports the operator’s subjectivity and provides more evidence to managers. Applications of the PaPRIKa highlight the:

• Vulnerability assessment validation tools;
• Sensitivity analysis of criteria.

The QGIS plugin will be updated in the future in order to include the single sensitivity analysis tool and further simplify the vulnerability assessment in karst aquifers. Additionally, an extra tool for the assessment of the unsaturated zone sub-parameter is under development. The aforementioned tools will update the plugin tool constituting the application process easier and more comprehensive.

In addition, standardization of intrinsic vulnerability assessment in this complex geological medium will also make it possible to improve pollution risk assessment.

6.2. The Intrinsic Vulnerability Map of the Fontaine de Vaucluse

At the regional scale, the vulnerability map generally conforms to morphological units. At the local scale, the vulnerability map is primarily affected by karst features that seem to drive local hydrodynamics. It is important to note that almost half of the catchment area is characterized by high vulnerability. Even though the Fontaine de Vaucluse catchment area is not very affected by human activity (villages, farms, and airport), which occupies about 1% of the area and an agricultural zone (15%), these activities are located within areas of high vulnerability. The vulnerability map provides relevant support for sustainable resource management and development strategies. The vulnerability map should be used at the regional scale to enable coherent and sustainable aquifer management.