Next Article in Journal
Mechanism Analysis of Floor Water Inrush Based on Criteria Importance though Intercrieria Correlation
Previous Article in Journal
Assessment of Heavy Metal Distribution and Health Risk of Vegetable Crops Grown on Soils Amended with Municipal Solid Waste Compost for Sustainable Urban Agriculture
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Groundwater Vulnerability and Delineation of Protection Zones in the Discharge Area of a Karstic Aquifer—Application in Agyia’s Karst System (Crete, Greece)

School of Mineral Resources Engineering, Technical University of Crete, 73100 Chania, Greece
*
Author to whom correspondence should be addressed.
Water 2023, 15(2), 231; https://doi.org/10.3390/w15020231
Submission received: 11 December 2022 / Revised: 27 December 2022 / Accepted: 30 December 2022 / Published: 5 January 2023
(This article belongs to the Section Hydrogeology)

Abstract

:
This work represents a contribution to the protection techniques of karst aquifers against groundwater pollution. The paper sets out the methodology being introduced for the protection of the karstic system that gives rise to five (5) major groups of springs and supplies fourteen (14) pumping wells near Agyia Chania (Crete, Greece). Starting from a geological and hydrogeological survey of the area, the work presents a vulnerability assessment of the karstic aquifer based on the application of three index-based methods (EPIK, PRESK and DRISTPI). The protection zones for the discharge area of the aquifer were delineated through an integrated geomorphological approach and groundwater flow modeling. At first, the risk of polluting substances migration from ground surface to groundwater was considered based on the spatial distribution of vulnerability. Following this, the vulnerability was evaluated in the saturated zone, where the attenuation mechanisms of contaminants were reducing due to the raised flow velocity. The groundwater flow and contaminant transport processes was considered using the MODFLOW code. Next, the data from the vulnerability mapping and the groundwater flow simulation were merged into an integrated assessment to delimit the protection zones for the water abstraction points. The vulnerability assessment outlines zones of high vulnerability in the SE part of the area, far away from the discharge zone of the aquifer and the water abstraction points. These zones are associated with an intensive infiltration process via carbonate formations. Protection Zone I was delineated 20 m around the water abstraction points, and it should be excluded from any anthropogenic activity. Protection Zone II coves part of the very high and high vulnerability zones defined by the DRISTPI method (located upwards of the water abstraction points), as well as an area downwards of springs and wells, where the flow path lines which demonstrate the subsurface travelling time of 50 days are projected to the ground surface. Protection Zone III extends outside Zone Ι and Zone ΙΙ, up to the limits of the hydrogeological or hydrological basin, whichever is larger. It includes the entire capture zone (i.e., the surface and underground catchment area) that feeds the water abstraction points. In this manner the protection zones include the entire contributing area to water abstraction points, not just the ground surface recharge zone.

1. Introduction

In order to prevent groundwater pollution phenomena and to preserve water resources for human consumption, protection zones around the water abstraction points are established. By placing some form of regulatory control on activities taking place in these zones, the anthropogenic impact on the quality of the abstracted water can be minimized [1,2].
Protection zones can be delineated around the abstraction points using various approaches and criteria that range from relatively simple methods based on fixed distances, through more complex methods based on the travel time of water through the saturated zone, to sophisticated modelling approaches using log reduction models and contaminant kinetics [2]. Geographic information systems (GIS) have been used widely for the preparation of groundwater vulnerability maps and the delineation of protection zones [3,4]. Masud and El Osta [5] proposed the combined use of the DRASTIC index-based method and the GIS techniques as an effective method for groundwater pollution risk assessment, and for mapping the areas that are prone to the deterioration of groundwater quality and quantity. Furthermore, El Osta et al. [6] evaluated the hydrogeological and hydro chemical aspects of groundwater evolution and vulnerability by the integration of Geographic Information Systems (GIS), hydro chemical modeling, and the DRASTIC model.
Commonly, three protection zones are delineated around the abstraction points [1,7] in order to achieve the following levels of protection:
  • Zone I is immediately adjacent to the site of the well or spring to prevent the rapid ingress of contaminants.
  • Zone II refers to the protection against microbial contamination. The delineation of this zone is based on the time expected to be needed for dilution, effective attenuation of slowly degrading substances, or the reduction in pathogen presence to an acceptable level.
  • Zone III refers to the remote protection area against persistent chemical contaminants. It sometimes covers the entire catchment area of a particular water abstraction point.
In the European Union, the Water Framework Directive (WFD) 2000/60/EC [8,9,10] proposes the establishment of safety zones around water abstraction points for human consumption, aiming to avoid further quality degradation and to reduce the level of treatment required for the production of drinking water. However, most of the European countries have not yet incorporated the Water Framework Directive in their national legislation with specific instructions. As a result, various empirical practices are followed in the member states of the European Union that differ depending on the type of aquifer and the specific conditions of the region [11].
The protection of the groundwater resources in karst aquifers constitutes an extremely important issue. The transport of the pollutants towards the aquifer is extremely rapid, and their attenuation downgrades as a consequence of the hydrogeological properties [12]. The absence of a sufficiently thick soil cover, the existence of the epikarstic zone and the preferential movement of groundwater through cracks and karstic features which allow the development of high flow velocities result in the easy movement of contaminants to the aquifer and then to water abstraction points [13].
Several methods for the vulnerability mapping and the delineation of protection zones in karst aquifers have been proposed [14,15,16,17], but the issue still remains a challenge due to heterogeneity and anisotropy of their hydraulic conductivities which makes the prediction of groundwater flow a difficult task [13,18].
The European approach to vulnerability and risk mapping for the protection of carbonate (karst) aquifers proposed by COST Action 620 [19] is based on the “hazard-pathway-target model”, which applies for both resource and groundwater source protection. As mentioned by Daly et al. [20] the “pathway” includes everything between the ground surface (the point of contaminants release) and the “target” refers to the water which has to be protected [19,20]. It should be pointed out that all factors influencing vulnerability have to be evaluated.
The objective of the current work is to propose a methodological approach for delineation protection zones of the water abstraction points in the discharge zone of a karst aquifer. The work is based on a geomorphological-structural criterion (index-based method) that takes into account the potential infiltration of water from surface to groundwater, and the numerical modelling that simulates the groundwater flow and contaminant transport processes within the saturated zone.
Mapping the groundwater vulnerability was based on the hydrogeological, geological and geomorphological evaluation of the area. Three index-based methods were used for the vulnerability assessment considering the site characteristics and the available data. The applicability and the results of the methods were evaluated in relation to site characteristics in order to reduce the uncertainty and to provide reasonable results.
The groundwater flow and contaminant transport processes within the saturated zone was analyzed through the groundwater flow model constructed in MODFLOW-2000 under steady state conditions [21]. The time needed for reduction in the presence of pathogens, and the effective attenuation of slowly degrading substances to an acceptable level were considered for the delineation of protective zones around the water abstraction points.

2. Study Area

The study area is located in the northern part of Chania Prefecture (Crete, Greece), between latitudes 35°18′48″ and 35°32′35″ N and longitudes 23°50′14″ and 24°08′42″ E (Figure 1). It extends from the foothills of the Lefka Ori Mountains in the South, to the lowland zone on the North side of the island (Figure 1), and includes the Agyia’s springs catchment area of 149 km2.
The karstic aquifer hosted in the study area is recharged mainly through the carbonate rocks located on the South and Southeastern side of the region, as well as through the infiltration of precipitation in the entire area.
Five major groups of springs (Kalamionas, Kolymba, Vrysidia, Platanos, Varypetro) near the Agyia’s settlement (Figure 1) comprise the main discharge outlets for the karst aquifer [21,22]. They constitute a front of springs with a total water flow that ranges from 1.7 to 3.4 m3/s, with an average value of 2.5 m3/s.
Beyond the exploitation of the aforementioned springs, a series of pumping wells has been constructed in the exploitation fields of Agyia and Myloniana (upstream of the discharge zone of the aquifer), (Figure 1) in order to cover the water demands in the wider region. North of Agyia’s springs, there is also a homonymous (man-made) lake with rich biodiversity.

2.1. Physiography and Geological Settings

The relief of the area is the result of strong tectonic movements and karstification. The terrain rises gradually from the sea level in the north to almost +2100 m a.s.l. in the south (Lefka Ori Mountains, Figure 1). The ground surface inclination ranges from zero in the lowlands, to 70% in the semi-mountainous region. Moreover, the Southern parts of the area show minimal vegetation that is mainly bushy, while in the North, at low altitudes, there is intense vegetation and cultivated regions.
The lithological and structural mapping of the region [23,24,25,26] has shown that the southern part of the area is covered by formations of the Plattenkalk series. It includes crystalline limestones and dolomites, interbedded with watertight horizons of chert and schists (Figure 2); they dip about 50–60° to the NW (55°/320 NW).
The Trypali carbonate unit significantly spreads in the South and central part of the area (Figure 2). It consists of limestones, dolomites and crystalline limestones locally associated with coarse cobblestone or cellular dolomites.
The carbonate formations of the Trypali unit are tectonically overlain by the Phyllite-Quartzite series that outcrops along the S-SW margins of the study area (Figure 2).
The uppermost Quaternary sediments, dominated by alluvial, marine terrace deposits, talus cones and scree, overlie the older geologic formations mainly throughout the Northern part of the study area (Figure 2).
The tectonic regime of the region is characterized by normal faults of E-W strike, which locally deviate to ENE-WSW and ESE-WNW (Figure 2).

2.2. Hydrologic Hydrogeological Setting

The climate of the area is Mediterranean, with mild and rainy winters and warm, dry summers [27]. The wintry period lasts from October to March, and the summer period lasts from April to September. The annual mean temperature is 18.7 °C. The lowest mean temperature is 11.5 °C in January, and the highest is 26.6 °C in July. Rainfall occurs mainly during winter.
Mourkakou [28] studied monthly rainfall data collected from seven stations located in the area for the period 2006–2018. The average annual precipitation at each station was calculated and plotted against the altitude of that station (Figure 3).
A regression line fitted to these data indicates that the relationship between the rainfall and the altitude is given by the Equation (1):
P = 2.065 h + 478.38,
where, P is the precipitation in mm, and h the altitude in m.
Based on the aforementioned relationship, the average annual precipitation in Agyia’s catchment area is estimated at 260 × 106 m3. Thirty seven percent (37%) of this quantity, which is equal to 96 × 106 m3, infiltrates, and twenty percent (20%), namely 53 × 106 m3, runs off. The mean annual evaporation ranges between 43–58% of total annual precipitation [28].
Plattenkalk formations are characterized by low permeability values (10−7 to 10−5 m/s), and where the geological conditions allow, they consist of the hydrogeological base layer of the overlying Trypali carbonate rocks [21,22,29].
The highly fractured nature of the Trypali unit, in addition to the dissolution and the erosion process, have created a great variety of karst features, from simple rillenkarren up to dolines and karst cavities. The surface water tends to concentrate into the aforementioned features and then infiltrates, recharging the groundwater [22,29].
The hydraulic properties of this formation were estimated by pumping tests in the Myloniana well-field [22]. The transmissivity of the carbonate unit ranges between 0.1 to 1 m2/s, which is equivalent to a permeability between 10−3 and 10−2 m/s.
Pertinent studies show that 55% of the precipitation infiltrates, recharging the groundwater, while surface run off is very limited (about 5%), and surface drainage occurs only after strong rainfalls [30].
Furthermore, the analysis of hydrographs of springs and wells within the karstic aquifer hosted in the Trypali unit presents a corresponding seasonal variation of the groundwater level [29] that suggests a homogeneous karst.
Based on the aforementioned annotations, it is assumed that the Trypali unit forms a separate hydrogeological entity (Figure 2 and Figure 4) that discharges through the Agyia’s springs while the Plattenkalk unit forms the accompanying base layer.
Phyllite-Quartzite formation is characterized as almost impermeable (k < 10−8 m/s) due to its mineralogical composition and structure and selective circulation (10−8 to 10−7 m/s), with low to very low water permeability, when interferences of quartzite or crystalline limestones exhibit higher growth [21]. The infiltration coefficient of this formation is estimated to be approximately 5%, and the runoff coefficient reaches an amount of 35% [30,31].
The uppermost Quaternary sediments are of medium-low permeability (10−7 to 10−5 m/s) and they host phreatic aquifers of not significant yield compared to the karst aquifer hosted in the Trypali unit. The infiltration coefficient in the alluvial sediments is estimated to be about 20%, while the runoff coefficient is estimated to be about 20–30% [30,31].
The tectonic structure determines the area that contributes water to the karst aquifer (Figure 2 and Figure 4). Specifically, the Phyllite-Quartzite unit in the west constitutes a regional aquiclude (no flow boundary); to the North, the stratigraphic (locally tectonic) limit between the Quaternary formations and the Phyllite-Quartzite unit determines an impermeable boundary that holds the aquifer and prevents the flow of the karst water towards the sea; in South the groundwater exchanges with the rest of the hydrogeological system of the Lefka Ori Mountains, while in the Eastern part of the area, data resulting from geological surveys and well drillings indicate that the Plattenkalk series (aquiclude) has a structural top at an altitude of +625 m a.s.l. (Figure 4C). It acts as a local watershed (Figure 2) and divides the groundwater supplying the Agyia’s springs in the West, from the groundwater supplying the Chania basin in the East. Based on this conceptual model, Agyia’s springs are fed by the groundwater flowing in the Eastern part of the catchment area.
The system is recharged mainly through the carbonate rocks exposed in the south-eastern and southern part of the area (Lefka Ori Mountains), as well as through precipitation infiltrating from the ground surface to the groundwater table, across the entire area. The watertight horizons of chert and schists in the Plattenkalk series restrict the percolation of water and strongly influence its movement and concentration in the aquifer hosted in the Trypali unit.
The results of the systematic piezometric surveys [30,32] confirm a regional groundwater flow direction from SE-S towards NW–N (Figure 2), with hydraulic gradients ranging from 3 to 8‰. In the North margin of the study area (Figure 2 and Figure 5), a fault zone running SW–NE approximately parallel to the coast, sinks the Phyllite-Quartzite nappe lower than the sea level and prevents the flow of the karst water to the sea. The direction of the fault identifies with the axis of occurrence of Agias’s group of springs (Figure 2), which comprise the main discharge of the aquifer. They are contact-overflow springs with a mean outflow rate of 60 × 106 m3 annually; they outflow at an altitude from +33.5 m a.s.l. (Kalamionas) to +40.6 m a.s.l. (Varypetro), [29].
The analyses of groundwater indicate good water quality that is safe for human consumption, with a temperature almost constant at 13 °C [33].
In addition to springs exploitation, fourteen (14) pumping wells have been constructed in the highland area (upstream of the aquifer’s discharge zone) (Figure 2 and Figure 4D), providing water for domestic and agricultural use in the wider area. The total amount of water exploitation is estimated at about 85 × 106 m3.

3. Vulnerability of the Karst Aquifer

The vulnerability of karstic systems depends mainly on the residence time of water in different sections of the aquifer (endokarst, epikarst and protective soil cover) [14]. Specifically, in endokarst or saturated karst, a well-developed network of conduits (where the flow velocity is high) implies high vulnerability. Also, the more directly the epikarst is connected to the karstic conduits network, the greater the vulnerability. Regarding the protective soil cover, the residence time of water and vulnerability depends on the permeability and thickness of the cover. For the vulnerability assessment of these aquifers, index-based methods are used, among others. They take into account the protective effect of the unsaturated zone, but they attach less importance to the groundwater transit times.
In the current work, the EPIK, PRESK and DRISTPI methods were applied. Initially, similar to Iván and Mádl-Szőnyi [34], the relationships between the methods were revealed in order to understand their innovations, advantages and drawbacks, as well as the data need. The results from the different methods were compared and the groundwater vulnerability was assessed, with emphasis on water infiltration from the ground surface to the saturated zone. The characteristics of the applied methods area is as follows:

3.1. The EPIK Method

The EPIK method was developed by Dörfliger and Zwahlen [35] for the study of karstic aquifers in Switzerland. Its name comes from the initials of the four parameters that are considered: the degree of development of the epikarst (E), the characteristics of the protective ground cover (P), the conditions of infiltration (I), and the development of the karst network (K).
Each of the aforementioned variables is graded with a value from 1 to 3 or 4, depending on the degree of protection it provides from pollution.
The protection factor (Fp) is calculated from Equation (2):
Fp = 3·E + 1·P + 3·I +2·K
The Fp value can range from 9 to 34. High Fp values indicate increased physical protection and therefore a low vulnerability value of the aquifer [35,36,37].

3.2. The PRESK Method

This method was developed by Koutsi and Stournaras [38], and it is an adaptation of the RISKE method in the Mediterranean region. It evaluates the aquifer protection (P), the rock type (R), the epikarstic zone (E), the soil cover (S), and the development of the karst network (K). Each factor is scored from 0 to 3, with values corresponding to very low, low, moderate and high vulnerability. The final value of the vulnerability factor (Ig) is calculated from Equation (3):
Ig = 0.263·P + 0.097·R + 0.160·E + 0.419·S + 0.062·K
Regarding the parameter P (variable of the aquifer protection), it is defined as the sum of the ground surface inclination (I) and vegetation cover (V) indices, multiplied by weight factors according to Equation (4):
P = 0.8·I + 0.2·V
The ground surface inclination (I) is directly related to the amount of water that infiltrates. With reference to the indicator of vegetation cover (V), it depends on the presence and density of the flora and affects the value of evapotranspiration.
High vulnerability factor values (Ig) indicate the reduced natural protection of the aquifer and correspond to a high vulnerability value.

3.3. The DRISTPI Method

The DRISTPI method was developed by Jiménez Madrid et al. [39,40] for karst aquifers in Spain and it is based on the older DRASTIC method [41,42]. It is worth mentioning that although the DRASTIC method was created for the evaluation of intrinsic vulnerability in all types of aquifers, it does not take into account the specific characteristics of the karst systems. Moreover, the major drawback of the DRASTIC model is the combination of the vulnerability of the unsaturated zone with the susceptibility of the groundwater once the pollutant has reached the water table [43].
In contrast to the DRASTIC method, the DRISTPI model is applicable to both karstic and non-karstic aquifers. The acronym DRISTPI is derived from the initials of the parameters evaluated in the method, namely: the depth of the potentiometric surface from ground surface (D), the recharge (R) of the aquifer, the unsaturated zone (I), the type of soil formation (S), the topography/ground surface inclination (T), and the privileged infiltration (PI). As is obvious, the parameters of the DRASTIC method referring to aquifer lithology and hydraulic conductivity have been removed. The introduced parameter PI concerns those areas where rapid infiltration of water into the aquifer occurs (e.g., dolines), making it more vulnerable to pollution.
Depending on the protection that each parameter provides to the system, the corresponding factor is evaluated from 1 (very high protection) to 10 (very low protection).
The vulnerability index (Id) is derived from the sum of the products of the aforementioned six parameters, multiplied with the appropriate weight factors. The final value of Id ranges from 17 to more than 162. It is noted that the parameters S, T, I and PI are estimated according to the principles of the DRASTIC method, while the parameter D (aquifer’s depth) differs depending on whether the rocks are karstified (scenario 1) or not (scenario 2). For scenario 1 (high karstified rocks), the vulnerability index (Id) is defined by Equation (5)
Id = 2·D + 4·R + 5·I + 2·S + T+ 5·PI,
and for scenario 2 (not karstified rocks, or detritic material) Equation (6) is used
Id = 5·D + 4·R + 5·I + 2·S + T+ 5·PI
High Id values (over 110 for scenario 1 and over 127 for scenario 2) indicate the reduced physical protection of the aquifer, and result in high vulnerability values.

3.4. Vulnerability Assessment in Study Area

The vulnerability assessment in the study area was based on the field survey, guided by the soil map (Figure 6), constructed using the geographical maps of the land area, published by the Greek Forest Service [44,45].
Based on the geomorphology, the soil depth, the erosion, the ground surface slope, the vegetation and the degree of anthropogenic impact in the study area, each parameter of the EPIK, PRESK and DRISTPI methods was evaluated separately. By applying Equations (2)–(6), the protection and vulnerability factors across the study area were calculated.
Following, the area was divided into five (5) vulnerability zones (Table 1 and Table 2), according to the requirement of the Greek environmental authorities [46]. The vulnerability indices evaluated by the EPIK, PRESK and DRISTPI methods are given in Table 1 and Table 2, and the distribution of vulnerability in study area is depicted in the maps presented in Figure 7.
The applied index-based methods show that the zones with the lowest degree of protection are located in the areas where water infiltration is favoured by fragmentation and intensive karstification (high and very high vulnerability zones). These areas concern the Northern foothills of the Lefka Ori Mountains, up to the exploitation field of Myloniana (Figure 2 and Figure 7). The vulnerability map produced through EPIK method was similar to those obtained from the PRESK and DRISTPI methods (Figure 7). However, the EPIK method (Figure 7a) indicates that vulnerability is affected mainly by the characteristics of the protective cover and to a lesser extent by the presence of fractured zones.
Moreover, the PRESK method (Figure 7b) differs from the EPIK one (Figure 7a), as the latter identifies higher vulnerability in the areas less sensitive to pollution. In detail, the T1 and T2 vulnerability class, according to EPIK, changes to the T2 and T3 class, respectively, according to the PRESK method.
In contrast, the DRISTPI method (Figure 7c) degrades the vulnerability of non-karstic formations (T2 class according to EPIC changes to T1 class according to DRISTPI), and upgrades the vulnerability class of the karst region (T3 and T4 classes change to T4 and T5, respectively).
Overall, the results of the DRISTPI method and their difference from those derived by the EPIK and PRESK methods are presented in Table 1 and Table 2 and Figure 7 and Figure 8. The extent of the high and very high vulnerability zones according to the DRISTPI method predominates with a percentage of 53% of the entire area, while according to the EPIK and PRESK methods it is limited to about 22% (Table 1 and Table 2, Figure 8).
Moreover, using the DRISTPI method, the low vulnerability class is wider than by using the other two methods.
From the comparison of the results (Figure 8), it is also concluded that the DRISTPI method emphasizes the vulnerability of karst formations by classifying them in high and very high vulnerability zones, while they are characterized as moderately vulnerable zones by the EPIK and PRESK methods.
Based on the aforementioned annotations and the methodologies proposed in COST Action 620 [19], which, as referred to by Jiménez-Madrid et al. [39] are in accordance to the philosophy of the WFD to protect water resources, the DRISTPI index was chosen as the proper method to assess the vulnerability to contamination of the study area.

4. Groundwater Flow Model in Aquifer’s Discharge Zone

Assuming that the local heterogeneities of the karst in the study area are smoothed out due to the karstification extent, the system was modelled as an equivalent porous medium.
The groundwater flow was simulated using the MODFLOW code [47] and the Visual MODFLOW v.2010.1 graphical interface [48]. The model was analyzed extensively in the discharge zone of the aquifer, assuming steady state flow conditions. The modelling domain extends to 20.7 km2 in the northern part of Agyia’s hydrological basin (Figure 2), at altitudes from +35 m a.s.l. to +400 m a.s.l. Alluvial, Diluvial and Neogene formations appear on the surface, while the underlying carbonate rocks of the Trypali unit outcrop in the south eastern part of the model area.
The aquifer hosted in Trypali carbonate formations was characterized as unconfined with the possibility of transition to confined conditions, depending on the position of the piezometric surface [21].
The appropriate set of boundary conditions were assigned according to hydrogeologic data of the region. In particular, the constant hydraulic head equal to +15 m a.s.l. was defined at Agyia’s lake location, while at the SE edge of the model where the main supply of the aquifer originates, the hydraulic head equal to +86 m a.s.l. was assigned (Figure 10). Furthermore, all significant groundwater intakes were incorporated in simulation modelling. The calibration procedure focused on choosing model input parameters (infiltration rate and hydraulic head in the boundaries) until the computed results (groundwater level and springs discharge) matched the observed data reasonably well.
Overall, the model’s performance concerning the measured and the estimated yield of the springs was acceptable, with a divergence lower than 1%. Moreover, the calculated mean percent error of the groundwater level in relation to the field data was equal to −0.51 m, indicating an acceptable degree of accuracy (Figure 9).
The produced potentiometric surface map (Figure 10) indicates that the main path flow of the groundwater is generally from the Southeast towards the Northwest. Agyia’s springs are supplied mainly by the groundwater that flows in the Trypali carbonate formation, without finding significant escapes north of the springs’ outlet. According to the results, the hydraulic gradients estimated along the flow path range from 1 to 5‰.
The numerical simulation of groundwater flow was also applied to define the contribution area and the “ransfer time zones” for each abstraction point. Initially “targets” (particles) were defined at the locations of the water points (pumping wells and springs). Next, the “Modpath” subroutine of MODFLOW was applied to generate flow path lines to the targets (Figure 11). The path lines delineate the groundwater flow path and the time required for the elementary fraction of water to reach the designated target.
Using the groundwater flow lines of 360-day transit time (shown in Figure 11) and by connecting points at which pollution arrives at the same time, different capture zones were delineated. The isochronal curves of 50-days (Figure 11) extend up to 2 km along the preferential groundwater flow paths (coming from the SSE to the NW).

5. Delineation of Protection Zones

According to Cost Action 65 [49], source protection zones should be defined around individual groundwater sources of supply, which ideally encompass the total catchment area, while the groundwater resources should be protected with a more extensive zoning that may cover the entire land surface of the aquifer (resource protection zones). The produced map divides the area into a number of groundwater protection zones according to the degree of protection required. Within these protection zones, a code of practice describes the recommended controls for both existing and new activities.
In fragmented and karstified aquifers, several methods have been proposed for the delineation of protection zones, but the process remains a challenge due to the heterogeneity and anisotropy of such systems [34,50,51].
As reported by Biava et al. [13] and Pochon et al. [51], the methodologies for the delineation of groundwater protection zones in such cases must generally be adapted to the characteristics of the site and the available data. Usually, for delineation of the protection zones, the isochronous groundwater movement in the saturated zone is considered, without taking into account the contaminant infiltration through the overlying layers. This limitation is important when combined with the great depth of the water table, because it may lead to unreliable results.
Based on the aforementioned views, in order to define the groundwater protective zones around the water abstraction points in the study area, the risk of pollution was investigated both:
  • in the unsaturated zone through which contaminants infiltration may occur, and
  • in the saturated zone, where the attenuation mechanisms of contaminants do not treat the polluted environment due to the rise in flow water velocity.
Firstly, the distribution of the aquifer’s vulnerability was assessed through proper index test methods, and the result was used to define the sensitive zones that must be protected from contaminant infiltration through the overlying layers of the aquifer.
Next, the isochrone method was applied through the numerical simulation of groundwater flow in order to determine the extent of the region that may be affected by contaminant transport through the subsurface flow in the saturated zone.
The results from the vulnerability estimation and the groundwater flow simulation were merged into an integrated assessment in order to delimit the protection zones for the groundwater abstraction points located in the discharge zone of the aquifer. In detail, the following protection zones were designed (Figure 12):
Protection Zone I: This represents an immediate protection zone against the direct introduction of contaminants to the groundwater. It was delimitated at 20 m around water abstraction points based on the introduced guidelines [2,16,46,51] and it should be excluded from any anthropogenic activity.
Protection Zone II: This is known as a controlled zone or a microbiological protection zone. It spans in the area around Zone I and it should be protected against pollution from any anthropogenic activity, such as crafts, oil mills, cesspools, cemeteries, etc. It is a sensitive area that is necessary to be guarded on a 24-h basis.
For delineating this zone, the water flow around the water abstraction points was simulated under steady flow conditions. It should be noted that the modeling domain extends up to the boundaries of the hydrogeological or hydrological basin, whichever is larger (Figure 12).
Like Lipfert et al. [50], the authors found it necessary to protect the entire contributing area, not just the underground recharge zone. For this reason, the projection area of the groundwater path lines to the surface was considered. As is shown in Figure 12, the zone includes part of the very high and high vulnerability zone defined by the DISTRPI method (located upwards of the water abstraction points), as well as the area extending downwards of the wells and springs, to which the flow path lines representing the 50 days travel time are projected to the ground surface (Figure 11). In this manner, the protection zone includes the entire contributing area to water abstraction points, not just the ground surface recharge zone.
Protection Zone III: This refers to the area of the unsaturated zone of the karstic system that can transport the pollution or the chemical contaminant to the downstream area where the springs and pumping wells are located (Figure 12).
This zone extends outside Zone Ι and Zone ΙΙ, up to the limits of the hydrogeological or hydrological basin, whichever is larger. It includes the entire capture zone (i.e., the surface and underground catchment area) that feeds the water abstraction points. In the study area, the divergence between the boundaries of the hydrologic and hydro geologic basin is depicted by the diagonal lines in Figure 12.

6. Discussion

The reliability of the produced vulnerability map was verified with a procedure similar to that proposed by Różkowski [52], based on the geological knowledge, the results of hydrogeological investigations [22,32,53], and the chemical and isotopic analyses of the groundwater [54]. The delineation of protection zones in the discharge area of the karst aquifer was conducted through a geomorphological-structural approach (index-based method) that takes into account the potential infiltration of water from surface to groundwater and a numerical modelling that simulates the groundwater flow and the contaminant transport processes within the saturated zone. The implementation of this approach provides a valuable tool to establish proper protective measures against groundwater pollution in karst areas.
Considering the contaminant infiltration through the overlying layers of the aquifer as well as the groundwater flow simulation in the saturated zone, three protection zones were designed to prevent groundwater pollution phenomena.
Particularly for the delineation of the microbiological protective zone (Protection Zone II), the risk of pollution was investigated both in the unsaturated as well in the saturated zone of the aquifer, where the attenuation mechanisms of contaminants do not treat the polluted environment satisfactorily as the water flow velocity increases. The data from vulnerability mapping and the groundwater flow simulation were merged in order to delimit the microbiological protection zone for each groundwater abstraction point.
In order to further improve the vulnerability assessment of the karst aquifer and to improve the confidence in the delimitation of protection zones, more data concerning the groundwater quality and isotopic investigations are required. However, the tracer test data are rarely available; and as Marín and Andreo [55] have mentioned the analyses of hydrochemical and hydrodynamic responses concerning the total organic carbon content is an alternative that may be of use for the further validation of the produced map.

7. Conclusions and Recommendations

The results from the used methods to assess the intrinsic vulnerability of the aquifer were compared and the groundwater vulnerability was assessed with an emphasis on water infiltration from the ground surface to the saturated zone.
Using the EPIK and the PRESK methods to assess the intrinsic vulnerability of the aquifer, it is concluded that the aforementioned methods only take into account the infiltration through overlying layers (soil and unsaturated zone), without considering the groundwater flow in the saturated zone. On the contrary, the DRISTPI method considers the water flow in deeper layers, even if only indirectly through the piezometric contours and the hydraulic parameters of the aquifer that are ignored in the other two methods. Based on the aforementioned annotations, the DRISTPI index is proposed for assessing the vulnerability to contamination in karst areas, as it emphasizes the behavior of the carbonate formations that constitute zones most vulnerable to pollution.
The vulnerability assessment in the study area outlines zones of high vulnerability in the SE part of the area, far away from the discharge zone of the aquifer and the water abstraction points. These zones are associated with an intensive infiltration process via carbonate formations.
Protection Zone I was delimited 20 m around the water abstraction points based on the introduced guidelines [2,7,16,46,51], and it should be excluded from any anthropogenic activity.
Protection Zone II includes part of the very high and high vulnerability zones defined by the DISTRPI method (located upwards of the water abstraction points), as well as an area extending downwards of the wells and springs, in which the flow lines of 50 days transit time are projected to the ground surface. In this manner, the protection zone includes the entire contributing area to water abstraction points, not just the ground surface recharge zone.
Finally, Zone ΙII was delimitated outside Zone Ι and Zone ΙΙ, up to the limits of the hydrogeological or hydrological basin, whichever is larger. It includes the entire capture zone (i.e., the surface and underground catchment area) that feeds the water abstraction points.
It is clear that for the proper delineation of the protection zones, coupling the results from the geological structure survey and the groundwater vulnerability assessment is necessary. In any case, the investigation should extend to the boundaries of the hydrogeological or hydrological basin, whichever is larger, in order to cover the entire area of the groundwater catchment.
The implementation of this approach that involves the investigation of pollution risk in the entire area of the groundwater catchment provides a valuable tool for the protection of the discharge area of a karst aquifer. The result enables the groundwater management agency to establish protective measures against pollution in the vulnerable zones and to preserve the high quality of water for human consumption.

Author Contributions

Conceptualization, E.S.; methodology, E.S. and D.V.; software, D.V.; validation, E.S., D.V. and O.M.; formal analysis, E.S.; data curation, D.V. and E.S.; supervision, E.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

Part of this work has been presented in the 12th Hydrogeological Conference: Groundwater in a changing environment, Nicosia, March 2022.

Conflicts of Interest

The authors declare that they have no conflict of interest.

References

  1. WHO (World Health Organization). Protecting Groundwater for Health: Managing the Quality of Drinking-Water Sources; Schmoll, O., Howard, G., Chilton, J., Chorus, I., Eds.; IWA Publishing: London, UK, 2006; ISBN 18433100795.
  2. Chave, P.; Howard, G.; Schijven, J.; Appleyard, S.; Fladerer, F.; Schimon, W. Groundwater protection zones. In Protecting Groundwater for Health: Managing the Quality of Drinking-Water Sources; Schmoll, O., Howard, G., Chilton, J., Chorus, I., Eds.; IWA Publishing: London, UK, 2006; pp. 465–492. [Google Scholar]
  3. Özdemir, A. Determination of protection zones in drinking water basins: A case study from Turkey, Sapanca Lake Basin. Environ. Earth Sci. 2020, 79, 178. [Google Scholar] [CrossRef]
  4. Deh, S.K.; Kouame, K.J.; Eba, A.E.; Djemin, J.E.; Kpan, A.; Jourda, J.P. Contribution of geographic information systems in protection zones delineation around a surface water resource in Adzope Region (Southeast of Cτte d’Ivoire). J. Environ. Prot. 2017, 8, 1652–1673. [Google Scholar] [CrossRef] [Green Version]
  5. Masoud, M.H.; El Osta, M.M. Evaluation of groundwater vulnerability by using modeling and GIS techniques in El-Bahariya Oasis—Western Desert—Egypt. J. Earth Syst. Sci. 2016, 125, 1139–1155. [Google Scholar] [CrossRef] [Green Version]
  6. El Osta, M.; Niyazi, B.; Masoud, M. Groundwater evolution and vulnerability in semi-arid regions using modeling and GIS tools for sustainable development: Case study of Wadi Fatimah, Saudi Arabia. Environ. Earth Sci. 2022, 81, 248. [Google Scholar] [CrossRef]
  7. Environment Agency-National Groundwater and Contaminated Land Centre. In Proceedings of the Protecting Groundwater: An International Conference on: Applying Policies and Decision-Making Tools to Land-Use Planning, Birmingham, UK, 4–5 October 2001.
  8. European Parliament; Council of the European Union. Directive 2000/60/EC of the European Parliament and of the Council of 23 October 2000 Establishing a Framework for Community Action in the Field of Water Policy; Official Journal of the European Communities: Brussels, Belgium, 2000.
  9. European Parliament; Council of the European Union. 2006/118/EC of the European Parliament and of the Council of the 12 December 2006 on the Protection of Groundwater Against Pollution and Deterioration; Official Journal of the European Communities: Brussels, Belgium, 2006.
  10. European Commission. Directorate-General for Environment, Guidance on Groundwater in Drinking Water Protected Areas, Guidance Document No 16; Office for Official Publications of the European Communities: Luxembourg, 2006.
  11. Voudouris, K.; Kazakis, N. Vulnerability assessment of Karst aquifers: The case of Mitsikeli, Ioannina. In Proceedings of the 14th Conference of the Hellenic Hydro Technical Union (EYE), Volos, Greece, 16–17 May 2019. [Google Scholar]
  12. Kacaroglu, F. Review of groundwater pollution and protection in karst areas. Water Air Soil Pollut. 1999, 113, 337–356. [Google Scholar] [CrossRef]
  13. Biava, F.; Consonni, M.; Francani, V.; Gattinoni, P.; Scesi, L. Delineation of protection zones for the main discharge area of the Gran Sasso Aquifer (Central Italy) through an integrated geomorphological and chronological approach. J. Water Resour. Prot. 2014, 6, 1816–1832. [Google Scholar] [CrossRef] [Green Version]
  14. Dörfliger, N.; Jeannin, P.Y.; Zwahlen, F. Water vulnerability assessment in karst environments: A new method of defining protection areas using multi-attribute approach and GIS tool (EPIK method). Environ. Geol. 1999, 39, 165–172. [Google Scholar] [CrossRef] [Green Version]
  15. Civita, M.V. An improved method for delineating source protection zones for karst springs based on analysis of recession curve data. Hydrogeol. J. 2008, 16, 855–869. [Google Scholar] [CrossRef]
  16. Bussard, T.; Tacher, L.; Parriaux, A.; Maître, V. Methodology for delineating groundwater protection areas against persistent contaminants. Q. J. Eng. Geol. Hydrogeol. 2006, 39, 97–109. [Google Scholar] [CrossRef]
  17. Nguyet, V.T.M.; Goldscheider, N. A simplified methodology for mapping groundwater vulnerability and contamination risk, and its first application in a tropical karst area, Vietnam. Hydrogeol. J. 2006, 14, 1666–1675. [Google Scholar] [CrossRef]
  18. Berkowitz, B. Characterizing flow and transport in fractured geological media: A review. Adv. Water Resour. 2002, 25, 861–884. [Google Scholar] [CrossRef]
  19. Zwahlen, F. Vulnerability and Risk Mapping for the Protection of Carbonate (Karst) Aquifers Final Report (COST Action 620); European Commission, Directorate-General XII Science, Research and Development: Brussels, Belgium, 2004. [Google Scholar]
  20. Daly, D.; Dassargues, A.; Drew, D.; Dunne, S.; Goldscheider, N.; Neale, S.; Popescu, I.C.; Zwahlen, F. Main concepts of the “European approach” to karst-groundwater-vulnerability assessment and mapping. Hydrogeol. J. 2002, 10, 340–345. [Google Scholar] [CrossRef] [Green Version]
  21. Steiakakis, E.; Vavadakis, D.; Kritsotakis, M. Simulation of springs discharge from a karstic aquifer (Crete, Greece), using limited data. Environ. Earth Sci. 2015, 74, 4303–4315. [Google Scholar] [CrossRef]
  22. Steiakakis, E.; Monopolis, D.; Vavadakis, D.; Manutsoglu, E. Hydrogeological Research in Trypali Carbonate Unit (NW Crete). In Advances in the Research of Aquatic Environment; Springer: Berlin/Heidelberg, Germany, 2011; pp. 561–567. [Google Scholar] [CrossRef]
  23. Martini, H.J. Geological Map of Greece, Platanias Sheet (Scale 1:50,000); IGME: Athens, Greece, 1956. (In Greek) [Google Scholar]
  24. Tataris, A.; Christodoulou, G.E. Geological Map of Greece, Alikianos Sheet (Scale 1:50,000); IGME: Athens, Greek, 1969. (In Greek) [Google Scholar]
  25. Karageorgiou, E. Geological Map of Greece, Chania Sheet (Scale 1:50,000); IGME: Athens, Greece, 1971. (In Greek) [Google Scholar]
  26. Vidakis, M.; Triantaphyllis, M.; Mylonakis, I. Geological Map of Greece, Vryses Sheet (Scale 1:50,000); IGME: Athens, Greece, 1993. (In Greek) [Google Scholar]
  27. Hellenic National Meteorological Service. The Climate of Greece. Available online: http://www.hnms.gr/emy/en/climatology/climatology (accessed on 10 December 2022).
  28. Mourkakou, O. Contribution to the Assessment of the Water Balance of the Karstic Aquifer of Agyia (Chania). Diploma Thesis, Technical University of Crete, Chania, Greece, 2018. (In Greek). [Google Scholar]
  29. Steiakakis, E.; Vavadakis, D.; Kritsotakis, M. Karstic springs of Agyia and attempted restoration: Possibilities and problems. In Proceedings of the 11th International Hydrogeological Conference, Athens, Greece, 4–6 October 2017; Volume 1, pp. 435–444. [Google Scholar]
  30. Lionis, M.; Perleros, B. Hydrogeological Study for Chania Area (KA 9481721); Ministry of Agricultural, Directorate of Hydrogeology: Athens, Greece, 2001. (In Greek) [Google Scholar]
  31. Perleros, V.; Papamastorakis, D.; Kritsotakis, M.; Drakopoulou, E.; Panagopoulos, A. Groundwater potential of the island of Crete. Problems and perspectives. Bull. Geol. Soc. Greece 2004, 36, 2048–2056. (In Greek) [Google Scholar] [CrossRef] [Green Version]
  32. Ydrogaia. Hydrogeological Study of the Area of Agyia (Chania); Ministry of Public Works: Athens, Greece, 1977. (In Greek) [Google Scholar]
  33. Pavlidou, S. GR 1303: Karstic aquifer system of the Lefka Ori Mountains. In 3rd CSF—Operational Program: Competitiveness, Project: Recording and Evaluation of the Hydrogeological Characteristics of the Country’s Groundwater and Aquifer Systems (K.E. 7.3.2.1); Institute of Geological and Mineral Research (IGME), Department of Water Resources and Environment: Rethimnon, Greece, 2009. (In Greek) [Google Scholar]
  34. Iván, V.; Mádl-Szőnyi, J. State of the art of karst vulnerability assessment: Overview, evaluation and outlook. Environ. Earth Sci. 2017, 76, 112. [Google Scholar] [CrossRef]
  35. Dörfliger, N.; Zwahlen, F. Groundwater Vulnerability Mapping in Karstic Regions (EPIK), Practical Guide; Swiss Agency for the Environment, Forests and Landscape (SAEFL): Bern, Switzerland, 1998; p. 56. [Google Scholar]
  36. Vogelbacher, A.; Kazakis, N.; Voudouris, K.; Bold, S. Groundwater vulnerability and risk assessment in a karst aquifer of Greece using EPIK method. Environments 2019, 6, 116. [Google Scholar] [CrossRef] [Green Version]
  37. Nekkoub, A.; Baali, F.; Hadji, R.; Hamed, Y. The EPIK multi-attribute method for intrinsic vulnerability assessment of karstic aquifer under semi-arid climatic conditions, case of Cheria Plateau, NE Algeria. Arab. J. Geosci. 2020, 13, 709. [Google Scholar] [CrossRef]
  38. Koutsi, R.; Stournaras, G. Groundwater vulnerability assessment in the Loussi polje area, N Peloponessus: The PRESK method. In Advances in the Research of Aquatic Environment; Lambrakis, N., Stournaras, G., Katsanou, K., Eds.; Environmental Earth Sciences; Springer: Berlin/Heidelberg, Germany, 2011; pp. 335–342. [Google Scholar] [CrossRef]
  39. Jiménez-Madrid, A.; Carrasco-Cantos, F.; Martínez-Navarrete, C. Protection of groundwater intended for human consumption: A proposed methodology for defining safeguard zones. Environ. Earth Sci. 2012, 65, 2391–2406. [Google Scholar] [CrossRef]
  40. Jiménez Madrid, A.; Carrasco, F.; Martinez, C.; Gogu, R. DRISTPI A new groundwater vulnerability mapping method for use in karstic and non-karstic aquifers. Q. J. Eng. Geol. Hydrogeol. 2013, 46, 245–255. [Google Scholar] [CrossRef]
  41. Aller, L.; Bennett, T.; Leher, J.H.; Petty, R.J.; Hackett, G. DRASTIC: A Standardized System for Evaluating Groundwater Pollution Potential Using Hydrogeologic Settings; EPA/600/2-85/018; US Environmental Protection Agency: Washington, DC, USA, 1987. [Google Scholar]
  42. Sinan, M.; Razack, M. An extension to the DRASTIC model to assess groundwater vulnerability to pollution: Application to the Haouz aquifer of Marrakech (Morocco). Environ. Geol. 2009, 57, 349–363. [Google Scholar] [CrossRef]
  43. Barbulescu, A. Assessing Groundwater Vulnerability: DRASTIC and DRASTIC-Like Methods: A Review. Water 2020, 12, 1356. [Google Scholar] [CrossRef]
  44. Anestis, G.; Ziagas, E.; Nakos, G. Geographical Map of Lands, Sheet Vryssai. Scale 1:50,000; Ministry of Agriculture, Forest Protection Service: Athens, Greece, 1997. (In Greek) [Google Scholar]
  45. Schinas, K.; Ziagas, E.; Nakos, G. Geographical Map of Lands, Sheets: Perivolia, Vatolakkos. Scale 1:50,000; Ministry of Agriculture, Forest Protection Service: Athens, Greece, 1997. (In Greek) [Google Scholar]
  46. YPEN (Ministry of Environment and Energy). River Basin Management Plans: Approved Management Plans & Methodologies. European Union, NSRF 2014–2020. 2021. Available online: http://wfdver.ypeka.gr/el/home-gr/ (accessed on 20 July 2022). (In Greek)
  47. McDonald, M.G.; Harbaugh, A.W. A Modular Three Dimensional Finite-Difference Ground-Water Flow Model, US Geological Survey Open-File Report; USGS: Washington, DC, USA, 1988; pp. 83–875. [Google Scholar]
  48. Fest, R. Visual MODFLOW 2011.1 User’s Manual: For Professional Applications in Three-Dimensional Groundwater Flow and Contaminant Transport Modeling; Schlumberger Water Services: Kitchener, ON, Canada, 2011. [Google Scholar]
  49. COST Action 65. Hydrogeological Aspects of Groundwater Protection in Karstic Areas, Final Report (COST Action 65). In Report EUR 16547 EN; European Commission, Directorate-General XII Science, Research and Development: Brussels, Belgium, 1995; ISBN 92-827-4682-8. [Google Scholar]
  50. Lipfert, G.; Tolman, A.L.; Loiselle, M.C. Methodology of delineating wellhead protection zones in crystalline bedrock in Maine. J. Am. Water Resour. Assoc. 2004, 40, 999–1010. [Google Scholar] [CrossRef]
  51. Pochon, A.; Tripet, J.P.; Kozel, R.; Meylan, B.; Sinreich, M.; Zwahlen, F. Groundwater protection in fractured media: A vulnerability-based approach for delineating protection zones in Switzerland. Hydrogeol. J. 2008, 16, 1267–1281. [Google Scholar] [CrossRef] [Green Version]
  52. Różkowski, J. Evaluation of intrinsic vulnerability of an Upper Jurassic karst-fissured aquifer in the Jura Krakowska (southern Poland) to anthropogenic pollution using the DRASTIC method. Geol. Q. Pol. Geol. Inst. Natl. Res. Inst. 2007, 51, 17–26. [Google Scholar]
  53. TUC (Technical University Crete, Chania, Greece). Development of Exploitation Methods of the Underground Water in West Crete. PENED 87ED53: Program for the Researchers Support; Ministry of Industry, Research and Technology/General Secretariat for Research and Technology GGET, Progress Report. Unpublished work. 1990. (In Greek) [Google Scholar]
  54. Dimitriou, E.; Tsintza, P. Hydrogeologic investigations in Western Crete by using isotopic analyses and GIS techniques. J. Water Resour. Protect. 2015, 7, 923–937. [Google Scholar] [CrossRef] [Green Version]
  55. Marín, A.I.; Andreo, B.B. Delineating Source Protection Zones of Karst Springs. The Case Study of Villanueva del Rosario Spring (Southern Spain). In Advances in Research in Karst Media; Andreo, B., Carrasco, F., Durán, J.J., LaMoreaux, J.W., Eds.; Springer: Berlin/Heidelberg, Germany, 2010; p. 317. [Google Scholar] [CrossRef]
Figure 1. Springs catchment area.
Figure 1. Springs catchment area.
Water 15 00231 g001
Figure 2. Geological map of the study area. The rectangle shows the area of the detailed flow model.
Figure 2. Geological map of the study area. The rectangle shows the area of the detailed flow model.
Water 15 00231 g002
Figure 3. Relationship between the rainfall and the latitudes of the rain gauge stations.
Figure 3. Relationship between the rainfall and the latitudes of the rain gauge stations.
Water 15 00231 g003
Figure 4. Schematic geological sections. Their locations are defined in Figure 2. (A) NW-SE geological cross section through Omalos Plateau, (B) SW-NE geological cross section from Omalos to Meskla settlement, (C) SW-NE geological cross section across Theriso Gorge, (D) Cross section through Myloniana well-field and Agyia’s group springs.
Figure 4. Schematic geological sections. Their locations are defined in Figure 2. (A) NW-SE geological cross section through Omalos Plateau, (B) SW-NE geological cross section from Omalos to Meskla settlement, (C) SW-NE geological cross section across Theriso Gorge, (D) Cross section through Myloniana well-field and Agyia’s group springs.
Water 15 00231 g004
Figure 5. Geological map of the aquifer discharge zone (fine scale modeling area).
Figure 5. Geological map of the aquifer discharge zone (fine scale modeling area).
Water 15 00231 g005
Figure 6. Soil map of the study area. It was constructed based on data gathered from the geographical maps of the land [45], and the field survey.
Figure 6. Soil map of the study area. It was constructed based on data gathered from the geographical maps of the land [45], and the field survey.
Water 15 00231 g006
Figure 7. Vulnerability map generated through the EPIK (a), PRESK (b) and DRISTPI (c) method.
Figure 7. Vulnerability map generated through the EPIK (a), PRESK (b) and DRISTPI (c) method.
Water 15 00231 g007
Figure 8. Vulnerability class distribution.
Figure 8. Vulnerability class distribution.
Water 15 00231 g008
Figure 9. Model verification. Calculated vs. observed groundwater head in observation wells.
Figure 9. Model verification. Calculated vs. observed groundwater head in observation wells.
Water 15 00231 g009
Figure 10. Piezometric map in the modeling area. Additional data from [21] contributed to this design.
Figure 10. Piezometric map in the modeling area. Additional data from [21] contributed to this design.
Water 15 00231 g010
Figure 11. Groundwater flow path lines to the springs and pumping wells. The isochronal curves of 50-days are also shown.
Figure 11. Groundwater flow path lines to the springs and pumping wells. The isochronal curves of 50-days are also shown.
Water 15 00231 g011
Figure 12. Protection zones in the discharge area of Agyia’s karstic system.
Figure 12. Protection zones in the discharge area of Agyia’s karstic system.
Water 15 00231 g012
Table 1. Zones and classes of groundwater vulnerability defined by the EPIK and PRESK method.
Table 1. Zones and classes of groundwater vulnerability defined by the EPIK and PRESK method.
Zones of Groundwater VulnerabilityEPIKPRESK
CodeVulnerability ClassFp
Index
Precentage (%) of the Study AreaIg
Index
Precentage (%) of the Study Area
Τ5Very high<134>2.44
Τ4High13–17181.9–2.418
Τ3Moderate18–22451.3–1.855
Τ2Low23–27240.7–1.216
Τ1Very low>279<0.77
Table 2. Zones and classes of groundwater vulnerability defined by DRISTPI method.
Table 2. Zones and classes of groundwater vulnerability defined by DRISTPI method.
Zones of Groundwater VulnerabilityDRISTPI
CodeVulnerability ClassId
Index
Precentage (%) of the Study Area
Τ5Very high140–1909
Τ4High110–13944
Τ3Moderate80–10923
Τ2Low50–792
Τ1Very low19–4922
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Steiakakis, E.; Vavadakis, D.; Mourkakou, O. Groundwater Vulnerability and Delineation of Protection Zones in the Discharge Area of a Karstic Aquifer—Application in Agyia’s Karst System (Crete, Greece). Water 2023, 15, 231. https://doi.org/10.3390/w15020231

AMA Style

Steiakakis E, Vavadakis D, Mourkakou O. Groundwater Vulnerability and Delineation of Protection Zones in the Discharge Area of a Karstic Aquifer—Application in Agyia’s Karst System (Crete, Greece). Water. 2023; 15(2):231. https://doi.org/10.3390/w15020231

Chicago/Turabian Style

Steiakakis, Emmanouil, Dionysios Vavadakis, and Ourania Mourkakou. 2023. "Groundwater Vulnerability and Delineation of Protection Zones in the Discharge Area of a Karstic Aquifer—Application in Agyia’s Karst System (Crete, Greece)" Water 15, no. 2: 231. https://doi.org/10.3390/w15020231

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop