Intrinsic vulnerability of the Isonzo/Soča high plain aquifer (NE Italy – W Slovenia)

ABSTRACT The paper presents the map of intrinsic groundwater vulnerability of the Isonzo/Soča High Plain, which is located between the Collio Hills and the Classical Karst Region and holds an aquifer shared between Italy and Slovenia. The map, produced at a scale of 1:25,000 and printed in A0 format, was obtained by means of the SINTACS method and shows the intrinsic vulnerability of the aquifer in terms of seven vulnerability classes, from extremely high to low. It is accompanied by four supplementary sketches that illustrate the geological framework, the bedrock top surface, the groundwater flow paths, the Hazard Index map and three diagrams that summarize the percentages of vulnerability classes and of Hazard Index classes of the study area.


Introduction
The Isonzo/Soča High Plain is located on the eastern side of the Friuli Venezia Giulia Region, straddling the border between Italy and Slovenia. It holds a significant phreatic aquifer that represents an important natural wealth, in terms of quantity, quality and ease of supply (Zini et al., 2013). The aquifer is used for drinking, household, industrial, and agricultural purposes. The increasing interest of the local administrations, private water suppliers and the resident population using this important resource gave rise to the 'GEP' and 'ASTIS' Projects (Bisaglia et al., 2014;Zini et al., 2014a). As part of these projects, this research focuses on the evaluation of the intrinsic vulnerability of the Isonzo/Soča aquifer that can be used as a tool to safeguard the groundwater resource and support environmental protection and management policies. The intrinsic vulnerability of aquifers is defined as the specific susceptibility of aquifer systems to ingest and diffuse fluid or/ and hydro-vectored contaminants, whose impact on the ground water quality is a function of space and time (Civita, 1994(Civita, , 2010. The vulnerability of an underground water body depends on the hydro-lithology and hydro-structure of the hydrogeological system, the nature of soil and overburden, its recharge, groundwater inflow-outflow processes, the physical and hydro-geochemical processes that produce natural quality of water and the attenuation of the contaminants affecting the system (Civita & De Maio, 2000).
Moreover, this work aims to improve the existing researches (Cucchi et al., 1999;Cucchi, Massari, Oberti, & Piano, 2002) with new and updated data and with the integration with the Slovenian part of the aquifer.

Study area
The Isonzo/Soča is an alpine river 138 km long that runs from the Julian Alps (source: southern side of Mt. Travnik, Slovenia, 990 m a.s.l.) to the North Adriatic Sea (mouth: near Monfalcone, Italy) (Figure 1(a)). Except along the sector between Tolmin and Gorizia, where the river flows over carbonate rocks (Figure 2 (a)), the Isonzo/Soča riverbed consists of a thick fluvioglacial and alluvial cover (D'Ambrosi & Mosetti, 1972).
The study area extends for about 150 km 2 , taking in the southern stretch of Isonzo/Soča, from the outlet on the plain, north of Gorizia, to the aqueduct pumping wells at San Pier d'Isonzo (Figure 2(b)).
The Isonzo/Soča High Plain represents the easternmost part of the Friuli Plain, which hosts several aquifers composed of thick Eocene-to Quaternary-aged gravel and sand beds interbedded with clay and silty layers (Cucchi, Franceschini, & Zini, 2008 and references therein) and is bounded to the north by the Collio (ITA) and Vrtojba (SLO) hills (Figure 2(b)), made up of marlstones and sandstones belonging to the Eocene Flysch Formation (Venturini & Tunis, 1991). To the South it is bounded by the limestones and dolomites belonging to the Cretaceous-Paleogene carbonate platform of the Classical Karst Region (Jurkovšek et al., 2016). Quaternary sediments constitute the Isonzo/Soča and Torre megafans, large fan-shaped landforms that were formed during the Last Glacial Maximum as a result, respectively, of the drainage of the Isonzo/Soča (Bavec, Tulaczyk, Mahan, & Stock, 2004) and Tagliamento glaciers (Fontana et al., 2014). The megafans show a longitudinal sedimentary differentiation with a strong contrast between the proximal sector, where gravels prevail and the distal one, where fine-grained sediments dominate (Fontana, Mozzi, & Bondesan, 2008).
The Resurgence Belt divides the Isonzo/Soča Plain into the High and the Low Plain. The High Plain consists mainly of coarse and very permeable sediments, related to the above-mentioned Isonzo/Soča and Torre rivers, as well as other minor streams such as Judrio and Versa. These sediments hold a phreatic aquifer. Rivers dissipate a great amount of water during their way across the High Plain and for this reason the Torre and Judrio rivers remain dry for most of the year while the Isonzo/Soča loses about 26% of its discharge (Zini et al., 2013). These river losses, together with the rain and runoff waters coming from the hills, recharge the phreatic aquifer of the High Plain. Towards the Low Plain, the phreatic aquifer joins with a multilayered aquifer system characterized by alternating gravel-sand and clay-silt deposits (Stefanini & Cucchi, 1976). The gradual decrease in permeability southwards constrains the outflow of the High Plain phreatic waters in correspondence to the Resurgence Belt (Figure 1(b)).
The study area includes 23 municipalities (20 in Italy and 3 in Slovenia). On the Italian side, the groundwater is tapped and distributed by two main aqueducts (Irisacqua and AcegasApsAmga) that supply more than 300.000 citizens of the provinces of Gorizia and Trieste. On the Slovenian side the groundwater is used for domestic and agricultural purposes.
an unknown use, 15% are industrial, 6.8% are used for agricultural and irrigation purposes, 5.4% supply drinking water and aqueducts, 0.5% are used for fish breeding and 3% are geognostic boreholes while the remaining 14% are for various other applications. The overall withdrawal amounts to about 2 m 3 /s.

Materials and methods
The first step for the production of the intrinsic vulnerability map was the assessment of the climatic and hydrogeological settings.
The average yearly effective precipitation (P e ) was calculated by means of the equation: where Er is the yearly average effective evapotranspiration (Turc, 1954) and P the yearly average precipitation. A 30-years period  was considered for an accurate assessment of the P e , taking into account daily precipitation and temperature for several meteorological stations, covering an area of 5817 km 2 straddling Italy and Slovenia. Rainfall and temperature data were mapped as a function of the orography, using a co-kriging interpolation method based on the stations' elevations.
A monitoring survey was carried out to determine depth to the water table and to identify groundwater flow paths in different hydrological conditions (Figure 2(d)). Depth to the water table was measured manually, biweekly, for 2 years (2012)(2013)(2014) in 58 wells (46 in Italy, 12 in Slovenia). Several isophreatic maps for different groundwater flow conditions were drawn up using the Natural Neighbor interpolation tool.
Moreover, the depth, temperature and specific electric conductivity of the groundwater were continuously monitored by means of diver data loggers, installed in seven wells (three in Italy and four in Slovenia). Seven Slug and Bail Test (Bouwer & Rice, 1976;Bouwer, 1989) and two Aquifer Pumping Tests were performed (Figure 2(e)). In addition, hydraulic conductivity was defined using the Specific Capacity Method (Theis, 1963) in 24 piezometers.
A 3D model of the bedrock was designed starting from the stratigraphic columns of 243 wells, 55 of which reach the bedrock, together with 197 profiles of Vertical Electrical Soundings, 47 of which reach the bedrock. Three seismic sections and three electrical tomography profiles were also considered (Berlasso & Cucchi, 1991;Accaino, Boehm, Busetti, Baradello, & Affatato, 2015). The integrated analysis of all these data provided the bedrock surface model with a grid cell of 10 × 10 m. Afterwards, by subtracting the bedrock surface from the DTM (10 × 10 m), the isopach map of the alluvial sediments was obtained.
The intrinsic vulnerability was evaluated using a point count system model, the SINTACS method (Civita & De Maio, 1997;Civita & De Maio, 2000). It takes into account seven parameters: the depth to groundwater table 'S', effective infiltration 'I', the unsaturated zone attenuation capacity 'N', the soil attenuation capacity 'T', the hydrogeological characteristics of the aquifer 'A', the hydraulic conductivity 'C' and the topographical slope 'S'.
An increasing score is assigned to each parameter from 1 to 10, where the lowest represents conditions that guarantee greater protection for the groundwater and thus low vulnerability, while a maximum score means low or no protection and therefore high vulnerability. Scores are attributed to each cell into which the area is divided, using special charts and reference tables. Afterwards it is necessary to identify the different hydrogeological contexts and each situation is quantified using a weight string made up of seven integers with values ranging from 1 to 5. Each number is a weight coefficient multiplier that serves to enhance or mitigate the effect of each parameter in the specific scenario encountered (see Section 4.3).
The SINTACS index calculated for each cell is the sum of the weighted scores of the seven parameters, while the intrinsic vulnerability map is obtained by normalizing the SINTACS index to 100.

Base maps
The achievements of the climatic and hydrogeological analysis that constitute the base maps of the vulnerability evaluation were the representation of precipitation and temperature as a function of the elevation, the isophreatic trend map and the reconstruction of the bedrock surface.
In particular, the analysis of long-term meteorological data was carried out for daily precipitation collected at 15 stations and air temperature at 10 stations. We computed mean annual precipitation (MAP) and mean annual temperature (MAT) for each weather station during the years 1981-2010. MAP varied between 965 and 2572 mm/yr, while the minimum MAT was 4.6°C and the maximum 15.0°C. Temperatures were strongly affected by topography and stations located in the flat areas closer to the Adriatic Sea had higher MATs, while the stations in mountainous sites had lower ones. Similarly, precipitation was influenced by topography and stations in flat areas had low MAPs compared to mountainous locations. We therefore analyzed the relations between the measured MAP and MAT with topography and found significant simple linear regressions between the elevation of the weather stations and MAT (MAT = −0.006*elevation + 13.781; R 2 = 0.93, n = 10) as well as between elevation and MAP (MAP = 1.206*elevation + 1100.746; R 2 = 0.83, n = 15). This strong effect of elevation on both precipitation and temperature was considered in the kriging approach used to interpolate MAP and MAT on the study area. Interpolated maps were then used as a basis for the estimation of Er and P e for each grid cell.
The Isophreatic map was drawn up using the depth to water table under high flow conditions during the first two weeks of February 2014 at a time of one of the highest groundwater levels ever recorded (lowerleft side of the Main Map).
The isophreatic trend highlights a NE-SW groundwater flow, turning southwards in correspondence to the residential area of Sagrado. In detail, the contribution coming from the Isonzo/Soča and Vipacco/ Vipava rivers and the surrounding hills is clearly visible in the northeastern part, together with the input from the Torre River in the western sector that brings about the southwards turn of the groundwater flow. South of Gorizia the groundwaters are conveyed towards the Classical Karst reliefs, where they partially supply the karst aquifer (Zini, Calligaris, & Zavagno, 2014b). At the same time, these waters, together with the contribution of precipitation, supply the Isonzo/ Soča High Plain aquifer between Fogliano Redipuglia and Ronchi dei Legionari (D'Ambrosi, 1972).
Finally, the Bedrock isobath map highlights a complex buried morphology characterized by valleys and reliefs formed as a result of past erosion and tectonic movements. The map shows a depression south from Gorizia, where the bedrock was detected at more than 100 m below mean sea level (b.s.l.). Another significant depression starts at 50 m b.s.l., in the municipality of Mariano del Friuli, and gradually reaches 350 m b.s.l., in the surroundings of Villesse (left side of the Main Map).

SINTACS parameters
On the base of the above-mentioned layers, the seven thematic maps required by SINTACS method were elaborated.

Depth to groundwater table 'S'
This parameter plays a very important role in the vulnerability assessment since, together with the unsaturated zone characteristics, it defines the travel time of a hydro-vectored or fluid contaminant and the duration of the attenuation process of the unsaturated thickness and, in particular, the oxidation process by atmospheric O 2 (Civita, 2010). The depth to groundwater table map (Figure 3(a)) was obtained by subtracting the piezometric surface from the topographic height acquired from the DTM. In particular, the water table elaborated in the high groundwater flow conditions recorded during the first weeks of February 2014, was used. Depth to groundwater values varied from few meters, in correspondence to streams and the most depressed areas, to more than 35 meters below the city of Gorizia. Consequently, the 'S' scores vary from 1 to 10 (Figure 3(b)).

Effective infiltration 'I'
The effective infiltration directly affects the dragging down of pollutant and its subsequent dilution during its passage through the unsaturated zone and then within the saturated zone (Civita & De Maio, 2000). The parameter values are calculated starting from the effective precipitation (P e ) coupled with the potential infiltration index (χ), which is a function of the surficial lithology (in the case of outcropping or under conditions of scarce soil cover) or of the hydraulic characteristics of the soil (where its thickness is more than 0.5 m).
The effective infiltration (I) for each cell was obtained by multiplying P e with χ. The results vary from 40 to over 700 mm/year, as a function of the soil permeability. The 'I' scores thus vary, from 2 to 10 (Figure 3(d)).

Unsaturated zone attenuation capacity 'N'
The unsaturated zone is composed by the subsoil included between the lower surface of the soil and the saturated zone. Its attenuation capacity is assessed starting from the hydro-lithological features such as texture, mineral composition, grain size, fracturing and karst development (Civita & De Maio, 2000).
In order to evaluate the 'N' ratings for the unsaturated zone, a score was assigned to the different hydro-lithological levels located between the ground and the water table for the 210 stratigraphic columns available (Table 1) and then a weighted average was calculated, based on their thickness. The weighted values varied from 5 to 9; in particular the highest values are observed where the unsaturated zone is richer in gravel and in correspondence of the riverbeds (Figure 4(a)).

Soil attenuation capacity 'T'
The type of overburden and particularly of the soil plays a significant role in the attenuation process for a contaminant travelling within a hydrogeological system. Indeed, within the soil, fundamental processes take place and provide some degree of attenuation (Civita & De Maio, 2000).
The physical (e.g. grain size, texture and bulk density) and chemical characteristics of the soil (e.g. pH, cationic exchange capacity and organic matter content) affect the soil's attenuation capacity. In order to calculate the 'T' scores, the soil (Michelutti, Barbieri, & Bianco, 2006;Stepančič et al., 1992) and the land use (RAFVG, 2003; TTN5) maps were categorized as a function of the textures (Figure 3(e)). Values of 'T' vary from 2, where clayey and/or silty soils outcrop and at the base of reliefs, to 10, where the soil is absent, such as in urban areas, and where clean gravel and sand prevail or constitute the riverbeds (Figure 3(f)).

Hydrogeological characteristics of the aquifer 'A'
The hydrogeological characteristics of the aquifer describe the processes that involve a contaminant below the piezometric level when it mixes with groundwater and loses part of its original concentration. These processes are molecular and kinematic dispersion, dilution, sorption and chemical reactions (Civita & De Maio, 2000).
Considering only the levels located below the water table of the 210 available stratigraphic columns and weighting them as a function of their thicknesses, the hydrogeological characteristics of the aquifer were defined and 'A' scores were assigned ( Table 1). The 'A' scores lie between 6 and 9; the highest value is correlated to the presence of coarse alluvial deposits, especially along the riverbed of the Isonzo/Soča (Figure 4(b)).

Hydraulic conductivity 'C'
Hydraulic conductivity represents the groundwater mobility capacity within the saturated aquifer and thus the mobility potential of a contaminant having density and viscosity almost the same as the groundwater itself (Civita & De Maio, 2000).
Hydraulic conductivity was evaluated for 33 wells through the above-mentioned hydrogeological tests and spatially interpolated to obtain a map of hydraulic conductivity of the whole study area. Thereafter, the map was critically reviewed to take into account the hydrogeological characteristics of the aquifer.
In general, the hydraulic conductivity of the Isonzo/ Soča High Plain is high, the values lying between 10 −2 and 10 −4 m/s and thus the 'C' scores for the alluvial area vary from 7 to 10 (Figure 4(c)).

Topographical slope 'S'
The topographic slope governs the amount of surface runoff and the velocity of water (a fluid and/or hydro-vectored contaminant) over the equipotential surface (Civita & De Maio, 2000).
The slope map was drawn up from the DTM. Slope values in the plain are generally low and, consequently, almost everywhere, the 'S' scores reach the maximum value of 10 ( Figure 4(d)).

Hydrogeological and impact situations
Four different hydrogeological and impact situations were identified: 'severe impact', 'seepage', 'fissured' and 'karstified' (Figure 5). The most widespread situation in the plain is that of 'severe impact' and is represented by scenarios of low-slope gradient and intense urbanization, diffuse agriculture and/or industrialization.  1. Scores of "N" (Unsaturated zone attenuation capacity) and "A" (Hydrogeological characteristics of the aquifer) assigned to the different hydro-lithological levels.

Hydro-lithological level
N score A score Clay 1 1 Clay and/or silt 1.5 2 Clay and/or silt with sand or gravel 2 3 Clay with sand or gravel 2.5 4 Silt with sand or gravel 3 5 Sand with clay or silt 4 6 Soil 4 6 Silty sand 4.5 6.5 Flysch 5 4 Sand with clay and gravel 5 7 Sand with gravel and silt 5.5 7.5 Partially cemented gravel with clay 6 6 Gravel with clay and silt 6 8 Partially cemented gravel with sand 6.5 6.5 Sand 6.5 8.5 Partially cemented gravel 7 7.5 Gravel with clay and sand 7 8 Gravel with sand and clay 7.5 8 Gravel and sand 8 8.5 Gravel 8.5 9 Gravel with pebbles 9 9 Limestone 9 9 Along the riverbeds delimited seepage areas occur, within which exist or could exist considerable water exchange between surface water bodies and groundwater. The outcropping rocks are classified either as 'fissured' (Flysch) or as 'karstified' (limestone). The weight strings used are shown in Table 2.

Hazard sources and hazard index map
An evaluation procedure of hazard sources, promoted by the European Union (Various authors, 2004), points out that a range of factors contribute, either individually or together, to render any human activity a 'hazard source'. Indeed, sources of pollution are associated with a wide range of industrial, agricultural, commercial and domestic activities. Hazard sources can be distinguished as punctiform (such as a petrol station) or diffuse (such as a cultivated area).
To integrate the intrinsic vulnerability map, the pollution sources were identified and digitized. The hazard of each point or area was evaluated using the rating system proposed by Civita and Zavatti (2006) for the anthropic and/or industrial activities and integrated with the hazard factors suggested by Trevisan, Padovani, Errera, Capri, and Del Re (1998) for the agricultural areas. Initially the types of hazard were mapped, separating the two classification systems. Figure 6 shows the potential hazard sources, according to Civita and Zavatti (2006). This classification associates a Hazard Index (HI) to each activity, as function of different hazard factors (e.g. the production of special or hazardous waste, the discharge of organic or inorganic pollutants). Punctiform and areal data were collected from different sources, such as the Chamber of Commerce of Gorizia for industrial and  manufacturing plants, the Regional Agency for Environmental Protection (ARPA) for liquid storage, AcegasApsAmga and Irisacqua for urban wastewater and municipal waste management and the S.I.AGRI.FVG (RAFVG, 2006) database for the location of animal husbandry activities. In order to cover all the categories, the data regarding traffic and transport, recreational facilities and quarrying were extracted from an online available database (Open Street Map), from the land use map (RAFVG, 2003), from the layers of the Regional Technical Map and from the Regional Geological-Technical Map (RAFVG, 2008). Most of the data were finally verified through orthophoto assessment or the Street View tool (Google Maps). Almost the 70% of the study area has no classified hazard, according to the categories proposed by Civita and Zavatti (2006). The main source of pollution hazard is derived from the urban wastewater and municipal solid waste management systems (26.31%, 8 ≤ HI ≤ 15) with all the other categories covering the remaining 3.88%. Figure 7 reports the main land use classes according to Trevisan et al. (1998). This system correlates a Hazard Factor (HF) value to each class, based on the probability that the pollution is caused by pesticides and/or fertilizers. Three RapidEye satellite images of 2012 were acquired to classify crops, which have different phenological phases in the course of the year. Specifically, 1st April (to detect the presence of vineyards, wheat and barley and grassland), 17th June (to classify maize, soybeans, meadows and vineyards) and 24th of October (to identify meadows and vineyards) were chosen. Contextual image classification methods using sequential maximum a posteriori (SMAP) estimation was used in GRASS 6.4 GIS. Due to the highly fragmented nature of the territory and the shift in crop cultivation (mainly alternation between barley and soya), each image was classified separately and then combined using a complex postclassification processing. Data belonging to 'Urbanization' and 'Traffic and transport' categories were derived from the existing land use maps (RAFVG, 2003;TTN5). More than 38% of the Isonzo/Soča High Plain is affected by shifting cultivations (HF = 5), more than 17% is occupied by urban fabric (HF = 1-2), about 11% by vineyard (HF = 9) and about 15% is covered by other vegetation (HF = 0-1) with the other categories representing less than 5% of the Figure 6. Hazard sources according to Civita and Zavatti (2006). total. In this case, the HI is obtained by adding to the HF other four control factors related to agronomic practices, climatic conditions, irrigations and slope (Trevisan et al., 1998).
In order to calculate a single HI, the two maps ( Figures 6 and 7) were overlapped and the worse condition, i.e. the higher HI, was taken as prevailing for creating a unique HI map (lower-right side of the Main Map).

Conclusions
The SINTACS intrinsic vulnerability map shows that the Isonzo/Soča High Plain aquifer is fairly vulnerable as a result of its hydrogeological characteristics, the most exposed areas occurring in correspondence to the riverbeds and, in the south, at the line of contact between High and Low Plain, along the Resurgence Belt.
Almost half of the study area (49%) has high vulnerability. 37% of the area is classified as of very high vulnerability, while 11% of territory has extremely high values. Only 2% and 1% of the study area are considered as having medium and low vulnerability, respectively. There is no territory with very low vulnerability. In general, the vulnerability is greater in the flat areas characterized by very high hydraulic conductivity and minimum depths to groundwater.
Conversely, the HI map indicates that the surface water bodies and their flat surroundings have minimum HI values. Medium values are identified around almost all of the towns. The higher HI values correspond to the industrial zones and to vineyards or shifting cultivation. However, classes of HI equal or greater than 11 constitute about the 20% of the study area. By combining the analysis of the Vulnerability and HI maps, the most vulnerable areas correspond to low values of HI (≤10) and in less than 3% of areas extremely high vulnerability coincide with high HI values (>10).
The Intrinsic vulnerability map highlights that the aqueduct wells of Gorizia and Farra d'Isonzo are located in extremely vulnerable areas. The wells located in the surroundings of San Pier d'Isonzo, that guarantee the supply of water to the city of Trieste (located about 35 km to the SE), are set within a slightly better context (high vulnerability), given the presence of finer and less permeable sediments.  Trevisan et al. (1998).

Software
The maps were produced by means of Esri ArcGIS v. 10.1 and GRASS 6.4 GIS while Adobe Illustrator CS6 and Adobe InDesign CS6 were used for graphic operations and layout building respectively.