**Geo-Environmental Approaches for the Analysis and Assessment of Groundwater Resources at Catchment-Scale**

Editors

**Evangelos Tziritis Andreas Panagopoulos**

MDPI ' Basel ' Beijing ' Wuhan ' Barcelona ' Belgrade ' Manchester ' Tokyo ' Cluj ' Tianjin

*Editors* Evangelos Tziritis Soil and Water Resources Institute ELGO-DIMITRA Thessaloniki Greece

Andreas Panagopoulos Soil and Water Resources Institute ELGO-DIMITRA Thessaloniki Greece

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Water* (ISSN 2073-4441) (available at: www.mdpi.com/journal/water/special issues/Geo-environmental groundwater catchment).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-4372-7 (Hbk) ISBN 978-3-0365-4371-0 (PDF)**

Cover image courtesy of Evangelos Tziritis

© 2022 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**

## **About the Editors**

#### **Evangelos Tziritis**

Dr Evangelos Tziritis is a senior research scientist at the Soil and Water Resources Institute of the Hellenic Agricultural Organization (ELGO-DIMITRA). His main research domain is focused on environmental hydrogeochemistry and relevant water resources aspects. He has significant experience in various methodological approaches for geo-environmental analysis at the catchment-scale level (e.g., hydrogeochemical modelling, environmental indicators, hydrological isotopes, groundwater vulnerability). He also has expert knowledge on environmental monitoring of water resources and groundwater salinization. His overall record of achievements includes several years of experience in geo-environmental projects of basic and applied research in liaison with private firms, stakeholders, and academia. He has participated in more than 20 projects and has published more than 60 scientific papers (h-index 15) in peer-reviewed journals and international conferences. He is a regular reviewer of more than 30 scientific journals of the Science Citation Index (SCI). Currently, he is the coordinator of the MEDSAL Project (PRIMA Programme), focusing on the multi-sourced salinization of coastal aquifers.

#### **Andreas Panagopoulos**

Dr. Andreas Panagopoulos is a Hydrogeologist, Research Director at the Soil & Water Resources Institute (SWRI) of the Hellenic Agricultural Organisaton. His main scientific interests concern water resources management, monitoring, modelling and policy, hydrogeology and climate change impact assessment in water bodies. With more than 25 years of working experience, he has participated in more than 40 national and international research projects, in most of which he was scientifically responsible for SWRI. He has worked as a consultant for European and Greek partnerships in the implementation of development studies in Greece and Turkey (EUROPAID) concerning the protection and management of water resources in the framework of the implementation of European Directives. He is a member of the National Water Council and a nominated national representative in the CIS WG groundwater of the EU. He is also the co-founder of Pinios Hydrologic Observatory, which is included in International Long Term Ecological Research (ILTER) sites and the European Network of Hydrological Observatories (ENOHA) network. He has been the Director of Land.

### *Editorial* **Geo-Environmental Approaches for the Analysis and Assessment of Groundwater Resources at the Catchment Scale**

**Evangelos Tziritis \* and Andreas Panagopoulos \***

Soil and Water Resources Institute, ELGO-DIMITRA, 57400 Thessaloniki, Greece **\*** Correspondence<EMAIL_ADDRESS>(E.T<EMAIL_ADDRESS>(A.P.)

Groundwater resources constitute nearly one-third of the globe's freshwater resources. They are widely used in national economies for various key purposes, including domestic water supply, industrial uses, irrigation, medicinal uses (balneotherapy), and energy applications (low- and high-enthalpy geothermal energy). Moreover, they often form the only supply to dependent ecosystems and give rise to spectacular underground karstic formations (lagoons, river canyons, etc.). In addition, they provide the necessary mechanical strength to lose formations that host phreatic aquifers, thus preventing matrix compaction and subsidence at the ground surface that causes devastating deformations and extensive damage to infrastructures. Hence, groundwater's quantitative sufficiency and high quality are paramount to secure environmental sustainability and socio-economic prosperity. Nevertheless, groundwater systems are subject to various factors that are adverse to or deteriorate their inherent characteristics. Such factors could be geogenic or anthropogenic; geogenic factors are mainly related to natural water–rock interaction processes, whereas anthropogenic factors are primarily associated with agricultural practices, industrial activities, urbanisation, landfills, domestic effluents, and aquifer overexploitation.

Another critical aspect of groundwater management is that the proper management of groundwater is difficult after contamination. On the other hand, environmental monitoring projects that provide critical ground-truth values for groundwater quality are not always feasible due to a lack of personnel, funding, and time. In this respect, scientists and decision makers seek alternative strategic tools to spatially identify threats and subsequently design and implement groundwater protection measures towards practicing sustainable groundwater management. These tools are often diverse but interlinked methodologies, such as hydrogeological and hydrogeochemical modelling, environmental isotopes, environmental indicators, geostatistics, and artificial intelligence.

The complexity of different hydrological and hydrogeological setups, hydrodynamic patterns, site specifications, and the wide variability of internal and external factors and/or processes at the level of catchment scale necessitates combined approaches that integrate robust methods, thus leading to more accurate and reliable outcomes towards sustainable groundwater management. Sound knowledge of a studied groundwater system may reduce uncertainties in the prediction of its future evolution, thus enabling better management and protection whilst limiting the need to hypothesise. In line with the above goal, this Special Issue aims to provide successful applications or new insights on the standalone or joint considerations of groundwater resources assessment and characterisation methods and explore new state-of-the-art methodological concepts in light of a rapidly changing environment.

This Special Issue of the journal *Water* comprises 12 papers with contributions of more than 50 authors originating from 10 countries, all of which deal with various geoenvironmental methods and tools. The papers include six feature papers and four editor's choice papers.

Two of these papers apply hydrogeochemical tools and methods to robustly evaluate groundwater resources quality and hydrodynamic regime. Specifically, Vrouhakis et al. [1]

**Citation:** Tziritis, E.; Panagopoulos, A. Geo-Environmental Approaches for the Analysis and Assessment of Groundwater Resources at the Catchment Scale. *Water* **2022**, *14*, 1085. https://doi.org/10.3390/w14071085

Received: 2 March 2022 Accepted: 16 March 2022 Published: 29 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

used combined hydrogeochemical and hydrodynamic characterisation to assess key aspects related to groundwater resources management in a highly productive agricultural basin of central Greece. A complementary array of tools and methods, including graphical processing, multivariate statistics, and environmental isotopes, was applied to a comprehensive dataset of physicochemical analyses and water level measurements. The outcomes proved valuable in the progression towards sustainable management of groundwater resources. The results provide spatial and temporal insights into significant parameters, sources, and processes that, as a methodological approach, could be adopted in similar cases of other catchments. Vasileiou et al. [2] investigated the hydrogeochemical processes and natural background levels (NBLs) of chromium in the ultramafic environment of Vermio Mountain, Western Macedonia, Greece. The holistic methodology proposed in this paper may be implemented in similar cases at the catchment scale to assess geogenic and anthropogenic Cr sources that degrade groundwater quality.

Falling within the broad category of groundwater quality, two papers deal with the significant problem of groundwater salinisation, which is crucial to address in arid and semiarid coastal areas. More specifically, Bonamico et al. [3] applied an integrated hydrogeologic and geochemical methods approach to describe freshwater–saltwater interactions in a coastal aquifer in the Ostia Antica archaeological park Roma, Italy. Their assessments were based on a water monitoring program that included the installation of multiparametric probes in wells, with continuous measurement of temperature, electrical conductivity, and water table level. Field surveys, water sampling, major elements, and bromide analyses were carried out to understand the detailed stratigraphic settings of the area. The authors used oxygen and carbon isotopic signatures of calcite from well sediments and evaluated major elements and Br to determine the salinisation sources and the processes of gas–waterrock interaction. The second paper focuses on assessing the drought index as an indicator of groundwater salinisation. The research by Alfio et al. [4] was performed on the Salento aquifer (southern Italy), where in recent decades, groundwater depletion and salinisation worsened because of the increased frequency of droughts, as revealed by the documented data derived from the analysis of the Standardised Precipitation Index (SPI) calculations during 1949–2011 based on monthly precipitation levels. Groundwater level series and chloride concentrations collected over the extreme drought period 1989–1990 allowed a qualitative assessment of groundwater behaviour, highlighting the concurrent groundwater drought and salinisation.

Two other papers deal with the conceptualisation and application of index and overlay methods to assess groundwater vulnerability, considered a significant, proactive measure towards successful decision making and rationale management. Kim et al. [5] evaluated the vulnerability to seawater intrusion by classifying the existing GALDIT method into static parameters (groundwater occurrence (G), aquifer hydraulic conductivity (A), and distance from shore (D)) and dynamic parameters (height-to-groundwater level above sea level (L), impact of existing status of seawater intrusion (I), and aquifer thickness (T)). Data indicating averages of measurements over a 10-year period for each month were used, representing the seasonal characteristics of local water cycles. To reflect subtle monthly variations, the range of scores was divided into deciles to capture the temporal dynamics of seawater intrusion. The proposed modified method can determine where to apply countermeasures to vulnerable coastal areas and develop water resources management plans considering vulnerable seasons. Vrouhakis et al. [6] assessed the intrinsic perspective of groundwater vulnerability in central Greece's highly productive agricultural area. A novel index-based method (RIVA) was applied to the Tirnavos basin to assess susceptibility to surface-released contamination. Data from field surveys, previous studies, and the relevant literature were used to calculate factors constituting the RIVA method, which was demonstrated to be a data-intensive and efficient method, thus a sound investment to reach highly accurate results. Overall, RIVA proved to be a robust tool for reliable groundwater vulnerability assessments and could be further exploited for risk assessment and decision-making processes in the context of groundwater resource management.

Fuentes-Arreazola et al. [7] estimated aquifer parameters in the Mexican wine-producing region Guadelupe Valley from fluctuations in levels of groundwater. They proposed an alternative tool with significant advantages in studying the groundwater-level response due to variations in pore pressure caused by internal deformation of the aquifer structure induced by barometric pressure and solid Earth tide. This analysis reveals helpful insights that can help to establish a framework to design and assess management strategies for groundwater resources in similar cases.

Two other papers provided insights into the application of modelling in groundwater resources. Lyra et al. [8] presented an integrated modelling system to evaluate the availability of water resources in coastal agricultural watersheds. Their modelling system was made from an ensemble of surface and groundwater hydrology models, crop growth/nitrate leaching, contaminant transport, and seawater models. Its efficacy to simulate the quantity and quality of water resources was tested at the Almyros basin in Thessaly, Greece. The proposed modelling system could be used as a tool for the simulation of water resources management and climate change scenarios. Further research on the plausibility of modelled nitrate concentrations in the leachate on the scale of federal states in Germany was performed by Wolters et al. [9]. This research aimed to model nitrate concentrations in leachates as a robust tool for water resources management, in line with the requirements enforced by the EU Water Framework Directive. The validity of simulations was checked against values from 1119 preselected monitoring stations from shallow springs and aquifers filtered near the surface with oxidising properties. The case study revealed that the applied model system (RAUMIS–mGROWA–DENUZ) can reliably represent interrelationships and influencing factors that determine simulated nitrate concentrations in the leachate. Moreover, it was demonstrated that observed nitrate concentrations in groundwater may provide a solid source of data for checking the plausibility of modelled nitrate concentrations in leachate in cases in which certain preselection criteria are applied.

Pisinaras et al. [10] studied the effects of agricultural water management in a Mediterranean coastal aquifer under current and projected climate change conditions. Their research focused on the coastal delta plain of River Pinios, central Greece. Such areas are significant for the Mediterranean region because of their high soil fertility and agricultural productivity. Nevertheless, they also constitute fragile systems in terms of water resource management due to the interaction of underlying aquifers with the sea. Soil and Water Assessment Tool (SWAT) and SEAWAT models were combined to simulate the impact of current practices in water resource management on the main groundwater budget components and groundwater salinisation of the shallow aquifer developed in the area. Moreover, the potential impact of climate change was investigated using projected data gathered from the Regional Climate Model for two periods (2021–2050 and 2071–2100) and two sea-level rise scenarios (increasing by 0.5 and 1 m).

The combined impact of the hydrological and socio-economic perspectives on the sustainability of groundwater resources was examined by Oke and Alowo [11]. Their research presented a spatial interpolation of the anticipated impact of the above factors on groundwater systems and predicted the sustainability of the Modder River catchment in South Africa. The results were presented with sustainability maps indicating areas with differing groundwater dynamics in the catchment. The key finding in this paper may assist groundwater managers and regulators to effectively plan groundwater resources utilisation, especially with regard to the prevention of licencing and overpumping practices.

Finally, a relatively new and challenging field of joint tools for hydrological sciences was addressed by Namous et al. [12]. Their research focused on the spatial prediction of groundwater potentiality in a large semi-arid karstic mountainous region by using a combination of machine learning models, such as random forest (RF), logistic regression (LR), decision tree (DT), and artificial neural networks (ANNs). A total of 24 groundwater influencing factors (GIFs) were selected based on a multicollinearity test and the information gain calculation. The results of the groundwater potentiality mapping were validated using statistical measures and the receiver operating characteristic curve (ROC) method. Compared with individual models, the combined models proved the most stable and suitable tools to map groundwater potentiality in mountainous aquifers, based on success and prediction rate.

The sustainability and environmental welfare of human civilization are based on utilising sufficient volumes of water with acceptable quality. As the impacts of climate change intensify, water resources safety becomes of primary concern. Groundwater resources impose extra challenges regarding their management, as these sources are invisible, thus needing proactive and carefully designed measures. World Water Day for the year 2022 is devoted to making the invisible, visible. This Special Issue provides a series of papers that propose state-of-the-art methodologies, technologies, and approaches that exactly contribute to this goal.

**Author Contributions:** E.T. and A.P. have equally contributed to the original draft preparation, review, editing and proof-reading. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors appreciate the efforts of the *Water Journal* editors and publication team at the MDPI and the anonymous reviewers for their invaluable comments.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

### *Article* **Groundwater Vulnerability Analysis of Tirnavos Basin, Central Greece: An Application of RIVA Method**

**Ioannis Vrouhakis 1,2, \* , Evangelos Tziritis 2 , Georgios Stamatis 1 and Andreas Panagopoulos 2**

**Abstract:** A novel index-based method (RIVA) for assessing intrinsic groundwater vulnerability was applied to Tirnavos basin (central Greece) to assess the susceptibility to surface-released contamination. Data from field surveys, previous studies, and literature were used to calculate the factors that compile the RIVA method. The aggregated results delineated the spatial distribution of groundwater vulnerability from very low to very high. The modelled results were successfully validated with ground-truth values of nitrates obtained from 43 boreholes. Overall, the modelled and the monitored values match more than 80%, indicating the successful application of the RIVA method. Few deviations were observed in areas dominantly affected by lateral crossflows and contamination from adjacent areas. RIVA proved an efficient method in terms of accuracy, data intensity, and investment to reach highly accurate results. Overall, RIVA proved to be a robust tool for reliable groundwater vulnerability assessments and could be further exploited for risk assessment and decision-making processes in the context of groundwater resource management.

**Keywords:** aquifer; groundwater; intrinsic vulnerability; RIVA method; index-overlay method; Tirnavos basin

### **1. Introduction**

Groundwater accounts for nearly 99% of the total volume of freshwater presently circulating on our planet [1]. This overwhelming percentage, the lower susceptibility to pollution, and the large storage capacity of groundwater compared to surface water highlight its paramount importance at a global socio-economic level. Its significance, however, requires the most outstanding possible effort to protect it. A significant task towards this goal is the assessment of groundwater intrinsic vulnerability, a practice that suggests areas where priority and special attention should be given to the protection and overall management of groundwater resources. However, assessing the vulnerability of groundwater to adverse effects of human impacts is one of the most critical problems in applied hydrogeology [2]. The anthropogenic agricultural activities are often responsible for overdraft, groundwater quality deterioration, and increasing vulnerability. Due to level decline and quality degradation, sustainable development plans are needed to protect these resources [3].

In general, in most parts of the world, groundwater vulnerability assessment is based on (i) process-based methods [4–6], (ii) statistical methods [7–9], and (iii) overlay and index methods [10–13]. The limitations of process-based methods are adequate data and quality to capture the physical, chemical, and biological reactions from the surface to the uppermost aquifer. The statistical methods focus on the uncertainty by minimizing the error and using parameters' coefficients instead of weights. A possible drawback of these methods is the required monitoring data, which is essential. These methods are only applicable to those

**Citation:** Vrouhakis, I.; Tziritis, E.; Stamatis, G.; Panagopoulos, A. Groundwater Vulnerability Analysis of Tirnavos Basin, Central Greece: An Application of RIVA Method. *Water* **2022**, *14*, 534. https://doi.org/ 10.3390/w14040534

Academic Editor: Pankaj Kumar

Received: 19 January 2022 Accepted: 9 February 2022 Published: 11 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

regions where similar factors govern groundwater contamination. The overlay and index methods are the most suitable for groundwater vulnerability assessment, overcoming all the limitations mentioned above [14]. They focus on critical factors potentially controlling contaminant transport, and they are relatively cost-effective and adaptable to on-site specific conditions. Moreover, they demand minimum data to produce outcomes that can directly facilitate the decision-making processes [15]. The basic steps of these methods include raw data analysis, ranking of features on maps, integration of maps, and classification of the integrated map based on an index. These methods can be applied from a regional to global scale; nevertheless, they should be supplemented with field visits and on-site validation to produce reliable results [16].

Due to the advantages above and their ability to be easily used as strategic tools, the scientific community has increased interest. Some examples of this type of groundwater vulnerability assessment are DRASTIC (D: aquifer depth, R: recharge rate, A: aquifer lithology, S: soil type, T: topography, I: impact of vadose zone, C: aquifer hydraulic conductivity) [17], GOD (Groundwater occurrence, Overall lithology of aquifer, and Depth to groundwater level) [18], AVI (Aquifer Vulnerability Index) [19], SEEPAGE (System of Early Evaluation of Pollution Potential of Agricultural Groundwater Environments) [20], EPIK (Epikarst, Protective cover, Infiltration conditions, and Karst network development) [21], GALDIT (G: Groundwater occurrence, A: Aquifer hydraulic conductivity, L: Height of groundwater level, D: Distance from the shore, I: Impact of existing status of seawater intrusion, T: Thickness of the aquifer) [22], RISKE (Rock of aquifer media, Infiltration, Soil media, Karst, and Epikarst) [23], and a global risk approach [24].

Two types of aquifer vulnerability are considered. The first refers to intrinsic vulnerability, which considers the system's inherent properties, such as the properties of the vadose zone or the recharge conditions, etc. The second and more complex one is specific vulnerability. In addition to the intrinsic, it considers the properties of a specific contaminant or group of contaminants [25]. This research focuses on intrinsic vulnerability using a relatively new method (RIVA) [26]. The selected method is based on the successful concept of the European approach [27] and incorporates additional elements that provide more realistic and representative results. The main advantage of RIVA over similar methods is that it can be applied to all types of groundwater bodies independently of the specific conditions, lithologic phases, and aquifer typology of each area. This way allows uniform evaluation and comparison between hydrogeological systems with different characteristics. To this aim, the Tirnavos alluvial basin area is considered an appropriate test site for its application, as it includes an intensively cultivated area with diverse geological and hydrogeological characteristics. The proximity to the adjacent karstic system, combined with the occurrence of significant tectonic structures, further increases the aquifer system's vulnerability potential.

#### **2. Case Study Area**

Thessaly Plain in central Greece is the largest alluvial basin of the country and is divided, through the mid-Thessaly hills, into two sub-basins, the western Thessaly basin and the eastern Thessaly basin, developed in an NW-SE direction as part of the broader tectonic trough. The case study area where the RIVA method was applied is located in the northwest part of the wider eastern Thessaly basin (Figure 1). The hydrological perspective includes parts of the Titarisios River basin to the north and the Pinios River basin to the south. The total area is estimated to be 688 km2 , including the alluvial Tirvanos basin and part of its surrounding geological formations (Figure 2), which are hydrogeologically interdependent. The perimeter of the study area is 105 km, and the mean altitude is 160 m. The smallest and largest morphological gradients recorded are 0 and 62.57%, respectively, with a mean slope of 10.5%.

**Figure 1.** Location map of the case study area where the RIVA method has been applied.

The Tirnavos area is characterised by a typical Mediterranean climate, with annual rainfall from 400 mm to 600 mm, distributed almost entirely during the wet hydrological period, without significant summer precipitation. Larissa station is the closest one to the study area regarding the meteorological data, with continuous and reliable data over several years. Adjacent meteorological stations, used later in the methodology, are in Elassona and Tirnavos.

It is primarily a rural area covered by agricultural land, where intensified agricultural activities, both cultivation and livestock, are a significant source of groundwater contamination by nitrogen compounds. Manure waste and the often excessive and improper use of nitrogen fertilizers, aiming to improve agricultural production, lead to the occurrence of elevated concentrations of nitrates in groundwater [28]. Figure 3 shows the percentage distribution of land use according to CORINE categorization [29], suggesting that about 80% of land use consists of agricultural areas (arable land and permanent crops). The most common irrigation methods are drip irrigation and sprinkler, mainly using groundwater resources.

**Figure 2.** Geological map of the case study area, based on [30,31].

**Figure 3.** Land-uses distribution at the study area.

Three local irrigation organisations (Tirnavou, Ambelona, Agia Sofia) cater for irrigation water to producers through collective networks within their districts. Nowadays, all operating networks are pressurised, while until recently, an extensive gravitational network operated on the eastern part of the Ampelona district utilizing the karstic spring of Mati Tirnavou until 1998, when it was abandoned because of diminishing spring discharge. Despite the collective irrigation works, several privately owned wells exist and operate to cover irrigation demands of land that does not fall into the jurisdiction of the aforementioned collective networks, mainly in the northern and eastern ends of the study area [32].

Quaternary alluvial formations fill the boundaries along the southwest part of the basin with Neogene marls and sandy-clay deposits, whilst the western margins consist of

karstified marbles of the middle-upper Cretaceous. The crystalline bedrock is composed of mica-schists and gneisses of the upper Palaeozoic and Paleozoic age, respectively, forming the northern boundary of the basin (Figure 2). Two major springs (Mati Tirnavou and Agia Anna, Figure 2) emerge at the contact of the karstified system with the alluvial deposits. Pinios and Titarisios rivers flow across the basin, which, as already mentioned, hydrologically is part of the wider Pinios River basin.

The Quaternary deposits host an unconfined aquifer near the talus cone of Titarisios at the northwest, which, towards the central parts of the basin, sinks under a sequence of clay layers that form an aquitard [32–34] and have a maximum thickness of over 550 m at their central parts [32,35]. Under this setup, confining conditions occur, while the phreatic aquifer has been seriously affected due to systematic over-abstraction and is almost depleted [28,33]. The marbles at the western margins of the basin host a karstic aquifer of great potential, which recharges the alluvial system by lateral crossflows [28,33,34,36].

In the western part of the basin, an extensive marginal cone exists, through which the alluvial system receives significant amounts of recharge as crossflow through the Titarisios River gorge sediments. To the south, a smaller volume recharges the aquifer system as crossflow from the Pinios River gorge sediments [32].

Crossflow from the crystalline bedrock at the northern margins of the basin also occurs but is of minor importance. The southern extent of the karst system and from the mid-Thessaly hills recharge the central plain parts by crossflows from the southwest and the southern part of the area; however, they have a lower potential than that of the northern parts due to the existence of marls.

Based on the above data, it is perceived that the specific area is characterised by increased complexity, as reflected by the geological and hydrogeological conditions, tectonics, hydrodynamics, and land and water use. In addition, the area lacks regional planning and management of natural resources, leading to considerable deterioration of existing groundwater reserves. For this reason, RIVA was selected as an appropriate method of intrinsic vulnerability assessment to be implemented in the Tirnavos basin.

#### **3. Materials and Methods**

RIVA is the acronym of the four key factors considered in assessing intrinsic vulnerability, each of which has a distinguished control over groundwater vulnerability. The first factor is the recharge (R), which refers to the overall assessment of the effect of recharge conditions. The second one is the infiltration (I factor), which considers the soil's infiltration, subsequently affecting deep percolation to the saturated zone. The third factor relates to the protective cover and specifically to the impact of the vadose zone (V factor). The last factor is focused on the aquifer characteristics (A factor) and relates to the potential contaminant migration within the saturated zone as assessed by the hydraulic conductivity of the aquifer. Each of the four factors affects groundwater vulnerability individually and independently of the potential interactions. The final assessment of each factor is related to five classes of vulnerability (Figure 4), from very low (VL) to very high (VH), respectively [26]. The cumulative effect of these four factors constitutes the final assessment of intrinsic groundwater vulnerability, according to the relationship

$$\mathbf{i} = \mathbf{R} + \mathbf{I} + \mathbf{V} + \mathbf{A} \tag{1}$$

The contribution of each factor in Equation (1) is not equivalent due to the different degrees of importance. Hence, each one is multiplied by a weighting factor, which reflects its significance to the result. The weighting factors (a, b, c, d), which were calculated with the use of Analytic Hierarchy Process (AHP) and proposed by [26], are a = 0.40 for V factor, recognizing that is the most important, followed by the I factor with b = 0.30, c = 0.15 for R factor, and d = 0.15 for A factor (a + b + c + d = 1). V factor is the most important one in shaping the overall vulnerability, as it assesses the buffering capacity of the geological medium in constraining deep percolation of the potential pollutants to the saturated zone, the flux of which is directly related to their leaching rate potential as defined by factor I, which holds the second highest weight. Weighting factors' values are thus deduced as a fusion of the scientific approaches and documented viewpoints proposed by a number of researchers in the international literature [26].

**Figure 4.** Schematic presentation of RIVA method (according to [26]).

In Table 1, the data used for the calculation of each factor and sub-factor of the RIVA method as well as their source of origin are presented. The tables of each method's classes and ratings can be downloaded from the link at the end of this paper as a supplementary document. The analytical description of RIVA's rationale and conceptualization per factor is beyond the scope of this research and may be reached through the original application of the method [26].

**Table 1.** Aggregate data used for RIVA method implementation in Tirnavos basin.

#### **Table 1.** *Cont.*

#### **4. Results**

#### *4.1. R Factor Description*

R factor refers to the assessment of the internal vulnerability of the groundwater body as a function of the total surface recharge it receives and which can reach the aquifer through percolation under certain conditions. It includes three sub-factors [26], which are expressed through Equation (2):

$$\mathbf{R} = \mathbf{P} + \mathbf{In} + \mathbf{Ir} \tag{2}$$

where

R = recharge factor P = precipitation sub-factor In = rainfall intensity sub-factor Ir = irrigation recharge sub-factor

As in the case of Equation (1), the contribution of each factor in Equation (2) is not equivalent. The adopted weights are 0.5, 0.3, and 0.2 for P, Ir, and In, respectively, which are embedded in the intermediate calculations of each factor and are not shown in the initial equation. The cumulative effect of these three sub-factors constitutes the final assessment of the R-factor (Equation (2)), which is linked to the intrinsic vulnerability of the groundwater system according to the classification shown in Figure 4.

#### 4.1.1. P Sub-Factor Calculation

Sub-factor P refers to the total precipitation depth received at the examined area over a year (y), in millimetres (mm). The representative precipitation value is deduced as the mean annual value over a sufficient time, e.g., over 20 years. If more than one station is available, the final figure is calculated by spatially integrating precipitation distribution of individual stations over the examined area, employing an appropriate interpolation scheme. The corresponding P factor values concerning precipitation ranges, according to the authors [26], are as illustrated in Table S1 in the supplementary documents.

In the case study area, there are three mereological stations (Figure 2) that were taken into consideration in the P sub-factor calculations: Larissa and Tirnavos stations located within the examined area, at 74- and 92-m altitude, respectively, and Elassona station at 314-m altitude, being representative of the higher elevation parts included in the examined area. For each station, precipitation data for a period of 30 years (1989–2018) were used to obtain a mean annual precipitation value for each of them. The calculated values were 422.2, 485.7, and 520 mm for Larissa, Tirnavos, and Elassona stations. Based on class values tabulated in Table S1, a P value of 2.5 is assigned to all three stations, and therefore, no spatial distribution is required.

#### 4.1.2. In Sub-Factor Calculation

Except for the total precipitation, aquifer recharge depends on the intensity of the precipitation, which represents the amount of water received in a specific period. The effect of this parameter is difficult to estimate, as various factors, such as soil texture, root system, moisture state of the soil, the discontinuities of the geological formations, etc., form different conditions [50]. Nevertheless, it is generally accepted that an increase in rainfall intensity eventually increases the recharge and consequently the infiltration [51,52], thus affecting the vulnerability of the groundwater system. Regarding its calculation, the approach of the RIVA method is a compilation of approaches followed in similar methods [53,54]. The In values are calculated through Equation (3) [26]:

$$\mathbf{In} = \Sigma \mathbf{P} / \Sigma \mathbf{d} \tag{3}$$

where

In = intensity value of rainfall ΣP = total precipitation at a given period (mm) Σd = total number of rainfall days in the same period (days)

The rain intensity values for all stations were in the same range according to Table S2 in the supplementary documents (5.8 for Larissa, 8.9 for Tirnavos, and 9.5 for Elassona). Hence, the same In value (0.4) and vulnerability class (II) are assigned to the stations, and in this case, no spatial distribution method needs to be applied.

#### 4.1.3. Ir Sub-Factor Calculation

The Ir sub-factor represents the effect of irrigation on groundwater vulnerability and constitutes a parameter that has not been considered in vulnerability by any other method. This is possibly due to the difficulty of estimating the volume of irrigation water used and its spatial distribution, as several variables are included that are not easy to quantify. However, the importance of the irrigation effect in groundwater vulnerability is recognized, especially in intensively cultivated regions, such as the examined area. Hence, RIVA, considering the irrigation impact through a qualitative approach based on the appraisal of the mean overall performance, as shown in Table S3 (Supplementary document) [26], is deemed appropriate for assessing vulnerability in the examined area.

The irrigated fields at the case study area are spatially distinguished based on data from the Greek Payment Authority of Common Agricultural Policy Aid Schemes [39]. For the categorization of the irrigated fields in accordance to the classification referenced in Table S3, data retrieved from the Local Land Reclamation Organizations of the region were used along with assessments carried out during in situ visits to the area, accounting also for the irrigation methods employed in each part of the basin, in cases where no specific data were available. Thus, drip-irrigated fields were assumed to apply nominal irrigation doses and result in minimal leaching rates due to the high efficiency of the method. This is in contrast to sprinkler irrigation, where due to low water use efficiency, significant water losses occur that lead to considerable leaching rates, and therefore, increased irrigation doses are applied to satisfy crop water demands.

Based on this spatial data distribution according to the classes discussed (Table S3), the Ir sub-factor spatial distribution map was extracted, as illustrated in Figure 5.

#### *4.2. I Factor Description*

The I factor refers to assessing the intrinsic vulnerability of a groundwater body as a function of the surface infiltration conditions upon which deep percolation to the geological layers and eventually groundwater recharge depends. Its calculation is based on Equation (4), and its classification is shown in Table S1.

$$\mathbf{I} = \mathbf{F} + \mathbf{S} + \mathbf{C} \tag{4}$$

where

I = infiltration factor

In Equation (4), F and S sub-factors have equal weighting factors (0.5), while sub-factor C does not contribute equally. Sub-factor C refers to the concentrated flow that could cause increased infiltration due to specific surface structures, such as epikarst, river beds, and tectonic contact. It is considered a critical aspect, which may induce maximum vulnerability to the groundwater body if existing.

#### 4.2.1. F Sub-Factor Calculation

The F sub-factor relates to the surface water flow conditions, which affect infiltration to the saturated zone and consequently vulnerability of the groundwater system. Its calculation is determined by two parameters, the topographic slope (s) and vegetation (v). Based on the method's developers [26], the classification of "s" is made according to the approaches of [55–57], as shown in Table S4 in the supplementary documents.

**Figure 5.** Spatial distribution maps of P, In, and Ir sub-factors in the Tirnavos basin and the R factor map compiled as a synthesis of the sub-factors above.

The vegetation (v) parameter is used as a slope correction parameter, as the denser it is, the more it includes plants with a well-developed root system, the more the surface runoff is prevented, and the more the surface infiltration is favoured, creating conditions of higher vulnerability. Vegetation is classified into three general groups: (a) forest areas (high vegetation), (b) cultivated and grassland areas (low vegetation), and (c) absence of vegetation or sparse vegetation.

Classification of vegetation was made according to the widely acceptable CORINE categorization [29] and shown in Table S5 (Supplementary document), where only Categories 2 (agricultural areas) and 3 (forests and semi-natural areas) related to vegetation are considered relevant. Based on the method, for the rest of the CORINE categories, i.e., 1 (artificial areas-artificial surfaces), 4 (wetlands), and 5 (water bodies), it is assumed that no surface flow occurs; thus, by definition, the sub-factor F takes the value zero (0).

Based on older studies [58,59], the effect of the vegetation parameter (v) was assessed comparatively for all considered land covers; regarding forest vegetation (forest), the latter assumed to result in the least total soil loss and runoff. High-vegetation areas (forest) may cause a decrease of surface runoff coefficient up to 88%, while low-vegetation areas (pastures) cause up to 44% decrease, respectively. The result constitutes the final vulnerability class of the sub-factor F, which takes values from 1 (very low vulnerability— V.L.) to 5 (very high vulnerability—V.H.) (Table S6, Supplementary document). Therefore, the final value of the sub-factor F to be used in Equation (4) will be deduced empirically from the combination of slope parameter (s) and vegetation parameter (v), according to Table S6 in the supplementary documents.

Based on the slope (s) parameter distribution (Table S4) as derived by the digital elevation model of the study area [40] and the vegetation (v) spatial distribution in accordance to the classes discussed (Table S5), the F sub-factor spatial distribution values were calculated regarding Table S6, as illustrated in Figure 6.

#### 4.2.2. S Sub-Factor Calculation

The S sub-factor accounts for the permeability of surface geological formations and is directly proportional to the vulnerability of groundwater systems. The higher the permeability, the greater the vertical infiltration (percolation) and therefore the more significant recharge that will potentially reach the saturated zone of the groundwater system. RIVA regards surface formations that occur up to 1.5 m below surface and control surface/subsurface flow [26].

Surface formations are classified into soils and consolidated geological formations. Soils are considered the upper loose earth horizons, including topsoil, developed over nonconsolidated geological formations (e.g., Neogene formations, Quaternary deposits). As consolidated geological formations, RIVA considers those that are lithified and constitute the underlying bedrock. The method assumes that soils of considerable thickness are not developed over the consolidated geological formations. Therefore, the surface hydrological conditions are controlled mainly by the permeability of the consolidated formations and not by the soil compartment.

S values for soils are shown in Table S7 (Supplementary document) based on the U.S. Department of Agriculture classification that is based on their texture [60], whereas, for the geological formations, S values are deduced based on their permeability, as proposed by the British Geological Survey and presented in Table S8 (Supplementary document) [61]. Permeability is a property that is not easily determined accurately, as it is affected by several factors (e.g., degree of fracturing, tectonic stress, karstification, alternations with horizons of different permeability, homogeneity, isotropy, etc.), which may more or less change the original nature of the formation. However, guidelines for a more general characterization framework are not limited and provide a range of values, as shown in Table S8. The lower values correspond to solid (unaffected) formations, while the higher reflect the factors above effect that eventually increases their permeability and thus their vulnerability class. The final calculation of the S sub-factor and its spatial distribution map results from compiling

F values for the soils and F values for the consolidated geological formations. Note that for each cell of the generated map, a value S of soil or formation is assigned and not an aggregated value of two individual values of soil and formation.

**Figure 6.** Spatial distribution maps of F, S, and C sub-factors in the Tirnavos basin. According to the I factor, the final map of vulnerability classes results from the composition of the three above mentioned thematic maps.

Regarding the Tirnavos basin, Quaternary deposits based on the geological map (Figure 1) were classified based on Table S7 categories, considering the soil study of the area [41] in which soil units were assigned with soil and geological and topographic criteria. In addition, Neogene formations at southwest and northeast margins and marbles at the west and the crystalline bedrock, composed of mica-schists and gneisses at the northern boundary of the basin, are classified based on Table S8. In addition, taking into account the degree of karstification of the marbles [28,32,34], vulnerability class III is assigned to them, and class II is assigned to the Neogene marls (resulting from their permeability based on the same studies). The above data processing composed the distribution map of the S sub-factor, as shown in Figure 6.

#### 4.2.3. C Sub-Factor Calculation

Sub-factor C addresses the special cases of concentrated flow in fractured/discontinuities media that may impact the vulnerability characterization attributing the maximum vulnerability (very high) value to factor I. This sub-factor refers to the spatially concentrated because of specific surface features, which results in increased infiltration and thus maximum aquifer vulnerability. These features may include (a) epikarst; (b) drainage patterns, which are documented, that is, in hydraulic relation with the aquifer; (c) sinkholes; and (d) tectonic structures (e.g., faults, overthrusts). Infiltration is significantly favoured within the influence zone of the above features, causing eventually very high vulnerability.

As the exact orientation of the impact zones is not possible to be defined due to variable influencing factors that require complicated modelling approaches, RIVA proposes using an approximate impact zone of 100 m around critical surface features, to which a maximum score of 10 (V.H. vulnerability) is attributed [26]. If such structures do not occur (therefore no influence zones), the C value is by definition negligible. Hence, the C sub-factor constitutes an on-off (0–10) feature.

In the case of the Tirnavos basin, there are no indications of the existence of developed epikarst or sinkholes. However, some significant fault zones [42,43] can increase the infiltration and, consequently, the aquifer vulnerability (C values map, Figure 6). Furthermore, in earlier studies [28,32], it is suggested that in certain parts of the Pinios and Titarisios riverbeds, the hydraulic relationship between rivers and the groundwater system occurs. Regarding the Pinios River, confirmed hydraulic relation exists from its entrance to the Tirnavos basin up to Larissa city. In contrast, for the Titarisios River, significant hydraulic interaction occurs along the western margins of the basin, where extensive talus cones are formed (C values map, Figure 6).

#### *4.3. V Factor Description*

Factor V (vadose zone) accounts for the protection provided by the vadose zone as a function of the nature of its geological formations and its total thickness, which is directly related to the piezometric level. Factor V takes values from I (very low vulnerability) to V (very high vulnerability) according to the grading of Figure 4. As mentioned, it differentiates from the I factor because it regards the part below 1.5 m from the surface (upper soil horizons). In this context, the V factor may include (a) the soil's underlying non-lithified geological formations and/or strongly weathered zones of bedrock and (b) the bedrock (lithified geological formations).

Calculation of the V factor is performed through the modification and compilation of previous approaches [14,61–63] as follows [26]:

From the above methodology, the critical point is to estimate the thickness of the formation (or formations) that shape the vadose zone. For this purpose, the following approach was adopted. For the phreatic aquifer zone within the alluvial basin, the piezometric level distribution as deduced from field measurements was considered [28,34,36], along

with the compiled digital elevation model (DEM) (subtracting piezometric map contours from DEM). When only a confined aquifer is active, the thickness of the vadose zone is based on the bottom of the confining layer as mapped by older foundation studies in the region [44], subtracted from the DEM. For the karstic domain considered, the piezometric levels provided by older studies were considered [45], under the assumption that no significant changes have occurred in the spatial distribution of the groundwater levels in this environment. If they have occurred, they would have related to groundwater level decline, leading to an increased thickness vadose zone and consequently lower vulnerability assessment. Therefore, this assumption provides a rather conservative calculation, favouring the area's environmental protection. No piezometric data exist for the crystalline bedrock at the northern boundary of the basin. However, a low-capacity aquifer is reported to occur in the weathering mantle of the metamorphic formations, the thickness of which does not exceed 30 meters [28]. Therefore, it is reasonable to assume that groundwater level will not occur at a depth greater than that. Hence, this value is used as a baseline for the vadose zone thickness in this environment.

The compiled map for vadose zone thickness, based on the above data, was multiplied by the map synthesized as the product of ly and f (ly × f), based on Tables S9 and S10, and the geological map of the area, in order to deduce a protection cover (pc) map.

#### *4.4. A Factor Description*

The A factor (aquifer) refers to the easiness with which a potential contaminant will travel within the saturated zone of an aquifer as a function of its hydraulic conductivity. Factor A takes values from 1 for an aquifer of very low hydraulic conductivity (very low vulnerability—V.L.) to 5 for an aquifer of very high hydraulic conductivity (very high—V.H.) according to the classification in Figure 4.

The value of factor A is calculated based on the quantitative or qualitative assessment of the hydraulic conductivity of the aquifer. In the case of measured values, the link between hydraulic conductivity and vulnerability class (A factor values) is shown in Table S12 in the supplementary documents, based on a modified, generic, yet widely accepted conductivity classification according to [64].

For the alluvial basin, the values of the hydraulic conductivity of the geological formations were derived based on the results from the pumping test analyses [33,44] and the geo-hydrogeological characteristics of the aquifer. Moreover, the hydraulic conductivity values regarding the marbles at the west and metamorphic formations (gneisses, micaschists) at the N-NE were based on bibliographic references [46–49].

#### *4.5. Compilation of Results*

The calculations of factors, sub-factors' values, and spatial distribution were performed in ArcPRO® software to extract the digital maps, employing the corresponding equations and attributing class values as referenced in the mentioned Tables above. All data were transformed to raster grids with cell size 100 m × 100 m, which covered the case study area of the Tirnavos basin. The derived maps for each sub-factor, including the weighting factors and the cumulative map for main factors, are depicted in Figures 5–8.

Two vulnerability classes are assigned for the R factor: class I for the greater part of the study area and class II at the central and southeast part of the study area, as shown in Figure 5.

Regarding factor I, all vulnerability classes emerged in the spatial distribution map (Figure 6) following the raster calculations of its three sub-factors. The highest vulnerability class (V) was assigned along the traces of significant fault zones and surface-ground water system interaction zones. The class (IV) occupies the highest area percentage, and its main distribution is at the alluvial basin and the eastern part of the case study area, excluding its central part. The central part of the alluvial basin and most of the karstified domain are characterised by medium vulnerability class (III), while at the north and southeast, the low class (II) prevails. The lowest vulnerability class (I) has only isolated occurrences at the central and eastern parts of the area.

In the part west of Agia Sofia local community, low- and high-class alternations are presented, while the intermediate medium class is absent, as observed at the I factor distribution map. This fact is probably due to Neogene formations in the specific area and the southern part. However, what differentiates the two parts, is the topographic slope, which is lower in the south than in the west of the Agia Sofia village, where a slope break exists from the mountainous terrain to the alluvial basin. That justifies the appearance of the medium class in the southern part in contrast to the area west of the Agia Sofia village.

**Figure 7.** Spatial distribution of vulnerability classes according to the V factor in the study area.

In order to compile the V factor distribution map, based on Table S11, the pc values were attributed to a vulnerability class/characterization, corresponding to a V factor value. The map produced after this exercise is illustrated in Figure 7. In this map, the full range of vulnerability classes is also displayed.

Following the approach mentioned in the A factor description paragraph, a value of hydraulic conductivity correlated to vulnerability class was attributed to Table S12. Ultimately, a value of A factor was assigned to each formation and known permeability value point, and subsequently, data values were spatially interpolated with Inverse Distance Weighting (IDW) method to obtain the spatial distribution map for factor A (Figure 8).

**Figure 8.** Spatial distribution of vulnerability classes according to the A factor in the study area.

Since vulnerability maps were created for all four factors (R,I,V,A), groundwater intrinsic vulnerability (i) was estimated based on Equation (1), considering the corresponding weighting factors. The final calculations between raster images were performed following a normal grid 100 × 100 m cell distribution with ArcPRO GIS software. The value of any cell of the final map was derived by summing up the individual values of corresponding cells of all factor maps, according to Equation (1). The final intrinsic vulnerability map, which was produced following this procedure, is illustrated in Figure 9. Based on this map, the central and southeast parts of the case study area, where alluvial deposits dominate, are characterised by low vulnerability (II), occupying 48.75% of the total area. Comparing this map with the V factor distribution map, it is concluded that protecting the vadose zone in this area is highly effective since factor F is assigned the highest weighting factor. The karstic area at the west and northwest parts of the alluvial basin shows medium (III) to high (IV) vulnerability, occupying 24.04% and 14.14% of the examined domain. This fact reflects the sensitivity of the nature of the formations in these areas, as they are karstified marbles at the west side and conglomerates and talus cone at the northwest margins of the alluvial basin. The dominance of the high (IV) class at the west and north boundaries of the plain areas is reasonable because of the lithology and topography of the existing transition zones. The north and northeast boundaries of the study area appear to have the lowest vulnerability class (I), occupying 13.01% of the total area and reflecting the crystalline bedrock, which is mica-schists and gneisses, respectively. It should be noted that the areas of very high (V) vulnerability occupy a minimum percentage (0.07%) of the total area and are placed on the tectonic structures between Pinios and Titarisios Rivers.

#### *4.6. Validation*

The validation procedure is critical for the initial vulnerability assessment [27]. The most common approach, particularly for verification of assessments done with overlay and index methods, is to compare the vulnerability map with the actual occurrence of some common contaminant in groundwater [65,66]. Considering the initial vulnerability concept and the specific characteristics of the study area, validation was performed by comparing the monitored values of nitrates with the modelled vulnerability, as defined by the spatial distribution of the final vulnerability map [26] (Figure 9).

Nitrate values from 43 boreholes were used to create the spatial distribution map of nitrates in the Tirnavos basin, based on previous surveys and sampling campaigns [28,34,36]. Each value is the average result of four sampling campaigns at wet and dry seasons from September 2016–April 2018. These average values of nitrate concentrations were classified into ranges as follows: 0–20 mg/L, class I (very low); 20–40 mg/L, class II (low); 40–50 mg/L, class III (medium); 50–70 mg/L, class IV (high); >70 mg/L, class V (very high). Based on the above classification, which constitutes a modified version of the groundwater nitrates concentration classification proposed in the framework of the Nitrates Directive [67], a spatial distribution map of nitrates' classes was compiled (Figure 10). Nitrates' classes distribution map was subsequently compared to the final vulnerability map, a result of summing RIVA factors (Figure 9). The subtraction of these two maps yields (modelled values—monitored values) the validation map (Figure 11).

**Figure 9.** Spatial distribution of vulnerability classes according to RIVA method in the study area.

**Figure 10.** Spatial distribution map of nitrates concentrations (average values of four periods, September 2016–April 2018).

**Figure 11.** Validation map of the RIVA application method in the Tirnavos basin. The class difference between modelled and monitored values, depicted.

In Table 2, the analytical results of the validation procedure are shown. Based on this table, 29.14% of the total area presents a perfect match between modelled and monitored

− −

−

−

−

values (difference class 0). Overall negative difference (underestimating, lower modelled than monitored values) presents 13.20% and positive difference (overestimating, higher modelled than monitored values) 57.65%. According to the original validation of the method [26], a difference of one class (−1 to 1) is accepted as a very good match between the modelled and monitored values. If that percentage is equal to or greater than 80%, then the modelled results are successful compared to the ground truth values.

**Table 2.** Difference of vulnerability class between modelled and monitored values.

According to the validation performed in the present research, the total area with a difference from −1 to +1 between the modelled and monitored values is 81.52%; thus, the validation may be regarded as successful. Continuing, 96.95% of the total area presents a difference of two classes (−2 to 2), and only 3.05% exhibits a difference of three classes (−3 or +3).

The largest class differences (−3) where RIVA underestimates the potential vulnerability appear in a few individual sites related either to point-source contamination or to the effect of migrating contaminants through lateral crossflows. The latter is probably the case for the deviations in the alluvial basin related to zones of lateral contaminant fluxes and especially so along the southwestern-most edge of the basin. By definition, aquifer vulnerability accounts for the vertical susceptibility of the system and does not account for any lateral crossflows of contaminant plumes from adjacent hydrogeological units. On the other hand, the infiltration of N-containing pollutants from surface water and the transport of nitrate contaminants through soil and groundwater occurs via a series of complex chemical and hydraulic phenomena [68]. As a result of these complex procedures, a horizontal migration of the contaminants is frequently dominant. Hence, practical validation with measured contaminant values at the saturated zone should only be performed at hydrologically "closed" systems without any hydraulic connections with other units (surface or underground); if not, then potential migrations of contaminant plume(s) must be taken into account prior to validation [69].

Indeed, based on the geometry of the groundwater system and its hydrodynamic evolution, as these were discussed in the earlier parts of this work, the area where the highest deviation between modelled and monitored values coincides with the parts of the system where contaminants' migration is influenced by a conurbation of factors where not only percolation is important, but lateral flows do play a significant role. Moreover, it needs to be taken into account that the results of water sampling on nitrate concentrations do not always reflect the regional background regarding contaminant's state, but instead, they may indicate local factors of N-bearing compounds mismanagement. Hence, deviations between the RIVA and nitrates distribution classes are sufficiently explained.

#### **5. Discussion**

#### **Insight to the factors used**

The effect of precipitation and irrigation on groundwater vulnerability was examined by calculating the R factor in the study area. The mean annual approach is more efficient for strategic planning at a regional scale, whilst the seasonal approach can flag vulnerability aspects that are significant for local or regional scale risk assessment uses [70,71]. Although data from three meteorological stations located at different altitudes were used, both the mean precipitation (P sub-factor) and its intensity (In sub-factor) did not form separate vulnerability classes based on the various ranges that have been defined.

Irrigation (Ir sub-factor) as a factor of vulnerability influence has not been considered in any other method yet. In rural areas where agriculture forms the key socio-economic activity, irrigation is indeed an overlooked parameter of paramount importance in assessing the vulnerability of a groundwater system, as clearly demonstrated through the performed analysis. Especially in the absence of multiple classes for the other 2 sub-factors shaping R factor (P + In), the irrigation sub-factor considerably influenced the spatial distribution of the R factor. Given that the average field size at the study area is rather small (ca 1 ha), it is thought that should irrigation system data occur at the field scale, a more detailed spatially diversified distribution of the R factor would have been yielded. In turn, such an approach would provide a higher spatial resolution assessment of the system's vulnerability concerning this particular sub-factor. On a regional scale, however, this differentiation would have not led to considerably different results apart from the potential for even better agreement between the vulnerability map and the ground-truth nitrate concentrations and, in particular, the ability to better explain hot spots of the nitrates concentration distribution map. However, in the absence of such detailed data, the dominant irrigation system per local irrigation organisation inherently incorporates the assumption of the least waterefficient irrigation system used. Hence, the higher irrigation water is applied on the field, which subsequently leads to a relative over-estimation of water use and, therefore, higher assessments of vulnerability due to this sub-factor. Overall, the approach followed results in a more conservative estimate of vulnerability, i.e., a more environmentally sound approach.

Assessing I factor incorporates evaluating the combined action of slope and vegetation, which constitutes a straightforward procedure based on the principle that the lower the slope and the denser the vegetation, the greater the class of attributed vulnerability. It also includes evaluating permeability of the top 1.5 m of soils and consolidated geological formations, as assigned with the relevant tables developed by the method. Last, it incorporates the effect of linear features, including tectonic lines [72,73] and hydrologically interacting river stretches [74,75]. As a whole, the subjectivity of this factor is rather limited due to the wide classes distinguished for the evaluation of each sub-factor and the overall sound principles it adopts in the definition of the considered sub-factors.

As already stressed, V factor constitutes the most influential factor in assessing the examined system's vulnerability, and this is acknowledged by assigning to the final vulnerability index calculation formula a very high weighting factor (a = 0.4). Hence, it is rather critical to evaluate this factor carefully and based on reliable information to avoid misconceptions on the final spatially distributed vulnerability product. Out of the three considered sub-factors, the dominant lithological type (ly) and the corresponding fracturing or karstification degree (f), which represent the bulk permeability values per formation and the alterations imposed due to secondary deformations, are safely assessed utilizing bibliographic references [28,61,63], geological maps [30,31], field observations [28,34,36], and the adopted classification scheme of RIVA, as earlier presented [26]. Calculation of the protective cap, which forms the third sub-factor (pc), is the most critical and presents the highest risk for erroneous calculations, especially when groundwater system vulnerability assessments are attempted for complex geometry environments, as the examined one is. However, the adopted approach and utilization of the detailed data available on the system's geometry for most of its parts ensures that performed assessments are reliable and, in fact, shifted towards the conservative margin of analysis. Thus, calculations tend to

accept the minimum thickness of protecting cap in the few parts of the basin where relevant data for this sub-factor are not detailed. As with the I factor, this provides a safety margin in favour of environmental protection, assuming the worst-case scenario on the geometry of the considered system.

Last, factor A reflects the influence of hydraulic conductivity on intrinsic vulnerability, and its spatial distribution is deduced based on pumping test analyses, augmented by region-specific and more generic but well-established bibliographic references [46–49] transposed to wide range classifications, as earlier presented. Hence, the margin of error in assessing the spatial distribution of this factor is limited.

Validation of the method was performed on the basis of a considerable population of real monitoring data that cover dry and wet hydrological conditions for two years and are typical of an agriculture-related pollutant. Considering the particularities of the studied region as it is reflected by geological structure, hydrogeological setup, and regional hydrodynamic evolution mechanisms, validation results are satisfactory, demonstrating the validity and efficacy of RIVA in similar complex environments. Spatially distributed zones of different vulnerability classes are delineated accurately and reliably even though the tested hydrogeological environment is characterised by a high degree of complexity.

Observed deviations relate mainly to intrinsic vulnerability assessment methods' inherent limitations in accounting for horizontal flow driven migration of pollutants. Moreover, it is expected that the produced results would have been further improved if monitored data utilized for validation were originating from a single aquifer of the examined system rather than reflecting an average concentration of the entire system. That is the result of the construction characteristics of most production wells in the studied region, which tap the entire groundwater system that consists of multiple aquifers, the hydrochemical evolution of each one of which is controlled by different mechanisms. Last, the highest deviations noted reflect hot spots areas that are most likely related to one or more of the following reasons: (a) N-compounds mismanagement at the vicinity of the well heads, (b) systematic over-irrigation leading to excessive leaching of contaminants, and (c) poor construction characteristics of production wells that enable migration of contaminants to deeper hydrogeological strata through the wells' gravelpack.

In total, implementation of RIVA in the studied region provides a valuable tool for regional planning and management of natural resources that the area lacks, leading to considerable deterioration of existing groundwater reserves. Compared to other approaches, the selected method enables uniform consideration and evaluation of the entire region structured of highly contrasting character lithological typologies, rendering vulnerability assessment easier and comparable amongst the existing geological media. Of particular importance is the consideration of irrigation water management, which leads to a more comprehensive assessment of the essence in rural areas where intensified agriculture is practised. Even though irrigation is an externality, being an established condition, it directly influences vulnerability and should be given a thorough, reliable, and pragmatic assessment. In several parts of the world, data scarcity is a fact. Still, reliable vulnerability assessments are imperative and may not await appropriate state-of-the-art data collection. RIVA is not a data-intensive method that can easily retrieve information based on sound geological and hydrogeological knowledge. Obviously, enhanced results may be obtained, and a higher spatial resolution achieved should spatially dense and reliable data-driven calculations be employed.

#### **Challenges and problems**

Groundwater systems vulnerability assessment is a valuable tool with profound applications in environmental science. One of them could potentially be its incorporation, as a factor in the delineation of Nitrate Vulnerable Zones (NVZs), in the framework of the Nitrates Directive [67]. Thessaly as a whole is one of the first regions of the country to be declared as an NVZ in the late 1990s [76], using as key criteria the surface and groundwater concentration of nitrates, the N-input, and the geomorphological characteristics of the region. Groundwater system vulnerability consideration would add to the reliability of the

NVZ assessment and perhaps even lead to a considerably differentiated delineation pattern, accounting for the essential properties that control the potential leaching of contaminants. Therefore, such an approach would lead to a more comprehensive and meaningful delineation of NVZs and greatly assist the design of appropriate and targeted, thus efficient, measures as part of the Action Plan in the framework of the Nitrates Directive. Such an approach is directly in line with the Common Agricultural Policy, which identifies the direct relationship between agriculture and groundwater resources and the need for effective protection and preservation of the latter from agricultural inputs [77,78].

Spatially distributed flow models have been developed for the studied area and the wider Thessaly region [79–81], which provide a reliable simulation of the flow domain accounting for the key hydrodynamic evolution mechanisms. Contamination transport and even hydrogeochemical modelling coupled with such flow models would enhance understanding of evolution mechanisms and improve validation of the results. Advective and dispersive migration elements of the pollutants could be quantified and isolated, thus enabling backwards analysis of nitrates concentration distribution deduced from the monitoring exercise, accounting for reducing or oxidizing reactions that influence groundwater chemistry and thus the nitrates concentrations in the collected and analysed water samples. In this way, leaching would be the sole pathway to be reflected on the nitrates distribution, thus directly relating to the evaluation provided by RIVA, which assesses vulnerability considering the vertical pathways of potential pollutants. On a different viewpoint, such models could be coupled to RIVA to account for the advective and dispersive migration of a potential pollutant, thus contributing to a more holistic assessment of groundwater systems' vulnerability. In such an approach, the vulnerability would be calculated as a function of the examined parameters in the framework of RIVA (for leaching potential assessments) and lateral crossflows from surrounding formations or groundwater systems.

As pointed out in previous parts of the paper, seasonal variations of precipitation patterns lead to differentiated calculations of key factors of the vulnerability index. Likewise, consideration of seasonal groundwater levels as opposed to mean annual or interannual levels may lead to considerable differences in the calculations of the protective cap sub-factor of V factor, which has the highest weighting factor in the calculation of the vulnerability index. Previous studies have demonstrated this in similar geo-climatic environments [82,83]. In the framework of climate change, considerable variations of key climate parameters, including precipitation depth and intensity, temperature, and evapotranspiration, are anticipated, especially for the region Thessaly, which is considered amongst the vulnerable regions along with the entire Mediterranean [84–86]. These changes will potentially affect the hydrological balance of the system in terms of anticipated natural recharge and abstractions given that irrigation water needs are expected to increase under the business as usual scenario. In turn, depth to groundwater is expected to be affected. Prolonged droughts duration and frequency along with the already experienced severe floods because of high rainfall intensity directly differentiate directly the values of several sub-factors considered in the vulnerability assessment under the RIVA method.

Likewise, improving irrigation water efficiency is also expected in response to climate change and resilience development. Finally, crop distribution changes are already being discussed as an adaptation measure. These are some of the key changes expected in the future that will influence the vulnerability of the considered groundwater system. The factors above will alter the spatial distribution of vulnerability. However, some of them will increase it, whilst others decrease it. Therefore, the overall prediction of the anticipated outcome may not be safely made empirically. Still, reliable assessments of the forecasted vulnerability under climate change environment are important to augment measures designed to increase groundwater systems' resilience and overall water sector safety.

#### **6. Conclusions**

This study assessed the intrinsic vulnerability of the Tirnavos basin (Central Greece) groundwater system using the RIVA method. RIVA has the advantage that it can be applied to all types of groundwater bodies independently of the specific conditions, lithologic phases, and aquifer typology of each area. Four main factors were used to represent the natural hydrogeological conditions of the specific area: recharge, infiltration, vadose zone, and aquifer.

Modelled results were validated with ground-truth values of nitrates obtained from 43 wells and proved to be quite successful, as they presented over 80% of matching (negligible or small deviations between modelled and monitored values). The few deviations are attributed to inherent uncertainty factors, such as the interpolation of the various factors used and the lateral contamination transport, which the index and overlay methods cannot assess without the help of a spatially distributed model.

The outcomes of the RIVA application to the Tirnavos basin can be further exploited as a preliminary tool for decision making and strategic assessment. At a regional scale, the results may be further valorised for the delineation of NVZs, considered as a key target for the trade-off between sustainable agriculture and water resources protection.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/w14040534/s1, Table S1: Correlation between annual precipitation (mm/y), P values, and vulnerability classes; Table S2: Correlation between rainfall intensity (mm/d), In values, and vulnerability classes; Table S3: Correlation between irrigation dose, In values, and vulnerability class; Table S4: Classification of topographic slope (s) and correlation with vulnerability classes; Table S5: Categorization of vegetation according to CORINE land-use classification (EEA, 2020); Table S6: Calculation of F values according to soil (s) and vegetation (v) parameters and link with vulnerability classes; Table S7: Calculation of S values for the soils and correlation with the vulnerability class (USDA, 1999); Table S8: Calculation of S values for consolidated geological formations and correlation with permeability and the vulnerability class (Lewis et al., 2006); Table S9: Classification of layer reference values (ly) for representative geological formations; Table S10: "f" factor values according to the assessed fracturing or karstification degree; Table S11: Protective cover (pc) values and corresponding vulnerability characterization and V factor values; Table S12: The suggested link between hydraulic conductivity and A factor values, corresponding to specific vulnerability class.

**Author Contributions:** I.V. conceived the methodological approach, performed the data processing and supervised the drafting and revision of the final text. E.T., A.P. and G.S. participated pin data processing, verified and optimised the methodology, contributed to the discussion parts, and revised the final text. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** Data is contained within the article.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Hydrogeochemical Processes and Natural Background Levels of Chromium in an Ultramafic Environment. The Case Study of Vermio Mountain, Western Macedonia, Greece**

**Eleni Vasileiou 1 , Panagiotis Papazotos 1 , Dimitrios Dimitrakopoulos 2 and Maria Perraki 1, \***

**Abstract:** The hydrogeochemical processes and natural background levels (NBLs) of chromium in the ultramafic environment of Vermio Mountain, Western Macedonia, Greece, were studied. Seventy groundwater samples were collected from 15 natural springs between 2014–2020, and an extensive set of physical and chemical parameters were determined. The ultramafic-dominated environment of western Vermio Mt. favors elevated groundwater concentrations of dissolved magnesium (Mg 2+ ), silicon (Si), nickel (Ni), and Cr in natural spring waters. Chromium was the principal environmental parameter that exhibited a wide range of concentrations, from 0.5 to 131.5 µg/L, systematically exceeding the permissible limit of 50 µg/L for drinking water. Statistical evaluation of hydrogeological, hydrochemical, and hydrological data highlighted the water-ultramafic rock process as the predominant contributor of Cr in groundwater. The NBL assessment for Cr and Cr(VI) was successfully applied to the typical ultramafic-dominated spring "Potistis" that satisfied all the methodology criteria. The NBLs of Cr and Cr(VI) were defined at 130 µg/L and 100 µg/L, respectively, revealing that a natural ultramafic-dominated environment exhibits the geochemical potential to contribute very high concentrations of geogenic Cr to groundwater. The holistic methodology, proposed herein, could be implemented in any catchment scale to assess geogenic and anthropogenic Cr-sources that degrade groundwater quality.

**Keywords:** chromium; ultramafic rocks; springs; water–rock interaction; natural background levels

#### **1. Introduction**

Natural background levels (NBLs) are defined as "the concentration of a substance or the value of an indicator in a body of groundwater corresponding to no, or only very minor, anthropogenic alterations to undisturbed conditions" according to the Groundwater Daughter Directive (GDD) (Directive 2006/118/EC) [1]. Broadly, the term of NBLs is synonymous with the terms of environmental geochemistry "natural/geochemical background values" or "geochemical baseline" used in the past [2]. The NBLs are a set of several varying hydrogeological (i.e., the residence time of groundwater in the saturated zone, recharge by precipitation, hydraulic connection with other aquifer systems) [3–5], and hydrogeochemical (i.e., water–rock interaction, pH/redox conditions, chemical, and biological processes in the unsaturated zone) [5–7] factors. The determination of NBLs requires in-depth knowledge of geological/hydrogeochemical processes [8] and the distinguishment of natural and anthropogenic factors that affect the groundwater systems [9]. The need to separate NBLs from the anthropogenic impacts (e.g., urbanization, industrialization, agricultural activity) is frequently satisfied through statistical and pre-selection (PS) methods [10]. Such methods were applied within the EU-Specific Targeted Research Project BRIDGE (Background cRiteria for the iDentification of Groundwater thrEsholds),

**Citation:** Vasileiou, E.; Papazotos, P.; Dimitrakopoulos, D.; Perraki, M. Hydrogeochemical Processes and Natural Background Levels of Chromium in an Ultramafic Environment. The Case Study of Vermio Mountain, Western Macedonia, Greece. *Water* **2021**, *13*, 2809. https://doi.org/10.3390/ w13202809

Academic Editors: Evangelos Tziritis and Andreas Panagopoulos

Received: 2 August 2021 Accepted: 2 October 2021 Published: 9 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the objective of which was to develop a comprehensive methodology to evaluate threshold values (TVs) and NBLs of various qualitative parameters in the groundwater resources [11]. The first stage of this approach includes the PS method, which assumes that the groundwater samples represent pristine groundwater not affected by anthropogenic pressures [9]. It constitutes the most frequent method to exclude samples influenced by anthropogenic activities based on specific criteria such as concentrations of Cl−, Na+ , NO3 −, NH4 + , and DO [6,10,12]. The PS method has been successfully applied to establish NBLs for different physical and chemical parameters, including EC, Cl−, SO4 2−, F−, As, Cr, Cr(VI), Mn, Ni, Fe, and V in many European water bodies [4,5,7,13–17]. The next stage contains statistical tools such as box plots and normality tests for assessing the NBLs of the target chemical parameter. An approach that incorporates both PS and statistical methods has been performed by many researchers [7–10,14,15], providing a comprehensive methodology to boost the validity of the assessment, mainly when the geochemical and geological features are adequately considered [16]. Thus, the challenging assessment of NBLs in an environment in which the prevailing geochemical conditions favor the occurrence and mobilization of naturally occurring chemical elements could provide essential information regarding the controversial geogenic and anthropogenic inputs in a complex environmental setting.

The water–ultramafic rock interaction is of great scientific interest due to the high content of the latter in Cr (1000–3000 mg/kg), and other potentially toxic elements (PTEs) such as As, Co, Fe, Mn, and Ni compared to the Earth's crust composition [8–21] and to other rock types [22]; it constitutes the principal geogenic source of Cr in the environment [18]. Chromium is mainly hosted in spinels (e.g., chromite and magnetite) and silicates (e.g., pyroxene, serpentine, chlorite, olivine, talc). Serpentine group minerals can be highly enriched in Cr because it substitutes for magnesium (Mg) and/or iron (Fe) [18]. In the crystal lattice of most minerals, Cr occurs in the trivalent valence state [Cr(III)]. However, the geochemically immobile Cr(III) is oxidized into the mobile and toxic for the living organisms hexavalent chromium [Cr(VI)] in the presence of natural manganese oxides (MnO2), specifically pyrolusite (b-MnO2), in the typical range of groundwater pH (6.5–8.5) and under oxidizing redox potential (Eh) conditions [23–28]. Although an increasing number of studies focus on the occurrence and fate of Cr in the environment [29–38] only a few have systematically examined the geochemical fingerprint of water–ultramafic rock interaction in natural springs [29,33,35,39]. Typical worldwide examples of ultramafic springs with elevated groundwater concentrations of Cr(VI) have been recorded in the Province of La Spezia, Italy (up to 73 µg/L) [29], the Pollino massif, Italy (up to 30 µg/L) [40], the Gerania springs, Greece (up to 17.2 µg/L) [33], the Euboea Island, Greece (up to 37 µg/L) [41], and the Lesvos and Rhodes Islands, Greece (10–15 µg/L) [42]. The water–rock interaction constitutes a crucial and controlling factor concerning groundwater evolution. The geochemical reactions between the recharging water and the minerals of the host rocks affect the groundwater quality [43,44]. Hydrogeological and hydrogeochemical conditions such as pH, Eh, dissolved oxygen (DO), and groundwater flow path play a significant role in elevated concentrations of PTEs, including Cr, in the aquifer systems. Geochemical reactions such as ion exchange, weathering, precipitation/dissolution, and sorption process control the groundwater's composition considerably. During chemical weathering, some major ions, PTEs, and other trace elements become mobile and release from the parent rocks to the groundwater along the flow path. In addition, the mobility and solubility of these elements are controlled by water–rock contact time, Eh–pH conditions, and chemical reactions with organic matter [36]. Ionic ratios, saturation indices (SIs), and geochemical bivariate plots are usually evaluated to determine the intensity of water–rock interaction and chemical reactions [43]. Hence, the primary target of studying the mechanism of water–rock interaction is to elucidate the indissoluble association between the geological environment and the qualitative characteristics of groundwater.

In this work, we study the geochemical fingerprint of the water–ultramafic rock interaction process in the western Vermio Mt., Western Macedonia, Greece, and determine the NBLs of Cr in groundwater from natural springs. At the catchment scale of the

Sarigkiol Basin, elevated groundwater concentrations of Cr (up to ~140 µg/L) have been recorded in irrigation wells in the lowland [45]. Based on geospatial and multivariate statistical analyses of data from selected natural springs, irrigation wells, and surface waters the increased concentrations of Cr were attributed mainly to geogenic origin with the synergistic contribution of anthropogenic factors [45]. Challenged by the leaching potential of Cr of the ultramafic rocks in the area, we focus, herein, exclusively on the natural springs located in the ultramafic environment of western Vermio Mt, assessing hydrogeochemical data of a 7-year monitoring period (2014–2020). The springs are ideal for setting the NBLs at the catchment scale of the Sarigkiol Basin, because: (a) they record a strong ultramafic footprint, (b) they are located at a high altitude (>1300 m), (c) they exhibit unique worldwide high to very high concentrations of Cr (up to ~130 µg/L) [45], and (d) they are not affected by anthropogenic activities. Defining the NBLs of Cr in western Vermio Mt., will facilitate the identification of Cr origin in groundwater in the Sarigkiol Basin. This is the first systematic study of the natural springs of western Vermio Mt. and provides important hydrogeochemical data for the geogenic footprint of a natural ultramafic environment on the groundwater quality.

#### **2. Materials and Methods**

This section contains basic information about the: (i) study area, (ii) geological and hydrogeological setting, and (iii) sampling, chemical analysis procedures, and data processing.

#### *2.1. Case Study*

The present study is focused on the western Vermio Mt., located in the eastern part of the Sarigkiol Basin, Western Macedonia, Greece. The altitude of Vermio Mt is 2025 m and the average altitude of the basin is 650 m. The study area lies between the latitudes 40◦25′00′′ and 40◦28′00′′ E and the longitudes 21◦56′00′′ and 22◦59′00′′ N (Figure 1). In this area, any extensive anthropogenic activities lack except for local livestock farming and sporadically logging.

#### *2.2. Geological and Hydrogeological Setting*

The western Vermio Mt. is composed of (Figures 1 and 2) [46,47]: (a) alluvial deposits, (b) clastic conglomerates, talus cones, and breccias, (c) upper Cretaceous flysch, (d) a complex of schists and cherts formations, (e) ultramafic rocks (serpentinites and peridotites), (f) Triassic–Jurassic limestones, (g) Cretaceous limestones.

The aquifer systems in Vermio Mt. are:

Different types of natural springs flow out in Vermio Mt. and specifically (Figure 1):

Contact springs formed where permeable formations (limestones, breccia, conglomerates) overlay formations of low permeability or impermeable (altered ultramafic rocks/ serpentinites). Contact springs studied here were: the springs S19, S10, S15 in the Agio Pnevma area, the "Potistis"-W13 spring, the spring S1 in the Agios Dimitrios area, the spring "Mouratidis"-S2, the springs S13, S14 in the Agios Panteleimonas area, and the springs S5 and S6 in the Vazelona area.

**Figure 1.** A simplified geological map of the western Vermio Mt.; the natural springs studied herein are marked.

‒

**Figure 2.** A simplified geological section of the western Vermio Mt.

Fault springs formed where impermeable rocks such as ultramafic rocks are in contact with an unconfined aquifer due to faulting; the springs "Elafakia"-W14 and W21 in the Agios Panteleimonas area belong to this type.

‒ ‒ ‒ ‒ ‒ The spring "Potistis"-W13, presents great interest, because of the very high concentrations of Cr it exhibits [45]. The spring "Potistis"-W13 flows out at an elevation of 1300 m in an ultramafic environment characterized by the absence of any anthropogenic activities. It constitutes a contact-type spring in the contact of conglomerates with ultramafic clastic material and limestones and ultramafic rocks. The aquifer, which discharges via the spring, flows through a weathered zone in serpentines. The natural recharge comes mainly from the seasonal precipitations via the permeable upper unsaturated zone (conglomerates, clastic material of ultramafic rocks, and limestones). The recharge water is mainly enriched with released PTEs (mainly Cr) from the ultramafic rocks, as the rainfall infiltrates through the weathered fractured ultramafic rocks. An additional lateral recharge takes place due to secondary porosity in the fractured ultramafic rocks. The high permeability of the unsaturated zone due to the presence of conglomerates in this area facilitates the direct recharge of the aquifer in a short time. The range of discharge was calculated from 205 L/h up to 1200 L/h, with an average value of 482 L/h. In Figure 3, the simplified hydrogeological section describes the natural recharge and the operation mechanism of the spring "Potistis"-W13.

In western Vermio Mt. ultramafic rocks, mainly serpentinites, carbonates, schists and cherts occur [48]. The main mineral phases of the ultramafic rocks, depending on the degree of serpentinization, are serpentine [(Mg, Mn, Fe, Co, Ni)3−xSiO2O5(OH)4], olivine [(Mg,Fe2+)2(SiO4)]), pyroxene [(Mg,Fe2+)(Si,Al)2O6], talc [Mg3Si4O10(OH)2], chlorite [(Mg,Fe2+)5Al(Si3Al)O10(OH)8], tremolite [Ca2(Mg,Fe2+)5Si8(OH)2O22(OH)2], magnetite (Fe2+Fe2 3+O4) and Cr-rich magnetite [Fe2+(Fe3+,Cr)2O4] and chromite (FeCr2O4).

‒

‒ **Figure 3.** Simplified hydrogeological section of the natural spring "Potistis"-W13 in western Vermio Mt.

#### *2.3. Sampling, Chemical Analyses, and Data Treatment*

− ‒ Springs to be systematically studied herein were selected based on their hydrogeochemical characteristics [45]. A total of 70 representative groundwater samples were collected from 15 natural springs during wet and dry periods from March 2014 to September 2020, following the groundwater sampling guidelines [49]. The 15 sampling sites (Figure 1) were classified into seven groups according to their location, lithology, and type: (i) twenty-three (23) samples were collected from the spring "Potistis"-W13, (ii) seven (7) samples from the Agio Pnevma area (S18, S16, S19, S10, S17, S15), (iii) twenty-three (23) samples from the spring "Elafakia"-W14, (iv) three (3) samples from the Agios Panteleimonas area (S13, S14, W21), (v) ten (10) samples from the spring "Mouratidis"-S2 (vi) two (2) samples (S5 and S6) from the Vazelona area, and (vii) two (2) samples from the spring S1 in the Agios Dimitrios area. Considering that the number of the water samples differs between the seven groups, each group has been treated and evaluated separately (the statistical and geochemical analysis), so the analyses are classified as reliable.

‒ ‒ The analytical methods for the determination of physical [i.e., temperature (T), pH, oxidation-reduction potential (ORP), DO, and electrical conductivity (EC)] and chemical parameters (i.e., major ions, PTEs, and other trace elements) are provided in detail in Papazotos et al. [33]. The calculations of Eh and the total dissolved solids (TDS) values were carried out by converting ORP measurements (i.e., adding 200 mV) and the summation of major ions in each collected water sample, respectively.

AquaChem 5.0 software was used to elaborate chemical analyses, develop a Piper diagram, and calculate alkalinity. The statistical analyses of the chemical data were performed with SPSS 22.0 software.

#### *2.4. Spearman's Rank Correlation Coefficient*

Spearman's rank correlation coefficient (ρ, also signified by rs) measures the strength and direction of association between two ranked variables, evaluating the degree of linear association or correlation between these independent variables. It presents many similarities to Pearson's coefficient except that it operates on the ranks of the data rather than the raw data [50].

The Spearman's rank correlation coefficient is calculated according to the following Equation (1) [51]:

$$\mathbf{r}\_{\mathbf{s}} = \frac{6\sum\_{\mathbf{i}=1}^{\mathbf{n}} \mathbf{d}\_{\mathbf{i}}^{2}}{\mathbf{n}(\mathbf{n}^{2} - 1)}\tag{1}$$

where di = difference in paired ranks, n = number of cases, xi and yi = data pair. The formula to use when there are tied ranks is Equation (2):

$$\rho = \frac{\sum\_{\mathbf{i}} \left(\mathbf{x\_i} - \bar{\mathbf{x}}\right) \left(\mathbf{y\_i} - \bar{\mathbf{y}}\right)}{\sqrt{\sum\_{\mathbf{i}} \left(\mathbf{x\_i} - \bar{\mathbf{x}}\right)^2 \sum\_{\mathbf{i}} \left(\mathbf{y\_i} - \bar{\mathbf{y}}\right)^2}} \tag{2}$$

The Spearman's rank correlation coefficient, rs, can get values from −1 to +1. The equation for the calculation is developed so that it gives rs = +1 when the data pairs have a perfect positive correlation (di = 0) and rs = −1 for the perfect negative correlation, whereas rs = 0 indicates no association between ranks. The closer rs is to zero the weaker the association between the ranks is.

The values of the correlation coefficient are classified as very strong (0.80–1), strong (0.60–0.79), moderate (0.40–0.59), weak (0.20–0.39), and very weak (0.00–0.19) [52]. The correlation coefficient is highly statistically significant, marginally statistically significant when the *p*-value is *p* < 0.01, *p* < 0.05.

#### *2.5. Shapiro-Wilks Test*

Shapiro–Wilks is a test of normality in frequentist statistics. The null hypothesis of this test is that the dataset is normally distributed. Thus, if the *p*-value is less than the chosen alpha level (0.05 in this case), then the null hypothesis is rejected and the data tested are not normally distributed. If the *p*-value is greater than the selected alpha level, then the null hypothesis cannot be rejected (Equation (3)) [53].

$$\mathbf{W} = \frac{\sum\_{\mathbf{i}=1}^{n} \left(\mathbf{a\_{i}x\_{(i)}}\right)^{2}}{\sum\_{\mathbf{i}=1}^{n} \left(\mathbf{x\_{i}} - \bar{\mathbf{x}}\right)^{2}} \tag{3}$$

where x(i) is the i-th largest order statistic, x- is the sample mean, and n is the number of observations.

#### *2.6. Quantile–Quantile Plot*

The quantile-quantile (q–q) plot is a graphical tool for defining if two datasets come from populations with a common distribution [54], basically tests the conformity between the empirical distribution and the given theoretical one. On a Q–Q plot normally distributed data, the points in a Q–Q plot will fit on a straight diagonal line.

#### *2.7. Geochemical modeling*

The geochemical software PHREEQC version 3.1.2 [55] coupling with the MINTEQ database was used to calculate the saturation indices (SIs) of natural spring samples. Mineral SIs employed to define mineral dissolution and precipitation processes in the natural springs of western Vermio Mt. Saturation index is calculated by the Equation (4):

$$\text{SI} = \text{Log}\frac{\text{IAP}}{\text{K}\_{\text{sp}}} = \text{Log}\text{IAP} - \text{Log}\text{K}\_{\text{sp}} \tag{4}$$

where IAP = ion activity and Ksp = solubility product constant.

A positive SI indicates that the mineral is oversaturated or supersaturated with respect to the solution [56]; thus, the mineral could precipitate. Conversely, a negative SI indicates that the solution is undersaturated with respect to the selected mineral, suggesting that the mineral is dissolved in groundwater to reach equilibrium.

#### *2.8. Calculation of NBLs of Cr*

The assessment of NBLs for the target parameter was implemented based on the BRIDGE methodology [11]. The applied modified multi-method was separated into three steps: (a) the hydrogeochemical (bivariate plots, Piper, SI), (b) the PS method, and (c) the statistical analysis for estimating the NBLs (box plots for outliers, Q–Q (quantile–quantile) plots, and normality tests).

The applied methodology for the assessment of NBLs of Cr is described in detail in Figure 4. The pre-selection (PS) method, which is widely applied worldwide, was employed to select the suitable spring water samples for the NBLs assessment [10,16,57–60]. The PS method constitutes the methodology geochemical approach to validate the dataset according to similar geochemical characteristics and recognize the water samples that are affected by anthropogenic activities. In the first stage, the hydrochemical facies were selected based on the DO concentrations and Eh (ORP) [61]. The first dataset group included the water samples with ORP > 100 mV and DO > 3 mg/L. All water samples from the natural springs satisfied this criterion. The next criterion included consideration of redox conditions; if the prevailing conditions were oxidizing, then the concentrations of NO3 − < 10 mg/L would be considered and if the conditions were reducing, then the NH4 + < 0.5 mg/L would be considered [9,61,62] to exclude the samples affected by anthropogenic activities [62,63]. Based on this criterion, the samples with NO3 − > 10 mg/L were considered to be affected by anthropogenic activities and thus, were excluded from the new dataset. The next criterion required eight measurements per year for two years or two measurements per year for at least four years to exist for each spring [16].

#### *2.9. Threshold Values (TVs) Derivation*

The assessment of TVs was based on three scenarios [11,58] (Figure 5). The reference value was set equal to the water drinking acceptable limit (i.e., World Health Organization (WHO) guideline value).

#### *2.10. Meteoric Genesis Index (MGI)*

The meteoric genesis index (MGI) was also employed to classify the groundwater sources based on the depth of the meteoric water [64]. This index was calculated using the following Equation (5):

$$\text{Fr}\_2 = \text{Na}^+ + \text{K}^+ - \frac{\text{Cl}^-}{\text{SO}\_4^{2-}}, \text{ all concentrations are expressed in meq/L} \tag{5}$$

when r2 < 1, the groundwater source is of deep meteoric water percolation type whereas when r2 > 1 the groundwater is of shallow meteoric water percolation type [65].

#### *2.11. Meteorological Data*

Daily and monthly rainfall data for the period 2014–2019 were evaluated from the meteorological station in the Ermakia village (40◦30′325′′ N, 21◦51′233′′ E) which is the most representative and the nearest one, located on Vermio Mt., at an elevation of 1100 m. The annual precipitation for the period 2014–2019 was estimated at 985 mm.

**Figure 4.** Flow chart of the modified conceptual model for the assessment NBLs of the target parameter.

**Figure 5.** Flow chart for the assessment TVs for the target parameter.

#### **3. Results**

#### *3.1. Chemometric Analysis*

ଶ = ା + ା − ష ௌைర మష ′ ″ ′ ″ All chemical analyses were grouped into the above-mentioned seven categories according to their location, lithology, and type to evaluate the results. The descriptive statistics (max, min, median) of the physical and chemical parameters, for the studied spring waters from March 2014 until October 2020, are summarized in Tables 1 and 2. Some concentrations above the detection limit (DL) and below the quantification limit (QL) were not excluded because their very low concentrations do not affect the data processing of this work. All concentrations determined by ICP–MS presented a QL equal to the DL. The parameters NO2 −, NH4 + , PO4 3−, Fe, Cd, Co were measured below the detection limit (BDL) in the majority of the samples. Silver, Au, Be, Bi, Cs, Pt, Re, Ga, Ge, Hf, Hg, In, Mo, Nb, Ta, Ti, Th, Tl, W, Zr were detected BDL in all water samples. Hence, these elements were excluded from the statistical analyses. In Table 3, the dataset for Cr-Cr(VI) and the geographical coordinates for the sampling points are given.

In the Agios Panteleimonas area, the natural spring "Elafakia"-W14 and in the Agio Pnevma area, the natural spring "Potistis"-W13 were examined separately from the other springs in these areas because of their elevated groundwater concentrations of Cr.

Box plots depicting the variation of eight groundwater quality parameters (pH, Eh, DO, EC, Ca2+, Mg2+, HCO3 − and Si) are given in Figures 6 and 7. The alkalinity of the water samples was calculated in a range of 1.95 × 10−3 meq/L (S18) up to 6.71 × 10−3 meq/L "Potistis"-W13, with an average value of 4.60 × 10−3 meq/L.

− − The abundance of major ions varied significantly among the natural springs (Table 4). Concentrations of Cr and Cr(VI) exhibited a wide range of values (Tables 1–3, Figures 8–11). Most physical and chemical parameters (e.g., EC, NO3 −, SO4 2−, As, B, Ba, Cd, Cu, Ni, Pb, Sb, Se, Zn, etc.) in the natural springs of Vermio Mt. exhibited lower concentrations than the desirable and permissible limits for drinking water. All samples were compared with the guideline of WHO for drinking water [66]. Exceedances were recorded in the natural springs "Potistis"-W13, "Elafakia"-W14, S1-Agios Dimitrios area, and S18- Agio Pnevma area. Specifically, in the spring S18-Agio Pnevma area, concentrations of K+ exceeded the guideline value of 12 mg/L for drinking water (WHO, [66]). Two exceedances of As, above the permissible limits of 10 µg/L, were recorded in the natural spring S1 of the Agios Dimitrios area. In the spring "Elafakia"-W14 seven water samples exceeded the WHO guideline value of 50 µg/L for Cr concentration for drinking water [66], while in the spring "Potistis"-W13 23 water samples exceeded this value. Finally, in total, three samples (one from the spring "Elafakia"-W14, and two from the spring "Potistis"-W13) exceeded the limit of 20 µg/L for the concentration of Ni for drinking water [66]. The systematic exceedances of Cr in the natural springs of Vermio Mt. were the principal reason for selecting Cr as the target parameter for NBLs assessment in this area.

**Table 1.**Maximum, minimum, and median values of physical and chemical parameters of natural springs in western Vermio Mt.

> BDL: Below the detection limit. DL: Detection limit. QL: Quantification limit.

**Table 2.** Maximum, minimum, and median values of physical and chemical parameters of the natural springs in western Vermio Mt.

BDL: Below the detection limit. DL: Detection limit. QL: Quantification limit.

**Table 3.**Geographical coordinates of the water sampling sites and the results of Cr-Cr(VI) for the natural springs in western Vermio Mt.

**Figure 6.** Box plots for the physical parameters of the natural springs of western Vermio Mt.

**Figure 7.** Box plots for the major ions and Si of the natural springs of western Vermio Mt.

− −

‒

‒ ‒ ‒

‒

μ

**Table 4.** The abundance of major ions among the natural springs in western Vermio Mt.

**Figure 8.** Concentrations of Cr and Cr(VI) in the natural springs of western Vermio Mt.

‒ **Figure 10.** Concentrations of Cr and Cr(VI) in the natural spring "Potistis"-W13 of western Vermio Mt.

**Figure 11.** Statistical analyses of concentrations of Cr and Cr(VI) of the natural springs in western Vermio Mt.

#### *3.2. Correlation Analysis of Water Samples*

− − − − − − − − − − − − − − − − − − − − − − − − Spearman's rank correlation coefficient was applied to selected parameters (Table 5); the rest of the measured parameters were excluded because most of their values were BDL. The most remarkable features of the Spearman's rank correlation coefficients are the following: Cr presented a statistically significant (*p* < 0.01) very strong positive correlation coefficient with Cr(VI) (rs = 0.98), strong positive correlation coefficients with Mg2+ (rs = 0.76), Si (rs = 0.75), EC (rs = 0.71) and Ni (0.61), and moderate positive correlation coefficients with HCO3 − (rs = 0.55) and alkalinity (rs = 0.55). Hexavalent chromium exhibited a statistically significant (*p* < 0.01) very strong positive correlation coefficient with Mg2+ (rs = 0.8), a strong positive correlation coefficient with Si (rs = 0.76) and Ni (rs = 0.67), and moderate positive correlation coefficients with HCO3 − (rs = 0.59) and alkalinity (rs = 0.59). Magnesium exhibited a statistically significant (*p* < 0.01) strong positive correlation coefficients with EC (rs = 0.68), HCO3 − (rs = 0.68), alkalinity (rs = 0.68) and Ni (rs = 0.60), while very strong correlation coefficient with Si (rs = 0.91). Bicarbonates presented a statistically significant (*p* < 0.01) strong positive correlation coefficient with EC (rs = 0.62), and moderate positive correlation coefficients with Ni (rs = 0.51), Si (rs = 0.48). Sulfates had a statistically significant (*p* < 0.01) very strong positive correlation coefficient with U (rs = 0.82), strong positive correlation coefficients with Sr (rs = 0.74), Br (rs = 0.70), Ba (rs = 0.63), and Ca2+ (rs = 0.67), and a moderate positive correlation coefficient with Na+ (rs = 0.52). Arsenic exhibited statistically significant (*p* < 0.01) strong positive correlation coefficients with U (rs = 0.72), Ba (rs = 0.67) and V (rs = 0.70), and a moderate positive correlation coefficient with K+ (rs = 0.49).

**Table 5.**Spearman's rank correlation matrix of selected major and trace elements of the natural springs of western Vermio Mt. **Table 5.** *Cont.*

\* Correlation is significant at the 0.05 level (2-tailed). \*\* Correlation is significant at the 0.01 level (2-tailed).

#### **4. Discussion**

#### *4.1. Hydrogeochemical Characterization of the Natural Springs of Western Vermio Mt. the Ultramafic Fingerprint*

The dominant hydrochemical types of the studied natural springs of western Vermio Mt. were Ca-Mg-HCO3 (40% of the water samples), Mg-Ca-HCO3 (33% of the water samples), and Ca-HCO3 (21%) (Figure 12). Other transitional water types in the study area comprised Ca-K-HCO3-Cl (3%) and Ca-HCO3-SO4 (3%). Based on the type and the geological environment of the springs, the Ca–HCO3 waters are considered to originate through the interaction of meteoric water with rocks containing Ca-bearing minerals, whereas water types enriched in Mg, were derived from the dissolution of ultramafic rocks [29]. The mixed Ca-Mg-HCO3 type indicated fresh recharge waters mainly related to carbonate rocks and less to ultramafic rocks. The Mg-Ca-HCO3 water type represents recharge waters related to Mg-rich rocks, suggesting the strong interaction with the ultramafic rocks of the area [67]. The springs "Potistis"-W13, S1, and W21 that belong to this type are associated with fissured aquifers in strongly serpentinised ultramafic rocks, or they are in hydraulic connection with them. ‒

**Figure 12.** Piper diagram of major ion chemistry for the natural spring samples.

The water–ultramafic rock interaction typically produces Mg-HCO3 water type [33,40,68] and slightly alkaline to strongly alkaline pH conditions [69] because of the absorption of dissolved CO2 from atmospheric water in the serpentine and olivine according to the Equations (6) and (7) [70]:

$$\mathrm{Mg\_3Si\_2O\_5(OH)\_4 + 6CO\_2 + 5H\_2O \to 3Mg^{2+} + 6HCO\_3^- + 2H\_4SiO\_4} \tag{6}$$

$$\rm Mg\_2SiO\_4 + 4H\_2O + 4CO\_2 \to 2Mg^{2+} + 4HCO\_3^- + H\_4SiO\_4 \tag{7}$$

The pH values that characterised the studied springs varied from 7.3 up to 8.5, indicating near-neutral to slightly alkaline conditions that are typical of groundwater interacting with ultramafic and carbonate rocks [40,71]. Redox potential conditions were oxidizing up to strong oxidizing, as indicated by Eh, ranging from 300 mV up to 410 mV. The pH and Eh conditions in the studied springs favoured the release and solubility of the Cr oxyanion in groundwater since the solubility of oxyanions such as HCrO4 −, CrO4 2−, Cr2O7 2−, H2AsO4 −, and HAsO4 2− is enhanced with increasing pH [72].

Various tools are usually employed to define and evaluate the water–rock interaction processes in an area [73]. Bivariate plots of major ions and ionic ratios were used to study the hydrogeochemical evolution processes in the studied springs (Figure 13). In the bivariate plot of Ca vs. Mg, the water samples were grouped into three classes based on their Ca/Mg ratio (Figure 13a). In the first class belong the water samples with a Ca/Mg ratio below the 1:3 line. This class included two seasonal water samples from the spring "Potistis"-W13 (W13\_10\_19 and W13\_4\_19) with a Mg-Ca-HCO3 water type, revealing that the flow path was mainly through serpentinites. The second class contained the water samples of the natural spring S1-Agios Dimitrios area and "Potistis"-W13. They are all of Mg-Ca-HCO3 type and characterized by Ca/Mg ratios plotted below the 1:1 and above the 1:3 lines; this suggests a mixture of Mg-HCO3 and Ca-HCO3, indicating that these waters were derived from interaction with serpentinites and Ca-rich rocks. The third class comprised the rest of natural springs, characterized by mixed water types and a Ca/Mg ratio above the 1:1 line, indicating a limited influence of serpentinites. The bivariate plot of Ca + Mg vs. HCO3 (Figure 13b) suggests an excess of (Ca + Mg) over HCO3 reflecting an additional non-carbonate source of Ca2+ and Mg2+ ions, such as the dissolution of silicate minerals [38,57]. Iron-Mg-silicates of ultramafic rocks, such as olivine, pyroxene, and amphibole are transformed to serpentine group minerals during the serpentinisation process. Dissolution reactions favour the Mg2+ and HCO3 − release of the Mg-rich minerals (Equations (6) and (7)) [68].

During water–rock interaction, various chemical processes (e.g., fluctuation of ionic concentrations, mobilization of the dissolved components, and change in pH) are fingerprinted on the groundwater quality [74,75]. Gibbs diagrams are generally used to identify the hydrogeochemical evolution of groundwater, which involves precipitation, water–rock interaction, and evaporation–crystallization processes, based on TDS vs. Na+/(Na+ + Ca2+), and TDS vs. Cl−/(Cl−+ HCO3 −) scatter diagrams [76]. Herein, Gibbs diagrams were employed to assess hydrogeochemical processes that affect the water chemistry in the natural springs of western Vermio Mt. Figure 14 illustrates that all samples from natural springs fall into the water–rock interaction field, suggesting weathering of carbonate and silicate minerals. Although the use of Gibbs plots for groundwater has been disputed [77], the case study discussed herein exhibits none of the characteristics that could result in misuse of these plots (e.g., high SO4 2− concentrations, salinity sources, evolutionary flow paths, etc.). The implications of the Gibbs diagrams are in accordance with the calculated MGI index, according to which the waters from the natural springs are characterized as deep percolation types.

#### *4.2. Hydrogeochemistry of Cr in Natural Ultramafic Springs*

To further study the hydrogeochemistry of Cr in the studied springs, the average concentration of Cr was plotted vs. the water type of each spring (Figure 15). As shown, each water type is characterised by a wide range of concentrations of Cr, attributed to the different operation mechanisms of the spring and the weathering degree of the host geological formations. The mixed Mg-Ca-HCO3 water type ranges from very high concentrations of Cr, in the spring "Potistis"-W13 (>100 µg/L), to much lower values (<20 µg/L) in the springs S1-Agios Dimitrios area and W21-Agios Panteleimonas area. The mixed Ca-Mg-HCO3 water type is related to a range of concentrations of Cr from 17 to 48 µg/L. On the other hand, all springs that are characterised by a Ca-HCO3 water type exhibit very low Cr concentrations (<5 µg/L) since mostly the carbonate rocks influence their hydrochemistry. In all springs, the dominant anion is HCO3 −, the principal source of which is the dissolution of carbonate and silicate minerals [33].

− − −

−

‒

‒ ‒

− −

−

**Figure 13.** Bivariate plots of: (**a**) Ca vs. Mg and (**b**) (Ca+Mg) vs. HCO3 for the natural springs of western Vermio Mt.

An interesting feature of the spring "Potistis"-W13, derived from the evaluation of hydrogeochemical, hydrological, and meteorological data, is the decrease in concentrations of Cr in a very short time after rainfall; this is further supported by the strong linear regression of Cr vs. discharge (coefficient of determination R2 = 0.85) (Figures 16 and 17). Low discharge results in increased water–ultramafic rock contact time and thus, in elevated concentrations of Cr.

− − −

−

**Figure 14.** Gibbs diagrams of the natural springs in western Vermio Mt.

Μ

**Figure 15.** Average concentrations of Cr vs. water types of the natural springs in western Vermio Mt.

**Figure 16.** Fluctuation diagrams of rainfall and Cr concentrations of the natural spring "Potistis"-W13 in western Vermio Mt.

μ

μ ‒ **Figure 17.** Plot of Cr (µg/L) vs. discharge (L/h) in the natural spring "Potistis"-W13 in western Vermio Mt.

μ μ ‒ μ ‒ In Figure 18a, the Mg/Si ratio vs. Cr in the springs is presented. The diagram is divided into three sub-groups according to the concentrations of Cr. The water samples with low concentrations of Cr (<30 µg/L) constituted 45.5% of the total samples, 73.3% of which exhibited a Mg/Si ratio lower than 2. Concentrations of Cr from 30 µg/L up to 50 µg/L, comprised 25% of the total water samples, 77.78% of which exhibited a Mg/Si ratio lower than 2.3; only the samples W13\_6\_18 and W13\_7\_18 which correspond to the lowest concentrations of Cr recorded in the spring "Potistis"-W13 exhibited a Mg/Si ratio higher than 2.3. Of the total water samples, 42% exceeded the permissible limit of 50 µg/L for drinking water [66], with most of them corresponding to samples from the spring "Potistis"-W13. Most samples presented a Mg/Si ratio higher than 2.3. Respective Mg/Si ratios have been reported for groundwater in other natural ultramafic environments [29,33].

‒

The strong fingerprint of the water–rock interaction on the spring water chemistry and the geogenic origin of Cr in groundwater are indicated by the statistically significant very strong positive correlation coefficient of Cr with Si, the strong positive correlation coefficients of Cr with Mg2+, EC, and Ni, and the moderate positive correlation coefficients of Cr with HCO3 and alkalinity. Magnesium and alkalinity are two parameters usually increased with increasing degree of weathering; the latter has been reported to relate to elevated concentrations of Cr in groundwater [78]. Nickel is derived from the dissolution of Ni-bearing silicates which are released to groundwater under morphological and geochemical conditions that do not favour the occurrence of Fe-hydroxides and other secondary minerals capable of adsorbing Ni [79]. Nickel exhibited statistically significant, moderate positive correlation coefficients with Mg2+, EC, and Si, further highlighting its geogenic origin. The two natural springs with high concentrations of Ni ("Potistis"-W13, "Elafakia"- W14), also exhibit high mean concentrations of dissolved Si, and are of Mg-Ca-HCO3 and Ca-Mg-HCO3 water type. A similar case of high concentrations of Cr and Ni in Mg-HCO3 groundwater has been reported by Margiotta et al. [40]. Unlike the spring waters, Cr in the irrigational wells in the lowland of the Sarigkiol Basin was reported to strongly correlate with NO3 − and P, indicating the synergistic role of the agricultural activities [45].

**Figure 18.** Plots of: (**a**) Mg/Si vs. Cr and (**b**) Cr(VI) vs. Cr (regression model) of the natural springs in western Vermio Mt.

The statistically significant, very strong positive correlation coefficients of Cr with Cr(VI) (Spearman's rank correlation coefficient rs = 0.986) was further proven by their linear regression with a strong linear relationship (coefficient of determination R2 = 0.99, Figure 18b). In the analysed water samples, the Cr(VI)/Cr ratio ranges from 20% up to 100%. Specifically: (a) 62.5–90.3% in the spring "Mouratidis"-S2, (b) 80–99% in the spring S1- Agios Dimitrios area, (c) 20–62.5% in the springs at the Vazelonas area, (d) 39–100% in the springs at the Agio Pnevma area, (e) 40–99% in the spring "Elafakia"-W14, (f) 76–100% in the spring "Potistis"-W13). The fluctuation in the Cr(VI)/Cr ratio depends on the prevailing geochemical conditions (redox reactions, pH), the presence of iron or manganese oxides, and competing anions in each area, suggesting that various processes take place [78]. Hexavalent chromium is the principal form of Cr in the natural water springs ("Potistis"- W13, S1-Agios Dimitrios area, "Mouratidis"-S2); several factors contribute to the high Cr(VI)/Cr ratio. Specifically, the geological environment, which is enriched in Ca and Mg-bearing minerals, enhances Cr(VI) to form complexes with Mg and Ca and inhibits Cr(VI) reduction [71]. The presence of manganese oxides enhances the Cr(III) oxidation to

Cr(VI) in ultramafic rocks, soils, and unsaturated zone releasing Cr(VI) to groundwater (Equation (8)) [24].

$$\text{Cr(III)} + 1.5\text{\textdegree } - \text{MnO}\_{2(s)} + \text{H}\_2\text{O} \rightarrow \text{HCrO}\_4^- + 1.5\text{Mn}^{2+} + \text{H}^+ \tag{8}$$

The pronounced role of minerals in the concentrations of Cr in natural springs was investigated via SIs of selected mineral phases present in the study area (Figure 19). The percentage distribution of the SIs for the selected mineral phases is given in Table 6 for all collected water samples.

**Figure 19.** Plots of saturation indices vs. Cr of the natural springs in western Vermio Mt.

μ ‒ ‒ The water samples from the natural springs in which concentrations of Cr exceeded 50 µg/L ("Potistis"-W13 and "Elafakia"-W14) were oversaturated in: (a) the carbonate mineral calcite (100%) and (b) the oxide minerals chromite (100%) and magnetite (100%). On the other hand, they are undersaturated in: (a) the serpentine group minerals lizardite (100%) and chrysotile (100%), (b) the pyroxenes enstatite (100%) and diopside (100%), and the amphibole tremolite (100%) and (c) the olivine (100%).

− In general, the mineralogical phases that appear undersaturated tend to dissolve in water. The dissolution reactions contribute major, minor elements and PTEs to the groundwater. Chromium-bearing silicate minerals (serpentine, amphibole, pyroxene, chlorite, talc) occurred mostly undersaturated in the water samples, whereas Cr-rich oxides (chromite and Cr-magnetite) were oversaturated in the water samples. Therefore, silicate minerals are the principal geogenic contributors of Cr and other major/minor elements (e.g., Mg2+, Ca2+, HCO3 −, Si) and PTEs (e.g., As, Ni) to the spring waters of western Vermio Mt.

**Table 6.** The percentage distribution of the saturation state with respect to the selected mineral phases in the collected water samples.

#### *4.3. NBLs of Cr in the Ultramafic Environment of Vermio Mt.*

The geochemical characteristics of the natural springs, the geological environment, and the water–ultramafic rock interaction are reflected in springs' water quality. Chromium constitutes the principal environmental component in groundwater of the Sarigkiol Basin, originating primarily from geogenic and incidentally from anthropogenic sources [45]. This paper aims to assess the NBLs of Cr, which is of great interest in the catchment scale of the Sarigkiol Basin.

Based on the above-discussed hydrogeochemical data (e.g., pH, DO, Eh, Mg2+, Si, Cr, alkalinity, etc.), the most representative natural springs, which flow through and interact with ultramafic rocks, are the S1-Agios Dimitrios area, "Mouratidis"-S2, "Potistis"-W13, and "Elafakia"-W14.