Hydrogeological and Hydrochemical Regime Evaluation in Flamouria Basin in Edessa (Northern Greece)

Abstract

Groundwater quality deterioration and overexploitation constitute two critical environmental issues worldwide. In this study, with the aim to achieve a groundwater sustainability purpose, a preliminary hydrogeochemical survey is conducted in the Flamouria basin, Pella prefecture, Northern Greece using available and collected data. For this purpose, chemical analyses of groundwater, springs, and surface water were collected and analyzed with three electrical resistivity tomographies (ERTs). A Groundwater Quality Index (GQI), along with a nitrate susceptibility assessment is applied within the porous aquifer. The water quality analysis along with GQI application showed excellent water quality for potable and irrigation use however highlighted future issue for irrigation utilization as the high alkalinity and total dissolved solid (TDS)could generate excessive soil salinization. Moreover, the application of a methodology for the identification of “Nitrate Vulnerable Zone” called the Protection from Natural and Anthropogenic sources (PNA) highlighted the natural susceptibility to nitrate pollution of the porous aquifer, especially in the central part of the area where most agricultural activity is localized. The work further confirmed how the proposed elaboration could represent an easy and widely applicable hydrological assessment where there is also limited data available.

1. Introduction

The growing population, matched with increasing industrialization and intensification of the agricultural sector, has led to the worldwide use of many chemical and organic fertilizers to satisfy the rising needs of food and products. This phenomenon is becoming more evident within the last century due to agricultural development [1]. However, these agricultural practices were often responsible for the qualitative and quantitative deterioration of the local ground and surface waters resources as they can introduce several pollutants inside an aquifer system, consequently reducing resource availability. Doubtless, anthropogenic activities, such as agricultural or industrial practices, have gradually introduced many new pollutants inside regional aquifers or, in several cases, they are responsible of an increased concentration of existing species (i.e., NO₃, F, and SO₄) over the threshold limits [2,3]. Among all chemical species that have been introduced in groundwater, nitrate (NO₃) is considered as the main pollutant worldwide [4] and consequently, the monitoring of its release and spreading’s processes has become a primary goal of global hydrogeological research [5]. Many studies on groundwater pollution have been conducted in many of the world’s regions, highlighting the global concern for this phenomenon, showing an emphasis on how the increasing concentrations of NO₃ in groundwater could also generate severe side effects on human health [6]. This issue acquires greater importance especially within large agricultural areas, where the continuous and sometimes irrational use of fertilizers and manure practices can slightly increase NO₃ concentration in groundwater, pushing it over the acceptable threshold limit [7,8,9,10]. Thus, accurate and proper identification of the various pollution’s sources is becoming mandatory to reduce and avoid groundwater quality deterioration so as to finally achieve the utilization of sustainable resources. Moreover, the concept of groundwater quality is strictly connected with its final utilization and according to its hydrochemical characteristics it can be classified as suitable for agricultural, potable use, or both respectively.

1.1. Groundwater Quality Evaluation

Reliable evaluation of the groundwater quality status is essential to define its proper and sustainable utilization. A valuable and easy tool applied worldwide and scientifically accepted for groundwater quality assessment is represented by Groundwater Quality Indices (GQIs). These indices analyze several hydrochemical factors (chemical, physical, and biological parameters) to compute average groundwater quality. They have been widely utilized due to its pliability in many studies to assess the qualitative status of groundwater or superficial water for different utilization purposes [11,12,13,14]. Along with GQIs, several other groundwater management tools have been proposed to handle and define groundwater pollution such as: (i) Bayesian belief network, to optimize the management of groundwater pollution [15,16], (ii) statistical analysis, to identify the main pollution patterns [7,17,18], and (iii) stable isotopes investigation to correctly define pollution sources [19,20,21,22]. All these tools can produce satisfactory and reliable pollution evaluation allowing local authority to consequently adapt the best water management actions [23,24]. Unfortunately, these methodologies have several drawbacks that could limit their utilization. Primarily, these methods often required a huge amount of field investigations and analysis that can sometimes be expensive and time consuming. Moreover, the identification of pollution sources only represented the tip of the iceberg, the main issue would be: (i) To contain and reduce pollution influx in groundwater and to (ii) steadily manage it especially within already polluted sites. As groundwater pollution remediation is highly expensive and difficult to achieve it becomes essential to avoid it, mainly in areas that are naturally susceptible and intrinsically vulnerable. In this context, vulnerability assessment methodologies have been identified as an easy and effective solution to evaluate specific groundwater vulnerability to general and specific pollution [25,26,27,28].

1.2. Aquifer Vulnerability Assessment

Among all the methodologies created for the aquifer vulnerability assessment purpose, the rating index methods have proven to be the most utilized worldwide due to their ease of use and wide applicability. While the GQIs requires water quality parameters for its computation, the groundwater vulnerability assessment is mainly based on the evaluation of geological and hydrogeological aquifer’s characteristics. Many methodologies have been proposed and further modified over years for vulnerability assessment like the DRASTIC method [29] or SINTACS [30] method, among others. These methods require an evaluation of those parameters directly connected with the intrinsic depuration property of the aquifer system, such as (i) the depth to the water table, (ii) net recharge, (iii) aquifer media, (iv) morphology, (v) vadose zone depuration capacity, and vi) aquifer thickness. Nevertheless, in regard to the methodologies previous described, their application can be highly limited by the lack of field data. Furthermore, they may require deep modification to correctly describe the specific vulnerability to some pollution species. A complete review and discussion of all the above mentioned methodology is available in Machiwa et al. [31]. The methodology entitled “Protection from Natural and Anthropogenic sources” (PNA) proposed by Busico et al. [32] has been created to overcome this issue as it can be applied using increasingly free and available global and regional datasets. Thus, thanks to the main advantages to be a friendly user methodology and to be specific for NO₃ pollution (as it incorporates some parameters like nitrogen soil content (%) and fertilizer utilization for certain kinds of crops) the PNA methodology was utilized to evaluate the intrinsic NO₃ susceptibility of the study area.

1.3. Study Aims

To summarize, the purpose of this study is to obtain a screening status of the hydrological and hydro-chemical regime of the aquifers finalized so as to (i) optimize the exploitation, (ii) to define its better utilization, and (iii) to avoid future pollution and enhance the management of the local groundwater in the considered area. According to the research concept, GQI along with vulnerability assessment and specific hydrogeochemical investigation will be utilized and described. The Flamouria basin has been chosen due to the importance of its groundwater resources used in human activities and due to the complete absence of previous study according to the current literature.

2. Study Area

The study area chosen for this investigation is located in the southern part of the Pella prefecture (Northern Greece) covering an area of 78.08 km² and with a perimeter of 40.97 km. It is mainly a rural irrigated area where all water needs are completely covered by the exploitation of ground and surface water resources. The local climate is generally characterized as continental [33] with an annual precipitation (PCP) amount greater than 800 mm, which is well distributed throughout the year. Widespread farming activities represent the main economical income for the local population and the water need of the area is completely met by the local resources, or through the network of water wells belonging to the Local Organizations of Land Reclamation (LOLR) of Edessa. From a lithological point of view, the area is mainly covered by the Almopian flysch, in a smaller extent by the Kerasia-Kedrona limestones in the eastern part of the area and by the volcano-sedimentary of the ophiolitic range in the southern part. The Almopian subzone dominates the area, with the Pelagonical formations taking the western part. The units of the Almopian subzone that are located within the study area are the middle (Lykoi, Mesimeriou) and western ones (Kerasia, Kedrona). The categorization of the aquifers was made based on the geological formation of the study area and was finally divided into porous, karstic, and fractured types (Figure 1). The groundwater flow is mainly oriented from West to East, feeding the Mega Rema river, which represents the main discharge passage of the basin. The study area is known for the fruit production and is characterized by high agricultural activity (mainly trees of cherry, apricot, and peach). Moreover, livestock units characterize the research area and could be responsible for groundwater quality deterioration and pollution. It should also be mentioned that small industries are located in the research area.

3. Materials and Methods

3.1. Data Collection and Analysis

3.1.1. Groundwater Sampling, Analysis and Electrical Tomography

To fulfill the study’s aims, a census was completed of the available water wells (municipal and agricultural wells) in the area, along with water level measurement during the wet periods of 2017 and 2018 (April–May) and the dry period of 2017 (September–October) to investigate the piezometric conditions and fluctuation in the Flamouria basin (Table S1). The water level measurements were possible in 18 of a total of 37 recorded water wells in the area. Most of the wells are located within the lowlands of Flamouria basin on torrential deposits, while 10 wells have been installed upon flysch formations and another two within ophiolitic formations. The groundwater quality monitoring of the study aquifers include the chemical analysis of 10 groundwater samples. These 10 samples were collected from 7 boreholes, 2 springs (in flysch and gabbro-pyroxenite formations), and 1 surface water sample in correspondence of the watershed’s outlet. The physicochemical parameters of temperature (T), pH, and electrical conductivity (EC) were measured in situ using a multi-parametric probe, HANNA (HI991300). All water samples were filtered using a 0.45 μm Millipore filter and stored in two 50 mL PE bottles (one acidified with ultrapure HNO₃⁻) for laboratory analysis (IC). The chemical analyses of ions concentrations (Na⁺, K⁺, Ca²⁺, Mg²⁺, Cl, SO₄²⁻, HCO₃⁻, and NO₃⁻) were determined via ion chromatography (IC) performed following standard chemical analysis methods [34]. The IC’s calibration curves were obtained using different calibration standards with known concentration and all sample analyses showed a charge balance error within 5%. Moreover, 3 electrical tomographies (ERTs) of 1000 m in length were performed in the Flamouria plain. Each ERT consisting of 21 steel electrodes were positioned with a distance of 50 m away. The dipole-dipole and multigradient arrays were used 3 times and the inductive polarization (IP) measurement was used once. Data were performed using DC2DPRO software extracting 2D images. The mean square error ranged between 4% and 7%, thus the measurements are considered as reliable.

3.1.2. Climatic Characteristics

The main climatic parameters such as T and PCP, were obtained and analyzed using the historical available data provided by the Edessa Weather Station, located 3 kilometers away from the research area. The calculated mean annual PCP considering a 15-year period (2000–2015) was calculated to be 869.9 mm, while the average annual T was 14.1 °C. The average real evapotranspiration (ET) was calculated according to the Thornthwaite method [35,36] following the formula:

∆E = 1.6 × (10 × T/I)^α

(1)

where (∆E) is the average monthly value of ET in mm, (T) represents the average monthly temperature in °C, (I) is the annual thermal index, and (α) indicates a calculated coefficient. All the formulas are available in the Supplementary Materials. The annual Et was calculated to be around 59% of the total PCP with a value of 514.9 mm. Accordingly to the results, the effective recharge was calculated multiplying the value obtained by the difference among P and ET for an infiltration coefficient typical of each geological formation (Table S2) [26,31].

3.2. Protection from Natural and Anthropogenic Sources (PNA)

The PNA methodology considers a variety of available data to produce a NO₃ vulnerability rating for a given region [32]. The parameters are divided into three categories: (i) P are the natural factors that influence NO₃ pollution (depth to water, protection of the unsaturated zone, and slope), (ii) N are those factors that mainly affect NO₃ infiltration into the groundwater (recharge + irrigation, nitrogen soil content), and finally (iii) A express the anthropogenic factors that can alter the rate of NO₃ leaching into the aquifers (land use, irrigation, and well density). The final ratings of the methodology are calculated using the equation:

PNA = (D + U + S) + (N × L) + (I + R + W)

(2)

where D is the depth to the water table, U is the unsaturated zone protection capacity, S is the slope values, N represents the nitrogen soil content (%), L is a score describing the amount of fertilizer suggested for kind of crops, I is the effective recharge, R is the water irrigation amount, andW indicates the well density of the study area. A table showing paramaters classes is available in ESM (Table S3).

3.3. Groundwater Quality Index (GQI)

The analysis of the GQI has been conducted following the methodology proposed by Babiker et al. [37]. In the first step, the contamination index (C) is computed for all water samples, which relate the concentration (X) of sample wells, with its World Health Organization (WHO), Food and Agricultural Organization (FAO), or regional standard (Y), as shown in Equation (3). All the major elements have been utilized to compute the index:

C = (X − Y) / (X + Y)

(3)

The produced C value can range between –1 and 1 and can be transformed in an appropriate rank value (R) via the polynomial Equation (4) to generate ranking score:

R = (0.5 × C ^ 2) + (4.5 × C) + 5

(4)

The R value scales between 1 and 10, with 1 being the lowest impact on the results and 10 being the highest. Following the R calculation for all computed parameters, the final index was calculated using Equation (5):

GQI = 100 − ((R₁W₁ + R₂W₂ +⋯+ R_nW_n) / N)

(5)

where W_n stands for the respective weight of each factor and correlates with the average rating value (Rn) of each rank, with a “mean R + 2” rating for factors that have potential health risks (i.e., NO₃ or F).

4. Results

4.1. Hydrochemical Assessment

The results of the chemical analyses have been investigated using the Durov [38] diagram (Figure 2). The localization of the water sample revealed the main calcium-bicarbonate (Ca-HCO₃⁻) hydrochemical facies for both sample’s belonging to the fractured and porous aquifers. The same chemical signature of the two aquifers probably indicate the strong interconection of the two systems. A general overview of groundwater sample is available in

Table 1. Physico-chemical characteristics of water samples. W indicate wells sample while S indicate springs.

		Water Samples
PARAMETERS	UNITS	ΜW18	MS2	MW40	MS3	MSF1	MW30	MW1	MW5	MW11	MW24
pH		8.00	8.15	7.90	7.70	8.45	7.90	7.80	7.70	8.20	7.60
EC	μS/cm	566	363	585	432	526	709	621	716	458	755
HCO3	mg/L	329.40	225.0	372.10	268.0	286.70	433.10	387.35	329.40	280.60	451.40
Cl	mg/L	4.50	ND	3.10	ND	4.50	6.00	7.00	8.00	ND	4.00
SO4	mg/L	16.00	13.00	11.00	13.00	31.00	36.00	11.00	56.00	8.00	40.00
NO3	mg/L	32.00	ND	9.00	ND	10.30	ND	23.00	53.00	4.00	ND
Na	mg/L	5.80	2.00	8.10	2.90	5.20	7.00	9.00	4.40	2.20	4.10
K	mg/L	1.70	0.30	1.70	1.00	1.90	0.16	1.60	1.10	0.70	3.00
Mg	mg/L	70.00	45.00	70.00	78.40	67.00	79.00	70.00	105.00	55.00	123.50
Ca	mg/L	30.00	19.40	30.80	6.20	25.00	44.00	37.00	27.00	24.40	21.90

. The groundwater samples are characterized as basic, with a mean pH value of 7.9. The electrical conductivity varies between 363 and 755 μS/cm, with the higher values being recorded inside the porous aquifers, while the lowest belong to the fractured ones. The highest values of sulfates were observed at the lowlands (8–56 mg/L with a mean value of 23.5 mg/L) and are probably related to the intensive use of fertilizers [39], while in the northern part they probably originate from the local ophiolite formations. The highest NO₃ concentrations were observed inside the shallow aquifers in the plain section, as mentioned before, where fertilizers are extensively used in agricultural activities, while the mountain part is characterized by lower concentrations. The concentrations of bicarbonates (HCO₃⁻), calcium (Ca²⁺), and magnesium (Mg²⁺) are ascribable to the presence of carbonate formations, dolomites, marbles, and ophiolites. The Mg/Ca, (Ca + Mg)/(Na + K), Cl/SO₄, and Cl/(CO₃ + HCO₃) ion ratios have been selected to determine the origin of the ions in the water solution (Table S3). The meq/L ratio of Mg/Ca varies from 0.13 to 0.92 with a mean value of 0.61 meq/L and a standard deviation of 0.24. The results related these elements to the presence of karst aquifers in the area, which are unloading through the fractured aquifers. Moreover, the increased values of the (Ca + Mg)/(Na + K) ratio with a mean value of 27.14 probably indicate that the aquifer is being constantly recharged (recharge zone). The meq/L ratio of Cl/SO₄ vary between 0.1 to 0.86 with a mean value of 0.28 and corresponds to a chloride–chlosulphide water type. Finally, the meq/L ratio of Cl/(CO₃ + HCO₃) range between 0.01 to 0.04 with a mean value of 0.02, indicating a good groundwater quality. Additionally, the concentrations of δD and δ18O were detected in the groundwater samples plots along the meteoric water line Hellas (HellasMWL). The relation of the Greek meteoric water line (Hellas MWL) is given by: δD = 7.2 × δ18O + 8.2 ‰. According to the results of the stable isotopes, the water is characterized by a meteoric origin that has not been subjected to any secondary evaporation, while the maximum recharge altitude goes up to 1194 m (Figure S1 and Supplementary Information 2). According to the water level measurements, artesian phenomena have been observed in the center of the lowlands, especially during the wet periods and a high fluctuation of the water levels was also noted between the dry and wet seasons, probably due to the elevation differences between the various aquifers. The water level measurements have been used to determine the piezometric regime of the area. The lowest hydraulic head values were registered in the lowlands of the basin, while the highest values were recorded at the eastern mountainous area. Despite the rapid increase in water usage during the summer months for the coverage of agricultural, stock raising, and general water needs, the groundwater reserves are being recharged to their original levels during the wet period and especially in the highland sector. For instance, an increase of 6 m in the water level during measurements for the 2018 wet period, compared to the 2017 measurements of the same period was registered

4.2. Groundwater Quality Evaluation

The overall water quality has been calculated using the GQI and graphic reoresentation. According to the Wilcox diagram (Figure 3a), the quality of the water samples is described from good to excellent. Moreover, the Richards diagram (Figure 3b) indicates that the groundwater is also suitable for agriculture purposes. To enforce the excellent groundwater status highlighted by the previous diagram, a GQI evaluation has been applied to assess the groundwater suitability for potable and irrigation use. The GQI has been calculated following the procedure proposed by Babiker et al. [37] utilizing the available chemical parameters. In

Table 2. Statistics of the parameters rank utilized to compute Groundwater Quality Index (GQI).

Parameters	WHO Limit (mg/L)	FAO Standard (mg/L)	Potable Use Average R	Irrigation Use Average R
Ca2+	300	-	2.49	-
Mg2+	300	-	1.58	1.58
Na+	200	69	1.17	1.48
Cl−	200	107	1.14	1.25
NO3−	50	5	2.30	5.21
SO42−	250	-	1.60	-
TDS	600	450	4.00	4.60
HCO3−	-	91.50	-	7.68

, the utilized threshold limits for potable (WHO) and irrigation use (FAO), along with the average calculated weights (R) of each parameters are shown.

The GQI for potable use further confirmed the previous results, characterizing all water samples as high-quality water for potable use with a GQI ranging from 91–96 (on a maximum of 100). In this case, the parameters with the higher negative impacts are the TDS (rank 4), Ca²⁺ (rank 2.5), and NO₃ (rank 2.30) reaching in some water samples concentratioins closer to WHO’s limit but never exceeding it. Concerning the irrigation use instead, the water samples always showed a medium-high quality but with a GQI ranging from 70 to 84, which is much lower when compared to the potable classification. In this latter case, the parametesr of HCO₃, TDS, and NO₃ showed the highest mean rank values. The more restricted values of FAO guidelines compared to WHO standards contributed in lowering the GQI for irrigation use.

4.3. Geoelectrical Data Evaluation

The geoelectrical measurements were performed to verify the geometrical characterization and to determine the main groudwater flow [40]. Three electric resisitivity tomographies (ERTs) were performed upon sedimentary formations (pebbles, sand, and clay), with a length of 1000 m, giving back 250 m depth of geoelectrical structure (red line in Figure 1). The evaluation of the various formations has been done accordingly to their electrical resistivity, while a great importance was given to the conductive formations, where the aquifers were detected. The vertical two-dimensional images of groundwater zones along with the possible faults presence and collection of the hydrochemical data, resulted in a solid conclusion regarding the hydrogeological and hydrochemical regime of the area. Specifically, a conductive body (25–37 Ohm-m) was observed in the section from 310 to 450 m of ERT1 with a width of 100 m and a depth between 65–150 m. The high resistivity was recorded (46–66 Ohm-m) due to the presence of a fault zone in the porous formation consisting of bedded tuffites. A second conductive formation is located at a depth of 150 m (20–37 Ohm-m) between 350–600 m, depicting the under-pressure karstic aquifer located within the limestones of the region. The presence of this karstic aquifer is confirmed in ERT2, where a conductive body (20–36 Ohm-m) was located at a depth of 160 m and a distance of 450 to 650 m from the start of the tomography. Moreover, a subsurface conductive formation was observed between 115 and 550 m from the start of the tomography, which is an unrestricted porous aquifer that resides between various bodies of pebbles. Lastly, an oversupply of conductive formations were observed in ERT3 tomography allocated as follows: (i) A surface conductive body between 200 and 450 m, characterized by an unrestricted porous aquifer located within pebbles, with a depth of up to 30 m, and (ii) a secondary conductive formation extends throughout the tomography between the depth of 40 to 130 m adescribed by a porous aquifer, resides within sand and gravel. Finally, a third conductive formation located at a depth greater than 170 m confirmed the presence of the pressurized karst aquifer.

4.4. “Nitrate Vulnerable Zone” Delineation

To evaluate the intrinsic susceptibility of the study area to NO₃ pollution, accordingly with its hydrological and socio-economical characteristics, an easy-to-use methodology called PNA to delineate “Nitrate Vulnerable Zone” was applied to the study area. All the information necessary for the evaluation has been collected through field data analysis, literature review, and free access dataset and then stored and elaborated in the GIS environment. Seven parameters have been evaluated for the final index and digitalized in raster format retaining a cell resolution of 10 × 10 m (Figure 4).

4.4.1. Depth to Water Table (D)

For the realization of the water depth evaluation map, an average value among the two sampling campaigns was utilized (∆GL column in Table S1). Afterwards the kriging function in the ArcGIS environment was utilized to forecast the depth to the water table for the entire alluvial basin. Using the previously described results, the study area is divided in three main classes of vulnerability (Figure 5a): (i) The low vulnerability assigned to the higher depth of the water level in the south area, (ii) the higher vulnerability, instead mainly spreads in the north-east sector on the plain characterized by a shallower water level, and finally (iii) the average vulnerability assigned to the center of the area. Summarizing, the deeper water table is located on the western part, with a depth to water ranging between 20 and 50 m (less vulnerable), while on the eastern side the groundwater level lies between 0 and 5 m (more vulnerable).

4.4.2. Slope Evaluation (S)

The morphology of the study area is evaluated starting from a Digital Elevation Model “DEM” with a cell resolution of 20 × 20 m (Figure 5b). The correspondence slope values in percentage (%) is calculated using the “slope” toolbox in ArcGIS 10.5. Those areas relative flat are naturally more prone to pollutant’s infiltration and consequently they suffer of a higher vulnerability compared to those areas with a bigger steep. As expected, the higher ratings are found within the plains between the Flamouria and Platani villages, while lower ratings were observed on the slopes of nearby mountainous terrains.

4.4.3. Vadose Zone Protection (U)

The unsaturated zone represents the main aquifer’s protection against pollutant leaching. The information about the unsaturated zone of the study area is obtained analyzing the results of electrical tomography study and from the available geological map. Accordingly, the lowland area is mainly characterized by reworked volcanic formations (tuffites) characterized by medium-high vulnerability (Figure 5c). Compact formations like metamorphic and old volcanic rocks usually offer better protection to the underlying vadose zone due to their low permeability indicating a lower vulnerability and they cover the southern part of the basin.

4.4.4. Wells Density

The wells density map is created by counting the amounts of wells present on one square kilometer of the study area. Fewer wells represent lower aquifer exploitation and so lower vulnerability, while the higher wells density indicating higher exploitation and anthropogenic impact, increasing the vulnerability. According to the index, in the study area, the lower rating is coincident with the center of the plain, where most agricultural activities take place (Figure 5d)

4.4.5. Recharge and Irrigation (R + I)

The recharge and irrigation map is created taking into consideration the recharge rate of the local aquifers and agricultural needs in irrigational water. Accordingly, for elaboration, the value of recharge upon the porous aquifer considering both precipitation and irrigation water is estimated to be around the 250 mm/y, retaining the higher vulnerability in the classification (Figure 5e).

4.4.6. Nitrogen Soil Content (N)

In these parameters the nitrogen soil content (%) represents the amount of nitrogen naturally occurring in the soil matrix. It is evaluated accordingly with the soil types identified using the Digital Soil World Map (DSWP, FAO). Despite the thinner soil thickness, it is characterized by a good amount of Nitrogen, (>1.0%) and accordingly to this evaluation, a rate of 0.96 has been utilized for the entire study area (Figure 5f).

4.4.7. Land Use (L)

Based on the Corine Land Cover (CLC) classification of 2018, most of the Flamouria basin is covered by forests (54.5%), followed by fruit trees plantation and fleshy fruits (18.24%), and generical agricultural lands (12.3%). The agricultural sector is the main source of income for the inhabitants of the area, mainly from the cherries and peach plantations (almost 90% of total agricultural usage), while the remaining 10% is represented by apples, kiwis, olives, chestnuts, and various micro-cultivations. The land cover evaluation for the analysis of “Nitrate Vulnerable Zone” was conducted analyzing the average nitrogen supply necessary for the optimum growth of a specific crops. The lowest ratings are allocated in lands covered by forests, pastures, or plain bare lands as they did not require any anthropogenic nitrogen input (Figure 5g). Plantations of trees that do not need that much care to grow, like olive groves, also show a low rating. On the other hand, orchards of fruit trees and fleshy fruits have higher nitrogen needs and are allocated at the mid-range of the ratings. The worst ratings of the land use evaluation accordingly with the table proposed by Busico et al. [32] are located at the urban areas

4.4.8. PNA Map Evaluation

Based on the discussed parameters, a final map was created combining all the created layers together through linear combination following Equation (2) to produce the final PNA map for the study area (Figure 6). The site has been divided into five classes of vulnerability going from very low to very high and the class ranges has been created using the geometrical interval. As expected, the areas with lower susceptibility to NO₃ pollution are those located in the hillsides around the plain, between the villages of Flamouria and Platani (green color in Figure 6), while the susceptibility greatly increases towards the urban areas and plain (red color in Figure 6), due to the presence of the agricultural and anthropogenic activities that take place in these areas. The higher vulnerability to NO₃⁻ located to the center of the plain is not only connected to the anthropogenic activities but also to the natural site characteristics as this zone, which showed the lowest depth to water table, the highest infiltration, and flat topography.

5. Discussions

All the main hydrological and hydro-chemical characteristics of the studied area have been established using all the available data and tools. According to the main water chemical characteristic and comparing the results obtained with ERT, the Edessa basin is characterized by a deeply interconnected multi-layer aquifer system. Their similar chemical composition, both retaining a HCO₃-Ca-facies could indicate that the porous aquifer obtains a constant recharge from the fractured one. The overall groundwater quality of the study area has been classified as “excellent” for potable use, according to the various applied methodologies. A different situation has been observed for irrigation use with a medium-high quality. The high amount of HCO₃⁻ and the total TDS are responsible of water quality deterioration for irrigation purpose. In particular, according to FAO guidelines: (i) The total salinity (TDS) can negatively influence the water availability for the crops, (ii) sodium and chloride can generate negative growth responses in some sensitive crops even at low concentrations, and (iii) NO₃ and HCO₃⁻ can cause a series of various negative effects such as excessive nutrients and precipitates on leaves, respectively. Thus, the continuous utilization of such water for irrigation, especially where the HCO₃⁻ concentration exceeds 5 meq/L, could bring to the necessity to neutralize the HCO₃⁻ for long-term irrigation purposes [41]. It is also clear how the main productive sector of the basin is represented by the agricultural activities and especially by cherry and orchards plantation. Many studies have showed the influence of land-use patterns, the type of aquifer, and soil-drainage capacity on nitrate pollution [4,25,42], as reported to the site characteristic and data availability the selection of the most suitable methodology for the delineation of groundwater vulnerability zones is critically important to define optimal management strategies [43]. For this purpose, the PNA methodology fits well with the study requirements as it includes in the vulnerability evaluation all the parameters connected to the land use patterns, aquifer, climatic, and soil capacity. Despite the high-quality level of the analyzed groundwater, computed with GQI, the PNA classification shows a completely different situation, characterizing the entire plain as highly susceptible to the nitrate pollution. The high susceptibility to the NO₃ pollution can be mainly ascribable to flat topography, high recharge, and to the medium high permeability of the vadose zone. It worth mentioning that areas characterized with the higher susceptibility to NO₃ pollution are not necessarily the same areas with higher nitrate concentrations. The results are also in agreement with previous elaboration in other regions of the world [32,44,45] where the more susceptible are to NO₃ pollution are those with a lower pollutant concentration, but that could become highly polluted in the future due to natural susceptibility, highlighting the importance of land use practice and their pressure or pollution problem. In our case a future increase of NO₃ concentration can be expected especially within porous aquifer, despite the relatively small nitrogen amounts necessary for orchards to grow. Summarizing, this can lead to a worsening scenario in the future if a sustainable land management planning is not applied along with a rigorous groundwater monitoring from a quantitative and qualitative point of view. Undoubtedly, the present work further confirmed the high flexibility of PNA methodology in correctly describing the aquifer’s vulnerability to NO₃ pollution, especially in highly cultivated agricultural lands. It is worth mentioning that it is necessary to gather more detailed analysis and information so as to obtain more detailed hydrogeological settings. A complete reconstruction of the underground system as well as a quantification of stream recharge could further improve the obtained hydrogeological framework. Anyway, the proposed elaboration represents an easy and widely-applicable assessment where limited data are also available.

6. Conclusions

The following conclusions can be exported from the analysis of the geological, hydrogeological, hydrochemical, and geophysical data of the Flamouria basin, Northern Greece. All these findings further remarked the importance of achieving an accurate hydrogeological survey to highlight all the main features and issues of an aquifer system. Summarizing five main findings have been highlighted by this study:

• The porous aquifer in granular rocks is being supplied laterally by the karst and fractured aquifers;
• The rational usage of fertilizers is essential, as the highest values of nitrates and sulfates were recorded inside the Flamouria plain, where intensive agricultural activities take place along with the presence of the boreholes for water supply;
• According to the PNA method results, the most vulnerable locations of the study area coincide with the main urban and agricultural zones of the Flamouria and Platani villages and the plain region between these two settlements;
• Based on the hydrochemical results and GQI application, the groundwater is deemed as of high quality for potable use and medium-high quality for irrigative purposes;
• The application of electrical tomographies and the processing of their data helped to locate underground aquifers as well as possible faults of the study area, which can explain the hydrological regime of the basin.

Lastly, future monitoring of the system is advised, both qualitatively and quantitatively, as well as the estimation of the vulnerability of the other aquifers (karst and fractured one) to better achieve a sustainable use of groundwater resources. Moreover, the assessment remarked how all the applied tools can be easily and widely applied around the world to evaluate actual and future issues related to groundwater sustainability.