Delineation of Protection Zones for the Main Discharge Area of the Gran Sasso Aquifer (Central Italy) through an Integrated Geomorphological and Chronological Approach

The paper proposes a study for the delineation of protection zones in the main discharge area of the Gran Sasso aquifer (Central Italy). At this aim, starting from a detailed geological and hydrogeological reconstruction, the study was divided into the following phases: 1) development of a conceptual model of water flow in the study area; 2) creation of a 3D numerical model in order to simulate the groundwater flow in saturated conditions, both at the basin and at fine-scale; 3) flow path analysis through deterministic and stochastic approaches; 4) assessment of the aquifer vulnerability based on a geomorphological analysis of the catchment area. Conceptual and numerical models were then used to delineate protection zones for wells and springs with chronological criterion and geomorphological-structural criterion (based on the EPIK method). The results show that with a chronological approach protection zones are located along the main flow directions, corresponding to the areas surrounding wells and springs with high hydraulic conductivity values (faults and thrusts) within the satured zone. On the contrary, the geomorphological method has found some important protection zones also quite far from wells and springs, in areas characterized by quick infiltration processes. The protection zones delineated with the stochastic method were finally intersected by the vulnerability map obtained with the EPIK method, to take into account both filtration and infiltration processes. The results show the local vulnerability of the groundwater to the pollution, with protection zones extending 1 to 5 km towards northeast from springs and wells.


Introduction
The problem of the supply and the protection of groundwater resources is currently a highly topical subject. In particular, karst aquifers often constitute important and strategic water resources in many mountain areas. Yet, these aquifers are highly vulnerable to degradation and pollution problems: geology (fractured carbonate rocks), morphology (presence of epikarst and network of karst cavities), and hydrogeology (rapid concentrated flow through fractures and conduits) strongly favor the movement of contaminants towards the water table and then to springs or wells [1] [2]. Consequently, several methods were developed to assess the vulnerability of karst groundwater to pollution and to define proper protection zones [3] [4], as required by law (i.e., [5]- [7]). In accordance with EU Directives, Italian laws [8] require a protection zoning system for groundwater sources, but its application is generally based on fixed-radio methods, especially for heterogeneous aquifers, whereas the use of more complex methods, such as isochrones methods and modeling approaches, is so far very limited and generally, it concerns wellhead protection zones in alluvial aquifers.
The overall objective of the present study has been to delineate protection zones for the main discharge area of the Gran Sasso aquifer (Figure 1). This wide fractured-karst aquifer has an average flow rate of about 25 m 3 /s, and therefore it is one of the most important water resources of Central Italy. In particular, the present report aims to update the prior knowledge about the main discharge area of the aquifer, where more than 75% of the groundwater flowing in the aquifer rises to the surface supplying springs and a large public-water service. The interest for this area is related not only to its great availability of groundwater resources but also to contamination issues. The aquifer has been the center of attention in recent years due to a local decrease in its qualitative characteristics because of the presence of a heavy industrial pollution, which resulted in the sudden closure of important public-water wells. New public water wells were then realized just a few kilometers upstream from the polluted site. Nowadays it is obviously a prime socio-economic issue to guarantee the protection of the water resource against further impact. To this aim, a fine-scale hydrogeological conceptual model was reconstructed and then used to develop a groundwater flow model and to delineate proper protection zones for wells and springs in the studied area with different methods: 1) a chronological method, based on flow simulation in saturated medium with both a deterministic and stochastic approach; 2) a vulnerability-based approach (EPIK method); 3) a combination of the two previous methods, in order to take into account both the filtration and infiltration processes.

Previous Knowledge on the Hydrogeological Setting at Basin Scale
The Gran Sasso aquifer (Central Italy, Figure 1(a)) can be identified with a carbonate karst hydrogeological system having an area of about 800 km 2 , and an altitude ranging between 2912 and 270 m a.s.l.. From a lithological point of view, the hydrogeological system of the Gran Sasso massif consists of carbonate rocks of Meso-Cenozoic units belonging to slope-to-basin lithofacies [9]. These carbonate rocks have a medium-high hydraulic conductivity, which determines the presence within the massif of a large regional aquifer having a flow rate of about 25 m 3 /s. The recharge of the aquifer is on average equal to 500 mm/year [10] and it arises from meteoric precipitations, mainly located in the northern tectonic-karstic area of Campo Imperatore (Figure 1(b)).
In the last decades, the Gran Sasso aquifer was the objective of a number of hydrogeological studies. Reference [11] was the first who outlined the geological and hydrogeological setting of the aquifer mainly based on the data arising from the construction of the motorway tunnels in the northern zone of the catchment area. Afterwards, several authors improved the knowledge on the Gran Sasso aquifer: references [12] and [13] defined the hydrodynamic scheme of the aquifer at catchment scale, based on hydrological, hydro chemical and hysotopic data. Reference [14] improved the hydrodynamic conceptual model, by analyzing the hydro chemical monitoring data on springs. Moreover, the availability of subsurface data arising from the construction of the motorway tunnels and scientific laboratories in the northern zone of the catchment area has favored the development of fine-scale studies in this zone [15]. On the contrary, there are few fine-scale studies for the southern zone of the aquifer. Among these studies, [16] proposed a detailed study of the Tirino River springs, achieving a hypothesis of the water table map of the zone, based on surface data. Previous studies at catchment scale pointed out that the groundwater flows from the central zone of the massif towards the marginal zones, and the structural setting, karst features, and depth of the dolomite bedrock control it. The tectonic setting of the massif also determines the permeability boundaries of the Gran Sasso aquifer system (Figure 1(b)). In particular:  [12]). The black rectangle in Figure 1(b) represents the study area of Figure 2.
• To the North and to the East: the main overthrust having E-W and therefore N-S direction, dipping respectively towards South and West, causes the tectonic overlapping of the carbonate units upon the clastic ones (Flysch della Laga), which constitute the regional aquiclude (no flow boundary); the main springs are located along this boundary, and they have an average total discharge equal to 1.6 m 3 /s along the northern boundary and equal to about 15 m 3 /s along the eastern boundary (Tirino Valley); • To SE: the discordant stratigraphic limit (locally tectonic) of the Quaternary units upon the carbonate ones, holding the main aquifer, determines a permeability boundary; • To SW: groundwater exchanges with the secondary hydrogeological system of the Sirente Mount.

Field Surveys and Investigations
In the present study the hydrogeological characterization at fine-scale was carried out in the main discharge zone of the Gran Sasso aquifer, located in the south-eastern part of the catchment area (Figure 1 and Figure 2), near the contact between two regional tectonic units, the Gran Sasso Unit and the Morrone Unit [17], [18]. In partic- Figure 2. Fine-scale geological map of the study area. Geognostic investigations and spring location are also shown. ular, the Gran Sasso Unit is overthrusted upon the Morrone Unit along a fault plane dipping from 10˚ to 60˚, named "Gran Sasso-Bussi Overthrust". The stratigraphic sequence in the study area consists mainly of carbonate rocks of the Gran Sasso Unit, ranging from Upper Jurassic to Paleogene. Quaternary deposits are present, which are mainly sands and pebbles, breccia, siltstone and travertine from Lower Pleistocene to the Present age. Carbonate rocks are characterized by widespread karstic landforms (e.g. surface micro-features and karrens in advanced stage of development).
A field survey and investigations were carried out in an area of about 20 km 2 ( Figure 2) and they included: • A detailed geological, hydrogeological and geostructural survey; • Subsurface surveys through the execution of: geophysical surveys (for a total length of 4.5 km, with both geoelectric and geoseismic measures), 11 drillings (depth ranging from 40 to 485 m, for a total length of more than 1.5 km), piezometric measures, and in bore-hole permeability tests; • The geochemical characterization and monitoring of springs and groundwater sampled from piezometers, with the definition of the standard chemical-physical parameters, electric conductivity, temperature and pH, as well as hysotopic analysis of 18 O/ 16 O and 2 H/H. The new data allowed improving the hydrogeological conceptual model of the studied area by: • The characterisation of the degree of fracturing of the rock masses and mapping the main karst and tectonic features, that represent preferential recharge and flow zones for groundwater; • The verification of the geochemical, hysotopic and hydrodynamic characteristics of groundwater in the studied area; • The reconstruction of the local trend of the impermeable bedrock, that constitutes the regional aquiclude of the whole aquifer; • The drawing of a water table map based on piezometric and springs data.
The following paragraphs synthetically describe the results of the geognostic investigations, and delineate a fine-scale hydrogeological conceptual model for the main discharge zone of the aquifer.

Hydrogeological Conceptual Model
The surface surveys and subsoil investigations provided a great number of data and information that were integrated and synthesized to achieve a fine-scale conceptual model.
From a geological and structural point of view, the interpretation of the data acquired from the geophysical surveys has been integrated with the surface geological survey to bring about the fine-scale geo-structural setting of the studied area. According to this scheme, a number of faults cause the displacement of the carbonate sequence, which is interrupted in depth by an important discontinuity: the main overthrust leads to the superimposition of the Gran Sasso tectonic unit over the Morrone tectonic unit. The latter consists in an impermeable bedrock, which acts as an aquiclude, determining a local uplift of the water table. Therefore, it favors the San Rocco Wells pumping (Figure 3, see cross-section A-A'), as well as the surfacing of groundwater along the Canestro Valley, where it supplies the alluvial aquifer and then the San Calisto springs (Figure 3, see cross-sec-tion B-B').
Data arising from drillings showed that the aquiclude has a structural top at an altitude of 475 m a.s.l. (Figure  3, see cross-section B-B'), and it therefore acts as a local aquifer divide (Figure 4): it divides the groundwater supplying the Tirino Valley (at east) from the groundwater supplying the Canestro Valley (at west). Based on this conceptual model, the San Calisto springs are supplied by the groundwater flowing in the western part of the catchment area, drawn by the tectonic structure of the Mt. Ospedalera, which is characterized by folds and faults. The San Rocco Wells are instead supplied by the groundwater flowing in the eastern part of the basin. After it passes the Parata Valley, it diverts its path from NW-SE to W-E, and it is stopped in its flow by the overthrust, which acts as an impermeable boundary.
Drilling and piezometers also gave data useful for drawing a water table map. Results of the available piezometric surveys confirmed a regional groundwater flow direction from NW towards SE, with high hydraulic gradients in the discharge zone, ranging from 3% and 8% (Figure 4).
As far as the hydrogeological balance is concerned, the new isotopic and chemical data confirmed that the water resources of the study zone are supplied by a wide catchment area located at an altitude much higher than the discharge zone: • The isotopic ratio of groundwater in the study area (altitude ranging from 230 and 700 m a.s.l.) corresponds to a recharge altitude higher than 1300 m s.l.m.; • Groundwater has the same calcium-magnesium-bicarbonate hydro chemical facies everywhere, with a very similar ionic content and medium-high electric conductivity values. Based on the groundwater balance of the whole Gran Sasso aquifer, the local contribution of the recharge in the studied area is quite small. Actually, the average infiltration in this area is equal to 150 mm/y [10], and it corresponds to a recharge less than 1.25 m 3 /s (equal to about 5% of the recharge in the whole basin). Therefore, the largest amount of groundwater (about 14 m 3 /s) inflows from the northern and western boundaries.
On the strength of geo-structural surveys, the Elementary Representative Volume (ERV [19]) ranges from 10 3 m 3 in the fractured superficial rock mass to a maximum of 10 6 m 3 in depth. Considering that the studied area is equal to about 20 km 2 and the aquifer has a thickness higher than 200 m, the ERV is significantly smaller than the domain size and the work scale can be considered a "very large field" [20]. The water flow occurs inside a fractured-karst medium that can be assimilated to a continuum. Moreover, at this work scale the fracture/conduit systems are uniform and well interconnected.
As a result, the hydraulic conductivity tensors and the corresponding equivalent values were calculated for the surface rock mass. Results show equivalent hydraulic conductivities ranging from 3 × 10 −3 to 6 × 10 −3 m/s, with maximum value of 10 −2 m/s and minimum value of 10 −4 m/s. The tensors have the major axis orientation parallel to the main tectonic structures. Based on the field surveys, along the fractures a local increase of about one order of magnitude in the hydraulic conductivity was assumed for a 50 m wide strip. In accordance with borehole permeability tests, the hydraulic conductivity generally decreases with increasing depth. The decrease is equal to about one order of magnitude in 100 m.

Groundwater Modeling
The groundwater flow model of the study area was developed in steady state by using the numerical model MODFLOW-2000, a modular finite difference groundwater model code [21]. Actually, on the strength of the previously discussed conceptual model, for a large flow scenario such as the present study, the treatment of a fractured-karst medium as homogeneous might lead to a reasonably good match with an equivalent porous medium.

Large-Scale Flow Model
The groundwater flow modeling was at first carried out on a large scale, by considering the whole catchment area of the Gran Sasso aquifer. The domain having a dimension equal to 50 × 50 km 2 was split into 115 × 117 cells (with size ranging from 200 m to 500 m), and the following boundary conditions were applied (Figure 5(a)): River and drain boundary conditions were defined in term of location and flow rate, according to data arising from previous studies (i.e., [12]). The only water inflow arises from the recharge, which according to previous studies ranges from 1000 mm/y to 100 mm/y [10] depending on the altitude and on the degree of karsification in the vadose zone (Figure 5(b)).
As far as the hydraulic conductivity is concerned, the regional aquifer was simulated as a single layer having a thickness ranging from 2000 m (in the northern recharge area) to 250 m (in the southern discharge zone). The hydraulic conductivity was considered generally uniform and equal to 6 × 10 −5 m/s with higher values (equal to 10 −3 m/s) just along the main tectonic features (see Figure 1(b) for their location). The values of the hydraulic conductivity were pointed out during the calibration process, which was carried out considering as a target the only available data at this scale of work: the flow rate of springs (Figure 6). The results of the large-scale flow model (Figure 7) confirm the flow path at basin scale, having a main flow direction from NW towards SE and a total discharge equal to 23 m 3 /s. These results were then used as input for the groundwater flow model at finescale.

Detailed Flow Model for the Aquifer Discharge Area
The fine-scale groundwater flow model involved the main discharge zone of the aquifer (Figure 7), with a domain having an extension equal to 6 × 10 km 2 . The modeling domain was split in 203 × 120 cells having a size ranging from 25 m to 50 m. On the strength of the detailed geological data available for this area (Figures 2-4),   In addition, constant head conditions (ranging from 565 m a.s.l. to 300 m a.s.l.) were applied along the northern, southern and western boundaries of the domain, based on the results of the large-scale flow model.
As far as the areal distribution of hydraulic conductivity is concerned, both the fractures and the presence of some diffuse karsification zones were considered. The surface hydraulic conductivity was assessed on the base of the geological-structural data arising from structural surveys, whereas the in-depth values were pointed out by the borehole tests ( Table 1). The values of the hydraulic conductivity along faults were later on calibrated, con-sidering the observed hydraulic heads measured in the piezometers. The results allowed reconstructing the water table map near the main discharge points and they were then used to delineate their protection zones.

Background
The protection and management of water resources is an issue of growing interest. In particular, karst aquifers  appear to be highly vulnerable to pollution, precisely because of their specific structure. For this reason, regulations in order to preserve water resources for human consumption have been introduced all over the world [5]- [7]. One of the main instruments introduced to prevent pollution phenomena is in the delineation of protection zones. Several methods for delineating protection zones in fractured and karst aquifers have been proposed [4] [22]- [27], but it still remains a challenge due to heterogeneity and anisotropy which often characterize fractured and karst aquifers. As a result, in these cases methodologies for delineating groundwater protection zones must generally be adapted to different site characteristics and to the type of data available.
In the present study, the following methods were used: • The isochrones method, with both a deterministic and stochastic approach, through the use of the previously developed flow model; • A geomorphological method, based on the fine-scale conceptual model of the study area.
Following sections describe the application of these methods, compare their results, and integrate them in order to find a more comprehensive map of protection zones.

Isochrones Method
The isochrones method is generally applied to the case of vulnerable springs and wells linked to slightly heterogeneous aquifers, as the case study is.
The previously described groundwater flow model (Figure 8) was then used to reconstruct the flow path and the particle tracking around wells and springs. Results (Figure 9) highlight, as was expected, a main flow path developed along tectonic features. According to Italian laws, the 60-day isochrones corresponds to the inner protection zone, where a number of activities are forbidden and strict safeguards for all surface water inflow points must be provided with a total prohibition of any effluent or spillage of polluting matter. The 365-day isochrones corresponds to the outer protection zone, where a number of restrictions for civil, productive, recrea-tional and agricultural land-uses are given. The inner protection zones obtained with this approach are only lo-cated around discharge points and they reach a maximum extension for the wells equal to about one kilometre, whereas the outer protection zone extends to about 3 km along the main flow direction.
It came out that the protection zones were controlled by the presence of tectonic features, and considering the great variability in the hydraulic conductivity values of these elements, the isochrones method was also applied by using a stochastic approach. To this aim, the hydraulic conductivity of fault zones was described as a log-normal distribution having the same average values used in the deterministic modeling ( Table 1) and a standard deviation equal to 1. In this way, the capture probability for the different time interval was pointed out (Figure 10).
As it is evident from the obtained results, the use of a stochastic approach characterised by a high dispersion of data brings about a significant enlargement of the capture zone, especially along the main fracture zones, which may constitute the main flow direction of the aquifer.

Geomorphological Method
Geomorphological methods take the protective effect of the unsaturated zone into account and assign less signi- ficance to transit times, which are included qualitatively in the parameter evaluation, but are not directly evaluated. These methods are very useful in highly heterogeneous media such as karst or highly permeable fractured aquifers, where the delineation of protection zones based exclusively on transit time in the saturated zone would result in the determination of excessively restrictive protection zones over large distances.
In the present study, the EPIK method [4] was used. It takes into account four factors: epikarst development (E), protective cover (P), infiltration conditions (I) and karst network development (K). Each factor is given a ranking index, according to the results of the field surveys. In particular, the P factor mainly includes the type and thickness of the sedimentary cover overlaying the karst rock. The weight coefficient is then attributed to each of the indexed factors according to their degree of protection. The final protection coefficient is therefore calculated and finally the protection zones were delineated (Figure 11, Table 2).
As it is quite evident, the most restrictive protection zones are located in those areas where infiltration processes are favored by the fracturing and karstic degree of outcropping rocks, and not necessarily near springs and wells. Figure 11. Protection zones delineated according to the EPIK method. Table 2. Groundwater protection zones according to the EPIK method [28].

Zone
Definition Characteristics S1 Wellhead protection zone It must prevent damage to the groundwater supply or artificial recharge facilities, as well as prevent pollution in its immediate surroundings S2 Inner protection zone It defines an area suitable for preventing pollution (i.e., germs, viruses, and contaminants arising from surface excavation or subsurface works) reaching the groundwater supplies S3 Outer protection zone It must provide sufficient space and time for remediation when accidental pollution threatens a groundwater supply

Discussion of Results
In the present case study, the main disadvantage in the use of the isochrones method for the delineation of protection zones is related to the significant depth of the groundwater table with respect to the land surface, which in some areas is higher than 200 m. Therefore, since the isochrones method is based on groundwater flow mod-elling within saturated soil, and it does not consider infiltration processes in the not saturated zone, this method could lead to an overestimation of protection zones. On the other hand, in the numerical flow model the recharge is applied to the highest active layer, which often does not correspond to the top layer where the main tectonic and karst features controlling infiltration are present. Therefore, the isochrones method could even underesti-mate the vulnerability of some areas, especially if these areas are not close to springs and wells.
To overcome these limitations, results of the two methods were finally integrated (Figure 12). At this aim, the 365-day probabilistic simulation was considered and the following criteria were chosen ( Table 3): • the immediate protection zone corresponds to areas where the capture probability is higher than 70% and contaminants can reach the groundwater (S 1 and S 2 of EPIK method); • the inner protection zone corresponds to areas where the capture probability is higher than 70% even if a good protection of the aquifer is guaranteed (S 3 zone of EPIK method) or areas having quite a lower capture probability (in between 40% and 70%) but with a high vulnerability of the aquifer (EPIK zone S 1 and S 2 ); • the outer protection zone corresponds to areas where the 60-day capture probability is lower than 40% or areas having where the capture probability is in between 40% and 70% by the aquifer is locally protected (S 3 of EPIK method).

Figure 12.
Protection zones delineated by integrating the isochrones methods with the geomophological one. Table 3. Criteria applied to integrate isochrones and geomorphological methods.

Conclusions
On the strength of a new fine-scale hydrogeological conceptual model of the discharge zone of the Gran Sasso aquifer, this study allowed the delineation of protection zones for one of the most important carbonate aquifers in Italy. In particular, the study demonstrates the necessity of combining detailed geomorphological analysis with groundwater modeling to delineate proper protection zones for springs and wells in the study area. Different methods were actually developed to delineate protection zones, both using groundwater flow mod-eling and geomorphological analysis. Different approaches obviously gave very different results in term of loca-tion and dimension of protection zones. These great differences arise from the assumptions of the applied me-thods, which do not seem to be completely suitable for the specific case study. To overcome these limits, an in-tegration between isochrones and geomorphological method is then pointed out allowing delineating protection zones, which take into account both filtration in saturated rock mass and infiltration processes in vadose zone.