Hydrogeochemical and Hydrodynamic Assessment of Tirnavos Basin, Central Greece

: A combined hydrogeochemical and hydrodynamic characterization for the assessment of key aspects related to groundwater resources management was performed in a highly productive agricultural basin of the Thessaly region in central Greece. A complementary suite of tools and methods—including graphical processing, hydrogeochemical modeling, multivariate statistics and environmental isotopes—have been applied to a comprehensive dataset of physicochemical analyses and water level measurements. Results revealed that the initial hydrogeochemistry of groundwater was progressively impacted by secondary phenomena (e.g., ion exchange and redox reactions) which were clearly delineated into distinct zones according to data processing. The progressive evolution of groundwater was further veriﬁed by the variation of the saturation indices of critical minerals. In addition, the combined use of water level measurements delineated the major pathways of groundwater ﬂow. Interestingly, the additional joint assessment of environmental isotopes revealed a new pathway from E–NE (which had never before been validated), thus highlighting the importance of the joint tools/methods application in complex scientiﬁc tasks. The application of multivariate statistics identiﬁed the dominant processes that control hydrogeochemistry and ﬁt well with identiﬁed hydrodynamic mechanisms. These included (as dominant factor) the salinization impact due to the combined use of irrigation water return and evaporitic mineral leaching, as well as the impact of the geogenic calcareous substrate (mainly karstic calcareous formations and dolostones). Secondary factors, acting as processes (e.g., redox and ion exchange), were identiﬁed and found to be in line with initial assessment, thus validating the overall characterization. Finally, the outcomes may prove to be valuable in the progression toward sustainable groundwater resources management. The results have provided spatial and temporal information for signiﬁcant parameters, sources, and processes—which, as a methodological approach, could be adopted in similar cases of other catchments.


Introduction
Groundwater is a critical natural resource that needs to be properly managed in order to sustain its paramount aspects of quantity and quality. In service of this goal, rational groundwater management requires accurate information and adequate knowledge about the processes affecting groundwater evolution in time and space. Therefore, scientists (and eventually stakeholders) need robust tools to efficiently evaluate the status of groundwater and decipher the governing factors that regulate its hydrogeochemical and hydrodynamic regimes. The acquisition and interpretation of such information is a demanding and complex task, which requires multidisciplinary approaches, different perspectives and varied methodologies [1]. Essentially, hydrogeochemical and hydrodynamic information needs to be dealt with holistically, with special emphasis put on their interactions. The methodological approach of the present work dictated the synergetic and combinational -Combine and test variable methodologies and tools through a specific proposed workflow, which may act a basic methodological array for assessing the hydrogeochemical and hydrodynamic conditions in similar cases. -Develop and optimize a conceptual model for the groundwater resources of the study area, based on previous literature and the newly applied combined methods/tools.
Performed analyses served as a basis for developing a deep understanding of the system's characteristics, in order to identify and adopt optimal solutions for rational groundwater management. The latter is important, especially in environmentally vulnerable arid or semiarid Mediterranean areas-like the Tirnavos basin-which are expected to face adverse impacts (e.g., increased temperatures, decreased precipitation, and reduced recharge of groundwater resources) in the future due to climate change [24]. Therefore, comprehensive planning in the framework of climate change is regarded as highly important for safeguarding food safety and overall sustainability of the water-ecosystem-food nexus on which the socioeconomic stability and welfare of the region depend.
It is apparent that, in complex aquifer systems that are strongly controlled by active tectonics, co-evaluation of numerous analysis methodologies of hydrological, hydrolithological, hydrogeochemical and isotopic data need to be conjunctively considered. That is the optimal way to get safe results and decipher controlling mechanisms that may prove to be critical in intensively stressed systems. The performed analysis confirmed known key hydrodynamic evolution mechanisms; it also revealed significant elements that were not apparent-having been masked by the controlling hydraulic and hydrogeochemical mechanisms.

Study Area
Thessaly plain is the largest alluvial basin in Greece, with a total area of 13,142 km 2 and an average altitude of 427 m. It is in central Greece and-through the Mid-Thessalic hills-it is divided into two sub-basins: the western Thessaly basin and the eastern Thessaly basin, each of which developed in a NW-SE direction as part of a wider tectonic trough. The Tirnavos subbasin, from a hydrological perspective, is essentially the northwest part of the eastern Thessaly basin. It includes sections of the Titarisios River basin to the north and the Pinios River basin to the south. Its area is estimated to 251 km 2 , or about 2.35% of the total area of the Pinios basin ( Figure 1). The perimeter of the study area was 106 km and the average altitude 76 m. The smallest and largest morphological gradients recorded are 0 and 45.93%, respectively, while the mean slope is 5.3%. The relief of the wider region is depicted in Figure 1. The climate is typical Mediterranean, with annual rainfall from 400 mm to 500 mm, distributed almost entirely during the wet hydrological period, without any significant precipitation during summer. Thus, several irrigation systems have been developed to support the cultivation of highly productive summer crops [25]. With regards to the meteorological data, Larissa station is the closest to the study area, and has gathered continuous and reliable data over several decades.
Agriculture is the dominant land use, covering approximately 165 km 2 or 65.73% of the Tirnavos sub-basin. Intensified agricultural activities, including both cultivation and livestock, are a major source of groundwater contamination by nitrogen compounds. Manure waste and the often excessive, improper use of nitrogen fertilizers aiming to improve agricultural production have led to the occurrence of elevated concentrations of nitrates in groundwater [26].

Geology-Hydrogeology
The Tirnavos subbasin forms the northwestern part of the eastern Thessaly plain of Central Greece. It is filled by Quaternary alluvial formations which are bounded along the southwest part of the basin by Neogene marls and sandy-clay deposits. The western margins consist of karstified marbles of middle-upper Cretaceous origin. The crystalline bedrock is composed of mica-schists and gneisses of upper Paleozoic and Paleozoic age, respectively, forming the northern boundary of the subbasin (Figure 2). Two major springs (Mati Tirnavou and Agia Anna) emerge at the contact of the karstified system with the alluvial deposits. The Pinios and Titarisios rivers flow across the subbasin which, as already mentioned, hydrologically part of the wider Pinios River basin. The climate is typical Mediterranean, with annual rainfall from 400 mm to 500 mm, distributed almost entirely during the wet hydrological period, without any significant precipitation during summer. Thus, several irrigation systems have been developed to support the cultivation of highly productive summer crops [25]. With regards to the meteorological data, Larissa station is the closest to the study area, and has gathered continuous and reliable data over several decades.
Agriculture is the dominant land use, covering approximately 165 km 2 or 65.73% of the Tirnavos sub-basin. Intensified agricultural activities, including both cultivation and livestock, are a major source of groundwater contamination by nitrogen compounds. Manure waste and the often excessive, improper use of nitrogen fertilizers aiming to improve agricultural production have led to the occurrence of elevated concentrations of nitrates in groundwater [26].

Geology-Hydrogeology
The Tirnavos subbasin forms the northwestern part of the eastern Thessaly plain of Central Greece. It is filled by Quaternary alluvial formations which are bounded along the southwest part of the basin by Neogene marls and sandy-clay deposits. The western margins consist of karstified marbles of middle-upper Cretaceous origin. The crystalline bedrock is composed of mica-schists and gneisses of upper Paleozoic and Paleozoic age, respectively, forming the northern boundary of the subbasin (Figure 2). Two major springs (Mati Tirnavou and Agia Anna) emerge at the contact of the karstified system with the alluvial deposits. The Pinios and Titarisios rivers flow across the subbasin which, as already mentioned, hydrologically part of the wider Pinios River basin.  [27,28] and [29], respectively.
The Quaternary deposits host an unconfined aquifer near the talus cones of Titarisios River in the NW, which converts to a phreatic and deeper confined aquifer toward the center of the basin. These are separated by a sequence of clay layers, which form an aquitard [29,30] and have a maximum thickness of over 550 m at their central parts [25,31]. Despite the limited potential of the phreatic aquifer nowadays, its importance is paramount as it forms an effective buffer zone that protects the confined units from potential surface-released contaminants of anthropogenic origin (e.g., due to agricultural activities). In addition, the marbles of the western margins host a karstic aquifer of great potential, which recharges the alluvial system by lateral crossflows.
Due to the occurrence of marls along the southern part of the karst, the inflow from this part is reduced compared to that of the northern parts. Crossflow from the crystalline bedrock at the northern margins of the basin also occurs but is of minor importance. In addition to the above, the southern extent of the karst system-and to a minor extent the Mid-Thessalic hills-recharges the central plain parts by crossflow from the southwest and southern parts of the area.
In the western part of the basin, there is an extensive marginal cone through which the alluvial system receives significant amounts of recharge, as crossflow through the Titarisios River gorge sediments. A smaller volume recharges the aquifer system as crossflow from the Pinios River gorge sediments to the south [29].

Methodological Array
A suite of different methodological tools was efficiently combined to provide a robust array of methodological workflows (Figure 3), having as their goal the construction of a comprehensive conceptual model for the evolution of groundwater resources at the study  [27,28] and [29], respectively.
The Quaternary deposits host an unconfined aquifer near the talus cones of Titarisios River in the NW, which converts to a phreatic and deeper confined aquifer toward the center of the basin. These are separated by a sequence of clay layers, which form an aquitard [29,30] and have a maximum thickness of over 550 m at their central parts [25,31]. Despite the limited potential of the phreatic aquifer nowadays, its importance is paramount as it forms an effective buffer zone that protects the confined units from potential surfacereleased contaminants of anthropogenic origin (e.g., due to agricultural activities). In addition, the marbles of the western margins host a karstic aquifer of great potential, which recharges the alluvial system by lateral crossflows.
Due to the occurrence of marls along the southern part of the karst, the inflow from this part is reduced compared to that of the northern parts. Crossflow from the crystalline bedrock at the northern margins of the basin also occurs but is of minor importance. In addition to the above, the southern extent of the karst system-and to a minor extent the Mid-Thessalic hills-recharges the central plain parts by crossflow from the southwest and southern parts of the area.
In the western part of the basin, there is an extensive marginal cone through which the alluvial system receives significant amounts of recharge, as crossflow through the Titarisios River gorge sediments. A smaller volume recharges the aquifer system as crossflow from the Pinios River gorge sediments to the south [29].
area. The bases of this workflow were the physicochemical analyses of groundwater samples, which were subsequently processed with the aid of (i) graphical methods (expanded Durov diagram), (ii) multivariate statistics (R-mode factor analysis and cluster analysis) and (iii) hydrochemical sections and bivariate plots. This combined approach was essential for grouping samples of similar hydrogeochemical identity (physicochemical fingerprint and processes). It identified the dominant previous and/or ongoing geochemical processes which strongly influenced the hydrogeochemical status of the study area. It also revealed the main sources of enrichment and/or contamination to/from specific parameters and allowed for the an assessment of their spatial evolution-progressively with the groundwater flow of the system. As a first result of this assessment, an evaluation of the groundwater quality was achieved, based on a comparison of the parametric values with relative standards and legislation. An additional set of raw data and field measurements-including groundwater level registrations-was then used for deciphering the basic hydrogeological characteristics (e.g., groundwater flow). This, combined with the results from stable isotope analysis (oxygen and deuterium) can provide insight into the recharge and hydrodynamic conditions of the system. Eventually, the combination and co-assessment of the hydrodynamic conditions with the dominant geochemical processes will constitute the fundamental aspect of the conceptual model for the groundwater resources of the area. Based on distinct steps and specified tools, this methodological array may be used as a model in similar cases where the groundwater conceptual model is the main scientific quest.

Groundwater Sampling, Analyses and Measurements
In total, 174 water samples were collected during 4 sampling periods (2 dry and 2 wet periods) between September 2016 and April 2018. Specifically, 153 samples were collected from boreholes within the study area and 21 surface water samples were taken from As a first result of this assessment, an evaluation of the groundwater quality was achieved, based on a comparison of the parametric values with relative standards and legislation. An additional set of raw data and field measurements-including groundwater level registrations-was then used for deciphering the basic hydrogeological characteristics (e.g., groundwater flow). This, combined with the results from stable isotope analysis (oxygen and deuterium) can provide insight into the recharge and hydrodynamic conditions of the system. Eventually, the combination and co-assessment of the hydrodynamic conditions with the dominant geochemical processes will constitute the fundamental aspect of the conceptual model for the groundwater resources of the area. Based on distinct steps and specified tools, this methodological array may be used as a model in similar cases where the groundwater conceptual model is the main scientific quest.

Groundwater Sampling, Analyses and Measurements
In total, 174 water samples were collected during 4 sampling periods (2 dry and 2 wet periods) between September 2016 and April 2018. Specifically, 153 samples were collected from boreholes within the study area and 21 surface water samples were taken from the rivers. Water samples were collected in double-capped polyethylene bottles of 1000 mL and stored in cool conditions during transport to a laboratory at The Soil and Water Resources Institute (SWRI) for analyses. Electrical conductivity (EC), pH, dissolved oxygen (DO) and oxidation-reduction potential (ORP) were determined in the field using the ProDSS (YSI inc., Yellow Springs, OH, USA) multiparametric probe. Laboratory analyses determined 27 parameters, including major and minor ions and trace elements; additional parameters were also calculated subsequently (Appendix A, Table A1). Calcium, magnesium and trace elements (except boron) were determined using an atomic absorption spectrophotometer by flame method. Sodium and potassium were analyzed through flame photometry while sulphate, nitrate, nitrite and ammonium were analyzed using spectrophotometer UV-VIS measurement. Bicarbonates, carbonate anions (neutralization by H 2 SO 4 ) and chloride (neutralization by AgNO 3 ) were determined volumetrically. The reliability of the results was determined by ionic balance error, which was found to be less than 10% for all samples, with a median value of −3%.
In addition, twenty-six (26) samples were collected for isotopic analyses from selected sites of the monitoring network in April (n = 15) and September (n = 11) 2018, respectively, according to the International Atomic Energy Agency's sampling specifications [32]. The first isotope sampling period took place in April 2018; 15 samples were collected. The second sampling period took place in September 2018; only 11 out of 15 samples were collected because some of the wells were not operating. Additionally, the Titarisios River was dry during that period. In detail, the following samples were collected: (a) 21 from groundwater sites (12 in the first period, 9 in the second); (b) 2 from the Mati Tirnavou spring (1 for each period); (c) 2 from the Pinios River (1 for each period); and (d) 1 from the Titarisios River (in the first period). Fully filled 100 mL double cap polyethylene bottles were used for sampling, and samples were analyzed for oxygen and hydrogen isotopic compositions ( 18 O and 2 H) at the Hydrology Laboratory of Lubeck Technical University in Germany. Analyses were conducted via off-axis integrated cavity output spectroscopy (OA-ICOS, DLT-100 Liquid Water Isotope Analyzer, Los Gatos Research Inc., Mountain View, CA, USA) and reported in per mil (‰). The analytical precision for δ 18 O and δ 2 H was 0.2‰ and 0.6‰, respectively.
In addition to water sampling, a groundwater level monitoring network comprising 46 wells was compiled ( Figure 4) and operated to facilitate the objectives of the research. At all groundwater level measuring stations (boreholes), the altitude was accurately recorded using the Real Time Kinematic (RTK)-GPS, Epoch 50 by Spectra Geospatial. Water level monitoring was performed in April and September 2017, to reflect representative conditions of the wet and dry hydrological periods. Based on the compiled piezometric maps (Figure 4), a relative drawdown is evident, ranging from a few to ten meters locally. This indicated a seasonal effect on water level and consequently volume of water stored in the aquifer system.

Data Processing
Groundwater chemistry data was processed with graphical approaches and multivariate statistics. An expanded Durov diagram [34] was used to identify the ongoing hydrogeochemical processes in the study area. In this diagram, the cation and anion triangles were recognized and distinguished along the 25% axes so that the main field was conveniently divided. The expanded Durov diagram has a distinct advantage over the Piper diagram [35] in that it provides a better display of hydrochemical water types [36] and has the potential to reveal geochemical processes that could affect groundwater evolution [37].

Data Processing
Groundwater chemistry data was processed with graphical approaches and multivariate statistics. An expanded Durov diagram [34] was used to identify the ongoing hydrogeochemical processes in the study area. In this diagram, the cation and anion triangles were recognized and distinguished along the 25% axes so that the main field was conveniently divided. The expanded Durov diagram has a distinct advantage over the Piper diagram [35] in that it provides a better display of hydrochemical water types [36] and has the potential to reveal geochemical processes that could affect groundwater evolution [37].
To get a better insight on the potential correlation between the examined parameters, the Pearson correlation coefficient (r) was calculated. If the correlation coefficient (r) is greater than 0.7, parameters are considered strongly correlated; if it ranges between 0.5 and 0.7, it indicates a moderate correlation at a significance level of p < 0.05. [38]. In addition, a classic robust multivariate statistical technique (R-mode factor analysis) was applied. Factor analysis (FA) has been widely used in environmental sciences and hydrogeochemical research [1,14,39,40]. It is a multivariate statistical technique that involves linear combinations of variables through a correlation-focused approach. FA seeks to reproduce To get a better insight on the potential correlation between the examined parameters, the Pearson correlation coefficient (r) was calculated. If the correlation coefficient (r) is greater than 0.7, parameters are considered strongly correlated; if it ranges between 0.5 and 0.7, it indicates a moderate correlation at a significance level of p < 0.05. [38]. In addition, a classic robust multivariate statistical technique (R-mode factor analysis) was applied. Factor analysis (FA) has been widely used in environmental sciences and hydrogeochemical research [1,14,39,40]. It is a multivariate statistical technique that involves linear combinations of variables through a correlation-focused approach. FA seeks to reproduce intercorrelations among variables, in which the factors represent their common variance. The relationship among several observed quantitative variables is represented in terms of a few underlying, independent variables, called factors, which may not be directly measured or even measurable [41]. The exact number of factors was chosen by Kaiser criterion [42] in which factors with eigenvalues smaller than 1 were eliminated. Prior to processing, initial data was standardized through log transformation and z-scores to eliminate the influence of different units between variables and acquire a normal or log normal distribution [43]. Variables that failed to meet that criterion were omitted from further data processing. Finally, 19 parameters compiled the correlation matrix that accounted for the degree of mutually shared variability between individual pairs of groundwater quality parameters. The higher the factor loading of a parameter, the greater its participation to the examined factor. Factor loadings above 0.750 were considered high, those between 0.500 and 0.750 were considered moderate, and those between 0.400 and 0.500 were considered weak [40]. To further assist the performed assessments, the saturation indices of critical minerals were calculated with the aid of PHREEQC software [44].

Physicochemcial Analyses
Electric conductivity (EC) varied from 253 to 1821 µS/cm (Figure 5a). The higher EC values occurred in the southeastern part of the study area, where most of the region's smallto-medium sized industrial units are located [45]. The lowest values of EC are noted in the eastern and central part of the study area, where no significant industrial activity occurs.
Groundwater temperature is known to be an important driver for water quality [46,47], and therefore a crucial parameter for groundwater quality management. In addition, groundwater temperature is one of the best environmental tracers for detecting water flux, because heat in aquifers is transported both by conduction and advection caused by groundwater flow water temperature [48]; hence it may be directly related to the hydrodynamic evolution of a system. Low groundwater temperatures normally indicate recharge waters while higher temperatures indicate greater distance from the recharge source. As shown in Figure 5b, the lower groundwater temperatures were observed in the western margins of the basin, at the contact of the alluvial with the carbonate formations, while more elevated temperatures occurred toward the central and southern parts of the basin. This spatial distribution of temperature-combined with the piezometric dataindicated that the karstic system to the west is the main recharge source. On the other hand, low temperatures were also observed in the eastern parts of the basin, where alluvial deposits exist between metamorphic formations of crystalline bedrock, probably denoting a potential recharge effect from the east. In the absence of monitoring points at that part of the basin, however, that hypothesis could not be verified by piezometric data.
The main recharge areas of the basin were also indicated by the spatial patterns of DO and pH. For DO, higher values suggested enrichment in oxygen; thus, recharge conditions were present in the western and northern parts, while DO values were progressively reduced toward the center of the basin (Figure 5c). Deviations from that pattern (e.g., NW part of the basin) may be attributed to potential secondary geochemical processes (e.g., redox reactions) which create locally reducing conditions. Regarding the pH values, it was evident that in the recharge areas (west and north), pH was circumneutral, but progressively more alkaline moving toward the central and eastern parts of the basin (Figure 5d).
Based on the analytical results (Appendix A, Table A1), groundwater samples were slightly alkaline (median value: med = 7.61) with relatively low electrical conductivity (med = 486 µS/cm). The order of abundance for cations is Ca 2+ > Mg 2+ > Na + > K + , and for anions HCO 3− > NO 3 − > SO 4 2− > Cl − . With regard to calcium and bicarbonates (which constitute the most abundant ions), it was evident that the main origin of their concentrations was the karstic system, cropping out at the western parts of the basin. That, as indicated in the previous section, is considered a major pathway for recharge of the system. Interestingly, nitrates are the second most abundant anion, clearly reflecting the anthropogenic impact due to widespread agricultural activities. In total, 72% of samples exhibited concentrations above 10 mg/L, which is an indicative threshold for nitrate contamination in natural systems [49]. Over 9% exceeded the Maximum Admissible Concentration (MAC) for drinking water (50 mg/L), according to [50,51], reaching up to 145 mg/L. However, their spatial distribution was scattered, without any profound pattern, elucidating that, apart from diffuse agricultural impact, contamination is point-source. This could be explained by agricultural malpractice and/or septic tanks, reflecting the potential influence of locally developed favorable conditions for percolation of pollutants through Water 2021, 13, 759 9 of 24 the vadose zone-which is in line with the above-described characteristics of the sediments that fill the basin. Other outlying values and cooccurrences of elevated ion concentrations (e.g., maximum values of Na + (287 mg/L) and SO 4 2− (604 mg/L)) were also related to local factors (e.g., soil amendments rich in Na + and SO 4 2− or dissolution of evaporitic minerals such as thenardite (Na 2 SO 4 ), gypsum (CaSO 4 ), and halite (NaCl)) which could have occurred due to previous paleoenvironmental evaporitic conditions [1]. Nevertheless, their distribution (and thus their impact) was very limited. Finally, concentrations of heavy metals and metalloids were generally low; the only exception was an outlier for B (1.25 mg/L), most likely attributed to local impact of boron-rich soil amendments.
influence of different units between variables and acquire a normal or log normal distri-bution [43]. Variables that failed to meet that criterion were omitted from further data processing. Finally, 19 parameters compiled the correlation matrix that accounted for the degree of mutually shared variability between individual pairs of groundwater quality parameters. The higher the factor loading of a parameter, the greater its participation to the examined factor. Factor loadings above 0.750 were considered high, those between 0.500 and 0.750 were considered moderate, and those between 0.400 and 0.500 were considered weak [40]. To further assist the performed assessments, the saturation indices of critical minerals were calculated with the aid of PHREEQC software [44].

Physicochemcial Analyses
Electric conductivity (EC) varied from 253 to 1821 μS/cm (Figure 5a). The higher EC values occurred in the southeastern part of the study area, where most of the region's small-to-medium sized industrial units are located [45]. The lowest values of EC are noted in the eastern and central part of the study area, where no significant industrial activity occurs. Groundwater temperature is known to be an important driver for water quality [46,47], and therefore a crucial parameter for groundwater quality management. In addition, groundwater temperature is one of the best environmental tracers for detecting water flux, because heat in aquifers is transported both by conduction and advection caused

Identification of Governing Processes and Spalital Evolution of Hydrogeochemistry
The results of chemical analyses were plotted in the expanded Durov diagram with the help of the DurovPwin application [52]. To ensure the representativeness of the results, the median values of the four periods were considered in the process. According to the expanded Durov diagram (Figure 6), two basic hydrochemical characters are identified, which can be further separated into subgroups: 1a, 1b, 2, 3, 5, 6. These characters are described below, while in Figure 6 the spatial distribution of the identified groups is illustrated.
(e.g., soil amendments rich in Na + and SO4 2− or dissolution of evaporitic minerals such as thenardite (Na2SO4), gypsum (CaSO4), and halite (NaCl)) which could have occurred due to previous paleoenvironmental evaporitic conditions [1]. Nevertheless, their distribution (and thus their impact) was very limited. Finally, concentrations of heavy metals and metalloids were generally low; the only exception was an outlier for B (1.25 mg/L), most likely attributed to local impact of boron-rich soil amendments.

Identification of Governing Processes and Spalital Evolution of Hydrogeochemistry
The results of chemical analyses were plotted in the expanded Durov diagram with the help of the DurovPwin application [52]. To ensure the representativeness of the results, the median values of the four periods were considered in the process. According to the expanded Durov diagram (Figure 6), two basic hydrochemical characters are identified, which can be further separated into subgroups: 1a, 1b, 2, 3, 5, 6. These characters are described below, while in Figure 6 the spatial distribution of the identified groups is illustrated. Figure 6. Expanded Durov plot illustrating the hydrogeochemical processes [36] and their spatial distribution in the study area. Figure 6. Expanded Durov plot illustrating the hydrogeochemical processes [36] and their spatial distribution in the study area.

•
Hydrochemical Type 1 (Ca-Mg-HCO 3 ): In Subgroup 1a, Ca 2+ and HCO 3 − were dominant, indicating that this was a recharge zone, which was verified by high dissolved oxygen and low water temperature values in this area (according to field measurements analyzed in previous paragraphs). Based on the dominant cations in the formation of the hydrochemical character, group 1b did not differ substantially from 1a. However, there was a slight increase in the concentrations of Na + and Mg 2+ and, with respect to the anions, higher SO 4 2− values. Spatially, this hydrochemical type occupied two zones on either side of the Titarisios River. The first began north of it and extended to the northwestern part of the basin, while the second began from the central part of the basin (south of Titarisios) and extended to the south and southwestern part of the basin. The dominance of Ca 2+ and HCO 3− in this subgroup also indicated a recharge zone, possibly lower than that of group 1a. The involvement of Mg 2+ in the formation of the hydrochemical type was consistent with the mild dolomitization process which characterized the carbonate system [27]. The reason for the observed differentiation between the two subgroups may have been due to the limited system recharge rate from the karst (which allows mixing of recharge water with the aquifer system and/or the development of a small degree of ion exchange processes). In particular, the development of this subgroup on its southwestern boundary may have been justified by the emergence of Neogene deposits in this area, which limited the rate of groundwater replenishment in this part of the aquifer. At the same time, in this area, according to earlier studies [53], lateral recharge from the Mid-Thessalic hills was reported-of limited extent and intensity and qualitatively inferior to that of the karst system. • Hydrochemical Type 2 (Ca-Mg-Na-HCO 3 -SO 4 ): The second hydrochemical type was characterized by higher Na + and SO 4 2− concentrations. The involvement of Ca 2+ and HCO 3 − in the formation of the hydrochemical type of subgroup 2 (Mg-Ca-HCO 3 -SO 4 ) was still significant; however, it was less than that of subgroups 1a and 1b. The projection position of subgroup 2 in the expanded Durov diagram indicated progressive involvement of ion exchange and perhaps also mixing mechanisms. At least locally and seasonally [27], the mixing of Pinios River water with aquifer water through filtration along its bed in the formation of the observed hydrochemical type cannot be excluded. It was also characterized by low analogues with subgroups 1a and 1b at Ca 2+ and HCO 3 − concentrations. Subgroup 2 (Ca-Mg-Na-HCO 3 ) can therefore be considered as representative of a transition zone from intense recharge zones (subgroups 1a and 1b) to a restricted recharge zone and longer groundwater residence time in the aquifer (subgroups 3, 5, 6). In subgroup 3 (Na-HCO 3 ), HCO 3 − and Na + dominated, and the ion exchange phenomenon was fully evolved, as indicated by the projection of its representative samples (wells 18 and 34) on the expanded Durov diagram. Subgroup 5 (Na-Mg-Ca-SO 4 -HCO 3 ) represented mixing or dissolution waters, where Mg 2+ (but mainly Na + from cations) and Cl -(but mainly SO 4 2− and HCO 3 − , to a lesser extent from anions) dominated, thus suggesting a contributing recharge mechanism in the form of lateral influx from surrounding formations that contributed to the aquifer balance and obviously affected its hydrochemical identity. Last, subgroup 6 (Na-SO 4 -HCO 3 ), which was only represented by a single but distinct sample, plotted on an uncommon part of the graph for groundwaters, which is often a product of mixing.
The above assessments may be further supported by the physicochemical evolution of groundwater along its flow path (Figures 7 and 8), from the identified recharge area to the end of the transition zone, as defined by the six (6) selected (representative) wells (A-A' axis) depicted in Figure 6. With regard to the physicochemical parameters (Figure 7), an increasing trend was identified for the EC (like Total Dissolved Solids -TDS-, as expected), following the typical ion enrichment of groundwater from the recharge areas to the central parts of the hydrological basin. Bicarbonates (HCO 3 − ) progressively decreased and then subsequently increased, reaching a higher value than the recharge area; this evolution probably denoted an external impact from CO 2 content. Specifically, the reaction of calcite from the karstic substrate in the recharge area (west) with the CO 2 derived from the respiration of organic matter, produced carbonic acid as an intermediate step and then bicarbonates-with simultaneous release of calcium, as shown in Equation (1).
Progressively, along the groundwater flow path, the CO 2 content increased (presumably using as its source the occurrence of organic matter in the sediments of the aquifer matrix due to the sedimentation process of the fluviolacustrine deposits under C rich paleogeographic conditions), leading to an increase of CaCO 3 solubility [54] and shifting the Saturation Index of calcite (SI CaCO3 ) values from supersaturated to slightly saturated. The latter was in accordance with Figure 8, indicating sharp SI CaCO3 reduction.
With respect to Cl − : a constant increase toward groundwater flow was identified. Bearing in mind that Cl − is a conservative ion, the increase of its concentration should be attributed to external factors which-based on the dominant land use activities and site characteristics-may be related to manure leaching and/or septic tanks. However, the dissolution of evaporitic minerals such as halite (NaCl) could not be excluded. Interestingly, the concentrations of SO 4 2− changed without a significant pattern, with notable variations.
The primary source of sulfates was geogenic, stemming from the weathering of S-bearing minerals (e.g., gypsum, pyrite, etc). Changes in their concentration were likely due to further impact from the S-bearing minerals, as well as from external factors, e.g., soil amendments. However, their concentration along the dominant groundwater flow path may also change due to variations in redox conditions. The latter was also evidenced by the similar trend in nitrate concentrations-which may act as electron donors. Nevertheless, in the case of nitrates, the changes in their concentrations were less sharp than sulfates. The redox impact was also evidenced in some wells through the elevated values of NH 4 + and the relatively low values of NO 3 − . Potential effects due to leaching of N-fertilizers could not be excluded. and then subsequently increased, reaching a higher value than the recharge area; this evolution probably denoted an external impact from CO2 content. Specifically, the reaction of calcite from the karstic substrate in the recharge area (west) with the CO2 derived from the respiration of organic matter, produced carbonic acid as an intermediate step and then bicarbonates-with simultaneous release of calcium, as shown in Equation (1).   and then subsequently increased, reaching a higher value than the recharge area; this evolution probably denoted an external impact from CO2 content. Specifically, the reaction of calcite from the karstic substrate in the recharge area (west) with the CO2 derived from the respiration of organic matter, produced carbonic acid as an intermediate step and then bicarbonates-with simultaneous release of calcium, as shown in Equation (1).   Regarding the major cations (Figure 8), there was an antithetic change in the concentrations of Ca 2+ and Na + , explained by the process of ion exchange. The fresh (recharge) groundwater (which was enriched in Ca 2+ due to the karstic substrate) is subject to ion exchange toward its flow with the Na + content of silicate minerals, leading progressively to a decrease in calcium and increase in sodium. This was also reflected in the shift of water types, from Ca-HCO 3 to Na-HCO 3 . The magnesium (Mg 2+ ) also showed a slight increase toward flow, probably enriched by the weathering of dolomite or other Mg-bearing aluminosilicate minerals. Potassium (K + ) was nearly constant, exhibiting only minor changes in concentration. Overall, the selected parameters' evolution along the regional groundwater flow direction (Figures 7 and 8) agreed with the findings and observations discussed above, supporting the concept of strong recharge from the karstic domain in the form of lateral crossflow which then mixes with resident groundwater and is impacted locally by leachates from anthropogenic activities. Moreover, it clearly denoted the progressive shift from mixing to ion exchange process, which was also affected by oxygen depletion and the relative dominance of reducing conditions. The dominance of the ion exchange process was further verified by the major ion relation analysis conducted and illustrated in Figure 9. Specifically, the relation between (Ca+Mg-SO 4 -HCO 3 ) and (Na+K-Cl) in meq/L was examined, as suggested by [55,56]. The product of (Na+K-Cl) represented the excessive sodium originating from sources other than halite dissolution (assuming that all chloride is derived from halite). In addition, the product of (Ca+Mg-SO 4 -HCO 3 ) represented the calcium and/or magnesium which originated from sources other than gypsum and carbonate dissolution (calcite and/or dolomite). If these processes were significant in defining the hydrogeochemistry of groundwater, the relation between the two products should have been linear with a slope of −1. As seen in Figure 9, many data points lie close to a straight line (r = +0.93) with a slope of −0.80, which clearly points to the existence of cation exchange [57,58].
According to [58] the plot of (Ca + Mg) vs. (SO 4 + HCO 3 ) could be a reliable indicator-helping to distinguish between ion exchange and reverse ion exchange processes, if these are active in a given study area. The points plotted below the 1:1 line toward the Mg + Ca axis suggest prevalence of the ion exchange process, while those plotted above the 1:1 line towards the SO 4 + HCO 3 axis denote prevalence of reverse ion exchange. The results of the plotted samples ( Figure 10) are in accordance with the outcomes of the expanded Durov diagram ( Figure 6). The relative positions of the points in the plot indicate that excessive calcium and magnesium in groundwater were exchanged with sodium from the aquifer matrix. Figure 10 shows the amount of Ca + Mg gained or lost relative to that provided by calcite, dolomite, and gypsum. When HCO 3 + SO 4 is low (< 5 meq/L) and the samples plot on 1:1 line, dissolution of calcite and dolomite is the major process influencing water chemistry; on the contrary, when HCO 3 + SO 4 is greater than 5 meq/L, in addition to calcite and dolomite, dissolution of gypsum is likely to occur [59].
According to the calculated Pearson coefficient (r), there were strong positive correlations between Na-Cl (+0.847), Na-SO 4 (+0.889), Cl-SO 4 (+0.855), B-Na (+0.899), and B-SO 4 (+0.904); there were moderate positive correlations between Ca-Mg (+0.532), Mg-Cl (+0.683), Mg-SO 4 (+0.612), B-Cl (+0.737) and Mg-NO 3 (+0.648). Most of the correlated parameters were indicative of a salinization impact, which, according to the dominant land use and site specs, can probably be attributed to leachates rich in high salinity fertilizers washed off through irrigation water return flow. It is interesting, though, that boron (B) was strongly correlated with Na, Cl and SO 4 , denoting a similar origin or enrichment process. Based on the hypothesis addressed above (occurrence of evaporitic minerals), B could potentially derive from the dissolution of buried evaporitic minerals such as borax (Na 2 B 4 O 7 ·10H 2 O), considering the paleoclimatic conditions of the area. Nevertheless, this cannot solely explain the correlation with Cl and SO 4 . Most likely, B is also related to fertilizers/soil amendments, whose leached products may contain boron, apart from other saline-related parameters. Water 2021, 13, x FOR PEER REVIEW 14 of 25 Figure 9. Plot of (Ca+Mg-SO4-HCO3) and (Na+K-Cl) in meq/L.
According to [58] the plot of (Ca + Mg) vs. (SO4 + HCO3) could be a reliable indicator-helping to distinguish between ion exchange and reverse ion exchange processes, if these are active in a given study area. The points plotted below the 1:1 line toward the Mg + Ca axis suggest prevalence of the ion exchange process, while those plotted above the 1:1 line towards the SO4 + HCO3 axis denote prevalence of reverse ion exchange. The results of the plotted samples ( Figure 10) are in accordance with the outcomes of the expanded Durov diagram ( Figure 6). The relative positions of the points in the plot indicate that excessive calcium and magnesium in groundwater were exchanged with sodium from the aquifer matrix. Figure 10 shows the amount of Ca + Mg gained or lost relative to that provided by calcite, dolomite, and gypsum. When HCO3 + SO4 is low (< 5 meq/L) and the samples plot on 1:1 line, dissolution of calcite and dolomite is the major process influencing water chemistry; on the contrary, when HCO3 + SO4 is greater than 5 meq/L, in addition to calcite and dolomite, dissolution of gypsum is likely to occur [59].

Multivariate Statistics
According to the calculated Pearson coefficient (r), there were strong positive correlations between Na-Cl (+0.847), Na-SO4 (+0.889), Cl-SO4 (+0.855), B-Na (+0.899), and B-SO4 (+0.904); there were moderate positive correlations between Ca-Mg (+0.532), Mg-Cl (+0.683), Mg-SO4 (+0.612), B-Cl (+0.737) and Mg-NO3 (+0.648). Most of the correlated parameters were indicative of a salinization impact, which, according to the dominant land use and site specs, can probably be attributed to leachates rich in high salinity fertilizers washed off through irrigation water return flow. It is interesting, though, that boron (B) was strongly correlated with Na, Cl and SO4, denoting a similar origin or enrichment process. Based on the hypothesis addressed above (occurrence of evaporitic minerals), B could potentially derive from the dissolution of buried evaporitic minerals such as borax (Na2B4O7·10H2O), considering the paleoclimatic conditions of the area. Nevertheless, this cannot solely explain the correlation with Cl and SO4. Most likely, B is also related to fer-

Multivariate Statistics
The results of the FA (Table 1) identified five (5) factors that define the hydrogeochemical regime and explained 83.3% of the total variance. The first factor (F1) accounted for the 32.8% of total variance and denoted the cumulative effect of groundwater salinization by variable sources. It included (with strong positive loadings) the parameters of EC, Na, Cl, B, and SI halite ; it also included (with medium loadings) the parameters of Mg and SI gypsum , reflecting the impact from irrigation water return flow and the dissolution of evaporitic minerals. It should be noted that higher values of SIs correlate with a greater state of saturation. Table 1. Results of Factor Analysis (FA) with rotated (varimax) factor loadings and communalities. Medium-to-high positive correlation is shown in bold, while medium-to-high negative correlation is shown in bold italics. The second factor (F2) accounted for the 19.6% of total variance and included: Ca and NO 3 with strong positive loadings; K and HCO 3 with medium ones; Mg and SI gypsum with weak ones; and pH with a strong negative loading. This factor may possibly interpret the hydrogeochemistry in the recharge areas, where Ca and HCO 3 contents are elevated (higher values), and pH (as shown previously) is circumneutral (lower values). This sufficiently explained the antithetic loading. The strong positive loading of nitrates possibly denoted an additional dominant process of nitrate enrichment in these areas. This was possibly due to existence of two additional criteria: the relative oxidizing conditions and the external source of NO 3 (e.g., agricultural activities and/or septic tanks); the latter being profoundly supported by the geometry of the aquifer system at that part of the basin, characterized by a phreatic unit of high hydraulic parameters.

Variable
The third factor (F3) explained 12.4% of total variance and included (with strong positive loadings) the SIs of calcite and dolomite. It probably denoted the areas which were supersaturated in the above minerals because of the calcareous substrate (limestones and dolostones). The weak correlation with the dominant cations of these formations (r = +0.35 for Ca-SI calcite , r = +0.45 for Mg-SI dolomite ) probably reflected the existence of additional sources for Ca and Mg (e.g., F1 and F2) which masked the direct covariance of these parameters due to impact from limestones and dolostones.
The fourth (F4) factor explained 11.9% of total variance, with strong negative loadings for Mn and NH 4 and a medium one for Fe. This factor clearly reflected the local reducing conditions, in which dominance of Fe 2+ and Mn 2+ prevailed along with ammonium (NH 4 ). These areas were probably affected by the increased organic content that creates reducing conditions and/or the anoxic environments caused by limited groundwater recharge (e.g., semiconfined and/or confined aquifers).
Finally, the fifth (F5) factor explained a minor percentage (6.6%) of total variance and included (with strong and medium positive loadings) the parameters of Cu and Fe, respectively. It probably denoted the occurrence of a weak local sulfide mineralization (e.g., chalcopyrite-CuFeS 2 or other) without excluding the potential impact from the use of Cu as a soil conditioner (mainly as a major constituent in several plant protection products, especially for orchards and vineyards).

Stable Isotopes
With regard to the stable isotopes, the Local Groundwater Isotope regression Line (LGIL) was compiled and plotted, along with the Global Meteoric Water Line (GMWL), Greek MWL and Thessalian MWL-as depicted in Figure 11. Isotope values of δ 2 H and δ 18 O were expressed as the difference between the measured ratios of the sample and reference divided by the measured ratio of the reference, which was in turn expressed as VSMOW values (Vienna-standard mean ocean water). The isotopic ratio of δ 18 O in the study area ranged between −9.8‰ and −6.9‰ with an average value of −8‰. This value almost coincided with the average value recorded from the spring waters of Thessaly (−8.29‰) [60]. Similarly, the isotopic ratio of δ 2 H ranged between −66.2‰ and −47.5‰ with an average value of −51.3‰ which was almost the same as the average value recorded in Thessaly (−51.1‰), based on the reference. Based on these observations, the isotope composition of groundwater in the Tirnavos subbasin seemed to fit perfectly with the obtained levels from a previous study on the Thessaly region. The function that expressed the LGIL (δ 2 H = 6.71 × δ 18 O + 2.09) compared to Thessaly's MWL (δ 2 H = 6.48 × δ 18 O + 1.7) presented a similar slope (6.71 and 6.48) and slightly higher value of the d-excess (2.09 and 1.7). Generally, deuterium excess is primarily controlled by kinetic effects associated with evaporation of water at the surface of the oceans or inland. It increases with an increase in the moisture deficit of oceanic air masses [61]. Differences in d-excess arise because of varying temperature, relative humidity, and wind speed at the sea surface, whereas global atmospheric moisture mainly originates from admixture of recycled continental vapor [19,60]. Hence, the value of d-excess (close to 2) in the study area, compared to the GMWL value (10) could probably be attributed to environmental conditions. Based on the results of the 18 O isotope analyses, a spatial distribution map of groundwater samples was constructed ( Figure 12), aiming to facilitate the identification and evolution of the active recharge mechanisms in the Tirnavos alluvial basin in support of hydrochemical and piezometric data analyses. For the construction of this map, the monitoring period of April 2018 was used because of the larger number of available samples.
Due to the influence of altitude and continentality on the isotopic composition of precipitation [64,65]-and based on the recharge mechanisms already identified and substantiated-it would be expected that the most negative values of δ 18 O were observed along the western part of the basin, where the alluvial aquifer system receives most of its natural recharge in the form of lateral crossflows from the karstic system of Tirnavos. This is the most distant boundary of the alluvial system from the sea, and receives water from direct infiltration of precipitation at altitudes that reach as high as 900 m. Likewise, moving to the downstream parts of the basin, towards the sea, the δ 18 O values would have been expected to become less negative. From the spatial distribution map (Figure 12), however, this only happened in the NW part of the basin while, as we move to its central and southeastern parts, the values were increasingly negative-with the minimum figure (the most negative value) found at monitoring point LB214, the easternmost of the compiled network.
(2.09 and 1.7). Generally, deuterium excess is primarily controlled by kinetic effects associated with evaporation of water at the surface of the oceans or inland. It increases with an increase in the moisture deficit of oceanic air masses [61]. Differences in d-excess arise because of varying temperature, relative humidity, and wind speed at the sea surface, whereas global atmospheric moisture mainly originates from admixture of recycled continental vapor [19,60]. Hence, the value of d-excess (close to 2) in the study area, compared to the GMWL value (10) could probably be attributed to environmental conditions. Based on the results of the 18 O isotope analyses, a spatial distribution map of groundwater samples was constructed (Figure 12), aiming to facilitate the identification and evo- This created a need for further study of the wider area to the east of the basin, in order to obtain a more satisfactory interpretation of this isotopic spatial distribution. About 18 km NE of monitoring point LB214 (where the δ 18 O value had the lowest negative value (−9.82‰)) is the Ossa mountain range. The Ossa range rises to 1978 m elevation. Ossa is mainly structured by pre alpine and alpine formations which, geotectonically, are integrated in three units ( Figure 12, [66]).
Regarding the tectonics of the area: according to [66], the most important transverse structures resulting from the statistical analysis of the Ossa tectonic elements are illustrated in Figure 12. These are: A, B, C, D, and E (the main fault zones); as well as F, and G (the secondary fault zones). The intense disruption observed along these zones (especially A, B and C, being the main fault zones), in combination with the anticlinal structure of Ossa, was a key factor in the formation of preferential groundwater flow paths to the eastern Thessaly basin (i.e., towards the study area). However, the rate of this recharge (and thus the exact contribution of this mechanism to the balance of the system) is yet to be determined.
The high altitude of Ossa carbonate formations (maximum altitude 1978 m) justified the existence of the lightest 18 O isotopes in the wells of the Pinios area in the eastern part of the basin, giving the spatial distribution of δ 18 O as shown in Figure 12. High negative values of δ 18 O, close to (−10), indicated water of meteoric origin. Considering the literature [67], for regions of Thessaly with similar characteristics, the composition of the precipitation indicated a recharge altitude greater than 1200 m. Hence, originating from higher elevations of the Ossa mountain, this recharge mound-in the form of lateral crossflows-does explain the observed δ 18 O spatial distribution anomalies and indeed the apparent anomaly in the spatial distribution of the hydrochemical water types that were earlier discussed.  Recharge: Talus cones and scree along the margins of the basin act as a favorable medium for the aquifer system's recharge, especially the extensive fluvial deposits at the exit of the Titarisios River to the basin. The aquifer system exhibits limited hydraulic interactions with the Titarisios River and the Pinios River. These var spatiotemporally in direction and magnitude [53]. The main recharge mechanism to the aquifer system is the lateral inflow across the western boundary with the karstic system of Tirnavos. Lateral

Conceptual Model of Groundwater Resources
Geometry: The Tirnavos subbasin is the northwesternmost part of the eastern Thessaly basin and is filled by alluvial sediments of varying thickness that reach a maximum of over 550 m towards its central parts [25,31]. The deposition environment progressively turns from fluvial/terrestrial to lacustrine toward the central parts of the basin. A sequence of confining units (which, on a regional scale, could be assumed to form a uniform layer of high clay content,) separated the alluvial sediments in a phreatic aquifer overlain by a thick and high-potential confined aquifer system ( Figure 13).
Recharge: Talus cones and scree along the margins of the basin act as a favorable medium for the aquifer system's recharge, especially the extensive fluvial deposits at the exit of the Titarisios River to the basin. The aquifer system exhibits limited hydraulic interactions with the Titarisios River and the Pinios River. These var spatiotemporally in direction and magnitude [53]. The main recharge mechanism to the aquifer system is the lateral inflow across the western boundary with the karstic system of Tirnavos. Lateral recharge from the karstic massif of Mt Ossa, at the northeastern extent of the eastern Thessaly basin, through preferential flow paths developed along the heavily tectonized and disrupted zones [66], is also a significant recharge source to the system. Deep percolation of precipitation forms a source of recharge to the system, and so do irrigation returns. The latter is of progressively reducing importance as more efficient and less water-consuming irrigation systems are being employed.
Discharge: Numerous production wells operate to cover irrigated agriculture and, to a lesser extent, domestic and industrial demands. These account for the majority of discharge from the aquifer system, which has led to the establishment of negative water budgets and the exhaustion of the phreatic aquifer. Natural discharge occurs in the form of lateral crossflows to the southeastern extension of the eastern Thessaly alluvial basin. Seasonal discharge occurs at temporally variable rates to the Pinios River upstream and the town of Larissa and the confluence of the Titarisios and Pinios Rivers downstream [53].  Recharge: Talus cones and scree along the margins of the basin act as a favorable medium for the aquifer system's recharge, especially the extensive fluvial deposits at the exit of the Titarisios River to the basin. The aquifer system exhibits limited hydraulic interactions with the Titarisios River and the Pinios River. These var spatiotemporally in direction and magnitude [53]. The main recharge mechanism to the aquifer system is the lateral inflow across the western boundary with the karstic system of Tirnavos. Lateral Flow domain: Dominant groundwater flow direction is from NW to SE, towards the extension of the eastern Thessaly basin. A key hydrodynamic feature in the basin is the development of a depression cone during the irrigation period, at the central parts of the subbasin, as a result of extensive groundwater abstractions. During hydrologically dry years, the depression cone expands and becomes sharper, characterized by higher hydraulic gradients [33]. Over prolonged droughts characterized by increased water demands, outflow to the SE extension of the basin reduces considerably and is even reported to be completely interrupted [25].
Quality characteristics & Land use impact: Excessive abstractions cause considerable groundwater heads decline and distort the spatial distribution and values of hydraulic gradients. Calcium carbonate is the dominant hydrochemical character attributed to the origin of recharge to the aquifer system. Regionally, residence times increase toward the central parts of the basin where ion exchange mechanisms become important, along with reducing redox conditions. Mixing of resident water with leachates of agricultural activityand possibly leaks of septic tanks and/or domestic effluent treatment plants-does bring in groundwater contamination issues. These predominantly take the form of elevated nutrients and sporadic findings of high Cu concentrations attributable mainly to orchards and cultivars. Small industrial unit activity is also reflected on groundwater quality in the form of elevated EC values.

Conclusions
Rational and sustainable management of groundwater resources requires knowledge of the geometric and hydraulic characteristics of the aquifer system that hosts the reserves, as well as a deep understanding of the hydrogeochemical and hydrodynamic mechanisms that control their evolution. This is even more true in cases of complex geotectonic environments, under intensive and prolonged exploitative conditions-especially within the climate change framework.
The alluvial system of the Tirnavos subbasin has a long record of exploitation, dating to the early 1970s, in support of the development of the region through irrigated agriculture. Anthropogenic activities have had a profound impact in shaping up groundwater reserves' availability and groundwater quality. Even though the Tirnavos subbasin is among the oldest and best-studied basins in Greece, this paper brings up new insights into its hydrodynamic evolution and controlling mechanisms. This enlightens and explains hydrogeochemical signatures and hydrodynamic behaviors which have been noticed in the past but disregarded.
The exercise performed in this research-characterizing and assessing the groundwater resources of the Tirnavos alluvial subbasin-consisted of a comprehensive overview of a suite of methodologies and approaches that acted complementarily to each other, leading to a clear and well-justified overview of flow and hydrochemical evolution mechanisms. Evidently, no single method is able to provide comprehensive, reliable answers or fully explain the observed spatiotemporal distribution of data. In this work, we have successfully made conjunctive use of hydrogeochemical, hydrodynamic, geological, structural, geotectonic and isotopic data-all analyzed in their spatiotemporal distribution, assessed in statistical indices, and evaluated in hydrogeochemical indices and tri-linear diagrams.
The key outcome of the present paper is the development of an integrated conceptual model for groundwater resource evolution, through the application of a specific methodological framework. This framework may be used globally as a generic model in other cases addressing similar scientific quests. The main findings of the research which compose this conceptual model, are briefly described below: - The dominant ions (Ca 2+ and HCO 3 − ) of groundwater were indicative for the main recharge mechanisms-which were related to the karstic substrate. - The recharge areas were delineated with the joint use of variable tools and concluded in two major axes in an E-W direction. The western direction was related to the recharge from the karstic system of Titarisios and the eastern one was related to the recharge through preferential flow paths from the karstic massif of Mt. Ossa. - The prevalence of carbonate formations was also reflected in the dominant hydrochemical type (Ca-Mg-HCO 3 ) which also encompassed the Mg content from relevant Mg-rich carbonate formations (dolostones). -Groundwater quality in the Tirnavos subbasin was in a generally good state, with few exceptions (elevated values of NO 3 and Cu) which indicate local impact due to anthropogenic activities. -Nitrates can be considered the main adverse environmental aspect, occurring locally in hot spots at the area, due to irrational agricultural activities. - The groundwater quality was also impacted by salinization due to the combined use of irrigation water returns (agricultural leachates) and evaporitic mineral leaching. - The governing hydrogeochemical process identified was ion exchange, which progressively alters the chemical composition of groundwater from west to east. -A secondary process which seems to affect hydrogeochemistry was redox, which locally controls speciation of groundwater solute parameters. - The physicochemical evolution mechanisms were jointly assessed and verified by the hydrochemical sections of major ions and saturation indices of critical mineralogical phases.
Even though the presented work sheds light on-and further proves the existence of-several major controlling evolution mechanisms, there is still room for further elaboration and in-depth analysis. This is important to capture potential alterations in the significance each identified mechanism might have on the system's evolution because of ongoing continuous changes in land and water use-and in the framework of climate change. Moreover, overexploitation of the phreatic aquifer-to the point of depletion and destruction or abandonment of shallow wells and boreholes-makes it imperative to seek out replacements for functional and aquifer unit-specific monitoring points, to ensure accurate, representative, and meaningful measurements are made. Moreover, the authors of the presented work propose further detailed analyses of redox processes and status. This will bring added value and knowledge and help to further decipher details of the controlling factors that shape the established hydrogeochemical signatures captured in the analyzed water samples.
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 in data processing, verified and optimized 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.

Acknowledgments:
The authors wish to acknowledge the valuable contribution of Hydrology Laboratory of Lubeck Technical University in analyzing the collected samples for isotopes-and especially C. Kuells for offering valuable guidance in the assessment of the results. This work was supported by the SWRI infrastructure, with the aid of which sampling and analysis of the collected samples was made possible. Special thanks are due to the laboratory staff of the Institute that performed the presented determinations. Last but not least, the technical and administrative staff of the Local Irrigation Organizations of Tirnavos, Ampelona and Agia Sofia are acknowledged for their support and guidance in the conducted sampling campaigns.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Table A1. Basic descriptive statistics of analyzed samples.