Vulnerability to Nitrate Occurrence in the Spring Waters of the Sila Massif (Calabria, Southern Italy)

Knowledge of spring waters’ chemical composition is paramount for both their use and their conservation. Vast surveys at the basin scale are required to define the nature and the location of the springs and to identify the hydrochemical facies of their aquifers. The present study aims to evaluate the hydrochemical facies and the vulnerability to nitrates of 59 springs falling in the Sila Massif in Calabria (southern Italy) and to identify their vulnerability through the analysis of physicochemical parameters and the use of the Langelier–Ludwig diagram. A spatial analysis was performed by the spline method. The results identified a mean value of 4.39 mg NO3−/L and a maximum value of 24 mg NO3−/L for nitrate pollution in the study area. Statistical analysis results showed that the increase in electrical conductivity follows the increase in alkalinity values, a correlation especially evident in the bicarbonate Ca-Mg waters and linked to the possibility of higher nitrate concentrations in springs. These analyses also showed that nitrate vulnerability is dependent on the geological setting of springs. Indeed, the Sila igneous–metamorphic batholith, often strongly affected by weathering processes, contributes to not buffering the nitrate impacts on aquifers. Conversely, anthropogenic activities, particularly fertilization practices, are key factors in groundwater vulnerability.


Introduction
In spring waters, nitrogen is one of the principal biogenous elements. It naturally occurs in the environment or derives from anthropogenic input, commonly found in fertilizers and animal and human wastes. Usually, human activities are altering all biogeochemical cycles [1,2] and tend to particularly affect the natural nitrogen cycle either directly (e.g., industrial, residential, agricultural, farming discharge) or indirectly, by altering soil degradation processes [3][4][5][6][7][8]. Furthermore, the dynamics of the hydrogeological systems and their water resource quality can be affected by climate change [9][10][11][12].
Once nitrogen is in the environment and converted to nitrate form, it will completely dissolve in water and move easily with water to aquatic ecosystems, where it can cause undesirable effects. Indeed, epidemiological studies evidenced that high levels of nitrate in water are a human health concern [13], considering that nitrate is classified as a probable human carcinogen [14]. Several studies have been conducted on the potential onset of stomach or brain cancer in people caused by exposure to nitrates through drinking water [7,[15][16][17]. More precisely, as a result of a reaction into an acid environment, when nitrates encounter the amines contained within food products, they produce nitrosamines, classified as carcinogenic substances [18]. Furthermore, nitrates have been demonstrated to be responsible for "blue baby syndrome" (methemoglobinemia) [19], especially in infants, whereby reduced intestinal acidity facilitates the proliferation of bacterial flora transforming nitrates The Sila Massif is part of a Paleogene orogenic belt (Calabrian Arc), thrusted over the Apennine Chain during Miocene times [27][28][29]. From a geological point of view, the area is a batholith formed by late Hercynian intrusive rocks and Paleozoic medium-to highgrade metamorphites, which are often strongly affected by weathering processes [30]: it underlies Mesozoic and Miocene to Quaternary sedimentary terranes [31,32]. The Paleozoic complex is composed mainly of paragneiss, biotite schists and grey phyllitic schists (with dominant quartz, chlorite and muscovite). The intense tectonics that molded the site, the subsequent and rapid uplifting affecting the area and an unfavorable set of climatic factors often reduced the crystalline rocks to sand [33,34].
The Quaternary is represented by Pliocene sediments, made up of light brown and red sands and gravels, blue-grey silty clays with silt interlayers, Pleistocene to Holocene alluvial sands and gravels and very small outcrops of Miocene carbonate rocks [35].
The dominant altitude of the central part of the Sila Massif is about 1300 m a.s.l., with the most representative peaks at almost 2000 m a.s.l.
Intense potato and cereal crops occupy the mountainous areas, while on the slopes, at the lower altitudes, there are pastures and forests of Pinus laricio (Pinus nigra subsp. laricio calabrica), many of which result from reforestation in the 1950s [36], and at the higher altitudes, beech (Fagus sylvatica L.) and silver fir (Abies alba) prevail.
Following the Köppen-Geiger classification [37], the climate of the study area can be identified as hot-summer Mediterranean, but due to its mountainous nature, colder winters with snow and cooler summer with some precipitation can be observed [38,39].
The study area is characterized by numerous springs with a density of 1.16/km 2 but of reduced flow rate (many of which have a flow rate lower than 1 L/s), typical of schistose-granitic settings which are almost impermeable.
In this survey, 59 spring waters were identified and sampled in the Sila Massif (Figure 1) together with 3 atmospheric deposition locations. For each site, an identification form was compiled accompanied by photographic material and contained all the initial The Sila Massif is part of a Paleogene orogenic belt (Calabrian Arc), thrusted over the Apennine Chain during Miocene times [27][28][29]. From a geological point of view, the area is a batholith formed by late Hercynian intrusive rocks and Paleozoic medium-to high-grade metamorphites, which are often strongly affected by weathering processes [30]: it underlies Mesozoic and Miocene to Quaternary sedimentary terranes [31,32]. The Paleozoic complex is composed mainly of paragneiss, biotite schists and grey phyllitic schists (with dominant quartz, chlorite and muscovite). The intense tectonics that molded the site, the subsequent and rapid uplifting affecting the area and an unfavorable set of climatic factors often reduced the crystalline rocks to sand [33,34].
The Quaternary is represented by Pliocene sediments, made up of light brown and red sands and gravels, blue-grey silty clays with silt interlayers, Pleistocene to Holocene alluvial sands and gravels and very small outcrops of Miocene carbonate rocks [35].
The dominant altitude of the central part of the Sila Massif is about 1300 m a.s.l., with the most representative peaks at almost 2000 m a.s.l.
Intense potato and cereal crops occupy the mountainous areas, while on the slopes, at the lower altitudes, there are pastures and forests of Pinus laricio (Pinus nigra subsp. laricio calabrica), many of which result from reforestation in the 1950s [36], and at the higher altitudes, beech (Fagus sylvatica L.) and silver fir (Abies alba) prevail.
Following the Köppen-Geiger classification [37], the climate of the study area can be identified as hot-summer Mediterranean, but due to its mountainous nature, colder winters with snow and cooler summer with some precipitation can be observed [38,39].
The study area is characterized by numerous springs with a density of 1.16/km 2 but of reduced flow rate (many of which have a flow rate lower than 1 L/s), typical of schistose-granitic settings which are almost impermeable.
In this survey, 59 spring waters were identified and sampled in the Sila Massif ( Figure 1) together with 3 atmospheric deposition locations. For each site, an identifi-Toxics 2022, 10, 137 4 of 18 cation form was compiled accompanied by photographic material and contained all the initial information collected. The sampling was carried out between the months of February and May 2016, a period in which the springs' maximum flow rates are generally recorded due to aquifer recharge following snow melting at high altitude and/or more simply following rainfall.

Laboratory Analysis
The samples collected were analyzed at the Genesis of Pollutants in the Water Cycle (GICA) laboratory at the University of Calabria. To evaluate the geochemical and hydrodynamic characteristics of the springs and in order to classify the investigated areas [40][41][42], a consolidated methodological procedure in geochemical approach was followed, analyzing the so-called labile parameters (temperature, pH and specific EC compensated at 20 • C) at the time of sampling by using portable apparatus including a pH meter and a conductivity meter and analyzing the main cations and anions once the samples reached the laboratory in clean high-density polyethylene (PE-HD) bottles with screw caps.
The major cations (Ca 2+ , Mg 2+ , Na + , K + ) were determined by using the atomic absorption spectroscopy (AAS) method, major anions (NO 3− , SO 4 2− , Cl − , HCO 3 − ) were determined by using the UV spectrophotometry method, alkalinity with determined by titration with 0.1 N HCl on unfiltered samples (expressed as mg HCO 3 − /L) and hardness was determined by using complexometric EDTA titration.
Total dissolved solutes (TDS) expressed in mg/L is calculated as the sum of all major anions (Cl − , NO 3 − , SO 4 2− and alkalinity as HCO 3 − ) and cations (Na + , K + , Mg 2+ and Ca 2+ ). A quality control analysis was performed evaluating sampling and analytical duplicates. Analytical precision was calculated by the 20% analysis (in duplicate) of randomly chosen samples and was found to be within accepted international standards (<5%).

Facies Classification Analysis
In literature, there are several methods to identify the dominant facies via a comparison of the most salient chemical characteristics, including by means of graphic analysis such as qualitative diagrams, which offer a good visual interpretation but no indication of the actual water mineralization, and quantitative diagrams, which, instead, provide an indication of mineralization [43].
The analyses carried out have allowed classifying the spring waters according to Italian Legislative Decree No. 31/2001, which categorizes waters as follows: (a) TDS (mg/L): • 0 < TDS < 50-very low mineral content water (or light mineral water); • 50 < TDS < 500-low mineral content water; • 500 < TDS < 1500-medium mineral content water; • TDS > 1500-waters rich in mineral salts. To determine the tendency of water to attack and solubilize certain minerals contained in rocks and soils, the aggressiveness index I A proposed by the American Water Works Association [44,45] was determined. This index, a simplified form of the Langelier Index, valid within the temperature range 4.5-26.5 • C, is expressed by the formula: where A is the total alkalinity (mg/L HCO 3 ) and H is the calcium hardness (mg/L CaCO 3 ). According to this index, the more aggressive the water is, the lower the value I A : • I A < 10-aggressive water; • 10 < I A <12-moderately aggressive water; • I A > 12-nonaggressive water.
In order to graphically visualize the water chemistry, various diagrams have been developed over the years [46,47]. Among others, in this paper, the hydrochemical facies were determined according to the Langelier-Ludwig (LL) square diagram [48], which has been largely applied in literature [49][50][51][52]. In particular, Langelier and Ludwig proposed a diagram in which rectangular coordinates are used to represent patterns and correlations between major cations and anions for multiple samples, thus allowing separation between Ca from Mg facies and Cl from SO 4 facies. More precisely, with the LL diagram, the water classification is ensured by positioning the water sample in one of the four sectors of which the diagram is composed: Chloride-sulfate Ca-Mg waters; II.
The identified ion concentrations are recalculated from mg/L to meq/L. Therefore, knowing all meq/L ion concentrations, the dimensionless value R can be calculated for each i-th cation C or anion A, according to the following relationships, and plotted on the LL diagram: The alignment of the sample points into the different sectors allows the identification of possible phenomena such as mixing or evolutionary processes of water [25].

Atmospheric Depositions Input Analysis
The presence of nitrogen compounds in the air is due both to natural causes (soil erosion, volcanic eruptions, etc.) and anthropogenic activities (fires, industrial plants, motor vehicles, heating, nitrogenous fertilizers). Nitric acid (HNO 3 ) in the atmosphere mainly results from the transformation of NO 2 reacting with free radicals. NO 2 in turn is formed by the reaction of NO with ozone. Nitric acid has a high deposition rate. In an aqueous medium, this compound turns into nitrite and nitrate ions. Nitric acid reacts with ammonia to form ammonium nitrate. Ammonium nitrate is in an equilibrium state with nitric acid and ammonia for relative humidity lower than 60%. Ammonium nitrate, which has low sedimentation, can be transported away from its originating source, and if humidity exceeds 60%, it turns back into nitrate and ammonia, contributing to the formation of acidity even in distant places.
To quantify the contribution of background nitrogen concentration from the atmosphere, total depositions at three different sites ( Figure 1

Data Analysis
A preliminary analysis of the dataset was performed through the box plot test in order to identify probable outliers.
Spatial interpolation is the procedure to predict the value of attributes at unobserved points within a study region using existing observations [53].
Spatial interpolation covers a variety of methods including trend surface models, Thiessen polygons, kernel density estimation, inverse distance weighting, splines and kriging. Among these several methods, mostly used in environmental analysis, in water Toxics 2022, 10, 137 6 of 18 resources management and in hydro-geochemical forecasting, the spline interpolation method was adopted for the analysis of the spatial variation of nitrate concentration in the spring waters of the Sila Massif. The spline interpolation model was chosen since the distribution of the sampled points did not allow a homogeneous coverage of the study area, and therefore the model that was least affected by the effect of inhomogeneity in the distribution and spacing between the points was applied. In fact, the spline method can generate sufficiently accurate surfaces from only a few sampled points, and the generated surfaces retain small features.
Spline interpolation is a deterministic interpolation method that fits a mathematical function through input data to create a smooth surface [54,55]. Splines are piecewise polynomial functions that are fitted together to provide a complete, yet complex, representation of the surface between the measurement points. Functions are fitted exactly to a small number of points while, at the same time, ensuring that the joins between different parts of the curve are continuous and have no disjunctions.
Spline interpolation creates a surface that passes through the control points and has the least possible changes in slope at all the points. This method uses a mathematical function with input points: where x and y are the coordinates of the point to be interpolated, and x i and y i are the coordinates of the control point. This function minimizes the overall surface curvature, resulting in a smooth curvature passing exactly through the input points.
A regularized spline yields a smooth surface and smooth first derivatives. With the regularized option, higher values used for the weight parameter produce smoother surfaces. The values entered for this parameter must be equal to or less than zero. Typical values that may be used are 0, 0.001, 0.01, 0.1 and 0.5. A tension spline is flatter than a regularized spline of the same sample points, forcing the estimate to stay closer to the sample data. It produces a surface more rigid according to the character of the modeled phenomenon. With the tension option, higher values entered for the weight parameter result in somewhat coarser surfaces, but surfaces that closely conform to the control points. The values entered must be equal to or greater than zero. The typical values are 0, 1, 5 and 10. The tension spline was evaluated as the best method for generating surfaces that vary gently, such as modeling of aquifer levels or groundwater pollution concentrations.
To determine the interpolation errors, the values of RMSE, %RMSE and MRE were assessed by Equation (6), Equation (7) and Equation (8), respectively. Values of RMSE and MRE close to zero indicate a good estimate of the model used to assess the unknown parameters: where x(P) is the estimated value of each component, x(m) is the measurement of water parameter, n is sample number, X is the mean of a measured parameter, z(x i ) is the observed value at location i, z*(x i ) is the interpolated value at location i and n is the sample size. The spatial analysis tool used for this analysis was ArcGIS 9.3.1 software.

Results and Discussion
A summary of physicochemical parameters of spring water samples and their statistics is provided in Table 1. The range of pH of spring waters is between 5.46 and 8.6 (mean = 6.86).
The mean water temperature of the analyzed springs is 10.51 • C with variability between 6-16 • C and a decreasing trend with altitude Z (m a.s.l.) according to the following correlation: The thermal gradient of −0.0052 • C/m indicates the presence of relatively shallow aquifers, influenced by the external air temperature. The spring temperature seems to be mainly influenced by vertical zonation because there are no statistically significant differences in temperature between springs located in different geological formations at the same altitude. In addition, the behavior of temperature vs. altitude is probably due to recharge from snow melting during springtime, which can promote water temperature decrease, especially in springs located at higher altitudes.
Electrical conductivity ranged between 51.3 and 710.61 µS/cm with a mean of 179.37 µS/cm.
The coefficient of variation (CV, %) representing the degree of scattering, showed a wide range among ions, with values for Na + (73.7%), Cl − (66.26%), SO 4 2− (216.62%), Ca 2+ (88.48%), Mg + (69.7%), K + (37.3%), HCO 3 − (90.34%) and NO 3 − (124.11%) ( Table 1). The major cation concentrations were detected for Ca 2+ and Na + , while the major anions ones were for HCO 3 − and Cl − . The nitrate (NO 3 − ) concentration had a mean value of 4.38 mg/L, and no sample exceeded the threshold value of 50 mg/L of WHO [21] and Italian Legislative Decree [22] recommendations. These results are in accordance with previous surveys, where the nitrate concentration did not generally exceed 2 mg/L in uncontaminated spring water aquifers [56,57]; anyway, in case of contamination, nitrates can reach extremely high levels. Nitrate pollution is increasing in many countries [58]. One example is in the Campania region (Italy), where nitrate values exceeding 200 mg/L in variable concentrations both in space and time have been found in spring waters and wells [59].
The spring water samples can be classified as very low mineral content (TDS < 50 mg/L) and low mineral content (50 < TDS < 500 mg/L), with a single sample that reaches the TDS value of 533 mg/L. These results indicate that the geochemistry of spring waters in the Sila Massif is strongly affected by the mineralogical composition of the local rocks and controlled by hydrolysis of sodium and/or calcium silicate minerals. The samples generally present an increase in low TDS with a decrease in altitude ( Figure 2a). As regards the hardness, a trend similar to that of TDS can be observed (Figure 2b). With few exceptions, the springs analyzed can be classified as light or soft waters, with a hardness value lower than 15 • f. There is a good correlation between hardness and TDS.
The spring water samples can be classified as very low mineral content (TDS < 50 mg/L) and low mineral content (50 < TDS < 500 mg/L), with a single sample that reaches the TDS value of 533 mg/L. These results indicate that the geochemistry of spring waters in the Sila Massif is strongly affected by the mineralogical composition of the local rocks and controlled by hydrolysis of sodium and/or calcium silicate minerals. The samples generally present an increase in low TDS with a decrease in altitude ( Figure 2a). As regards the hardness, a trend similar to that of TDS can be observed (Figure 2b). With few exceptions, the springs analyzed can be classified as light or soft waters, with a hardness value lower than 15 °f. There is a good correlation between hardness and TDS. According to the results of the aggressiveness index I A , it was possible to classify 78% of the analyzed springs as aggressive (I A < 10) and the remaining 22% as moderately aggressive (10 < I A < 12), the latter having a higher conductivity than the former. Figure 3 shows the I A classification of the spring waters. According to the results of the aggressiveness index IA, it was possible to classify 78% of the analyzed springs as aggressive (IA < 10) and the remaining 22% as moderately aggressive (10 < IA < 12), the latter having a higher conductivity than the former. Figure 3 shows the IA classification of the spring waters.  Aggressiveness is mainly due to the presence of carbon dioxide (CO 2 ), in the form of carbonic acid, as well as the presence of sulfates and nitrates which reduce the water pH. The acids, not being buffered by the acid rocks of the aquifers, act as a solubilizing agent on some minerals, which are washed out.
A correlation matrix of hydrochemical parameters was constructed and is shown in Table 2. Coefficients are significant at the 0.05 level (2-tailed), and those higher than 0.5 are shown in bold.
Because geological properties are closely related to characteristics of natural springs, including discharge, chemical composition, and temperature, hydrochemical parameter analysis between the occurrence of spring water and geological factors can be useful. Indeed, extensive monitoring campaigns are necessary, which can determine the natural cycle of spring waters resulting from the interaction between rocks, atmosphere, mixing with older aquifers and the effects of human activities. Groundwaters, because of their interaction with various matrices, become rich in gases such as carbon dioxide, oxygen and nitrogen; carbonic acid salts and strong bases, e.g., calcium, magnesium, sodium and potassium carbonates; and strong acid salts and strong bases, such as sulfates, chlorides, calcium, magnesium, sodium and potassium nitrates.
Accordingly, in this survey, major cations including Na + and Mg 2+ showed a linear relationship with alkalinity. This can be explained by the acidic hydrolysis of minerals. The hydrolysis reaction consumes water and acid that might have originated from CO 2 and increases the pH, alkalinity and cation concentration of water.
Further evaluations can be detected by observation of electrical conductivity values according to altitude for each zone (Figure 2c): although the altitude decreasing is combined with the conductivity increasing (with the exception of the three spring waters in which higher NO 3− concentration was detected), this last one remains low, thus evidencing a typical behavior of the igneous lithology and the presence of relatively shallow aquifers affected by atmospheric depositions [60].
The Langelier-Ludwig (LL) square plot, obtained by positioning the representative points of the individual springs within the four standard quadrants, is represented in Figure 4a.
to these facies, falling in metamorphic geological units located at the extreme offshoots of the investigated area, arising from metamorphic rocks. An exception is a spring in another setting (the municipality of Acri) classifiable within these facies due to its high content of sodium, chlorides, etc., although the aquifer rocks are granite like those of the IV quadrant. The waters are rich in sodium and potassium ions and show a low calcium and magnesium content. Alkalinity does not differ from that described for the I quadrant. According to their chemical compositions, the samples fall into three of the four quadrants of the LL diagram: • Chloride-sulfate Ca-Mg waters (I quadrant). The waters belonging to these facies mainly spring from the slopes forming the crown of the Sila Massif (Figure 4b). These springs fall into geological units made up of acid rocks (acid granulites, biotic gneisses). They have a low alkalinity which is almost independent of altitude. For them, good correlations between conductivity EC (µS/cm) and alkalinity A in terms of HCO 3 (mg/L) and between chloride and sodium ions were observed.
• Chloride-sulfate alkaline waters (II quadrant). A limited number of springs belong to these facies, falling in metamorphic geological units located at the extreme offshoots of the investigated area, arising from metamorphic rocks. An exception is a spring in another setting (the municipality of Acri) classifiable within these facies due to its high content of sodium, chlorides, etc., although the aquifer rocks are granite like those of the IV quadrant. The waters are rich in sodium and potassium ions and show a low calcium and magnesium content. Alkalinity does not differ from that described for the I quadrant. • Bicarbonate Ca-Mg waters (IV quadrant). They are the springs arising in the central part of the Sila Massif, from acid rocks (acid granulites, biotic gneisses, granites, granodiorites, magmatites). Alkalinity A increases with decrease in altitude Z. For these waters, an excellent correlation between conductivity EC (µS/cm) and alkalinity A was found (Figure 5a). • Bicarbonate Ca-Mg waters (IV quadrant). They are the springs arising in the central part of the Sila Massif, from acid rocks (acid granulites, biotic gneisses, granites, granodiorites, magmatites). Alkalinity A increases with decrease in altitude Z. For these waters, an excellent correlation between conductivity EC (µS/cm) and alkalinity A was found (Figure 5a). A good correlation between chloride and sodium ions was observed (Figure 5b). Waters with alkalinity close to 40 mg/L show a change in their gradient, becoming very similar to that of the I quadrant. A good correlation between chloride and sodium ions was observed (Figure 5b). Waters with alkalinity close to 40 mg/L show a change in their gradient, becoming very similar to that of the I quadrant.
The sampled waters located in the north and east of the Sila Massif, which fall in the facies of chloride-sulfate alkaline waters, show a NO 3 − concentration above 10 mg/L. This value is considerably higher than the other samples of about 2 mg/L found in springs unaffected by anthropogenic activities.
The increase in EC follows the increase in alkalinity values, a correlation especially evident in the bicarbonate Ca-Mg waters and linked to the possibility of higher nitrate concentrations in springs (Figure 5a).
Generally, the analyzed waters have a nitrate content lower than 5 mg/L, especially the bicarbonate Ca-Mg waters (IV quadrant), which show a mean NO 3 − content of 2.45 mg/L, as well as the chloride-sulfate alkaline waters (II quadrant) with a mean NO 3 − content of 1.8 mg/L. Instead, the chloride-sulfate Ca-Mg waters (I quadrant) are generally richer with a mean NO 3 − content of 4.8 mg/L with a maximum value of 13 mg/L. The results of statistical data analysis have detected the absence of outliers. As shown in Figure 5c, three water springs arising at hilly altitudes, falling in the municipalities of Acri and San Demetrio Corone, exhibit an anomalous behavior compared to the other analyzed water springs in terms of conductivity and presence of ions (Cl − , Na + , etc.). These three anomalous spring waters have nitrate concentrations ranging between 20 and 24 mg/L. These samples also present an anomalous content of the other analyzed parameters, such as a high conductivity, indicating probable aquifer contamination mainly due to widespread wild farming and pastures in the area, most practiced in the municipalities in which they fall. This factor could induce a greater exposure to anthropogenic contaminants and result in higher nitrate concentration. Table 3 shows the mean concentrations of NH 4 + , NO 2 − and NO 3 − recorded in the three sampling sites for the atmospheric deposition analysis. They allow us to define the extent of the atmospheric contribution into soil, considering that the concentration trend is seasonal with the maximum in spring. From previous studies, the nitrate concentration in the Bonis area (Sila Massif) is about 2 mg/L [61]; it shows lower concentrations than those found in the Crati valley (Sant'Antonello and Settimo) in Montalto Uffugo. Factors that may account for the higher deposition rates in these sites include higher rainfall amounts, local atmospheric inputs of N from volatilization of ammonia from fertilizers and animal wastes and N from the combustion of fossil fuels.
Cross-validation results for predicted versus measured nitrate values evidenced a good performance of the spline model in the study area. In fact, both MRE and RMSE values showed that the spline model can predict the data accurately. Figure 6 shows the spatial distribution of nitrates in the study area mapped through the spline method in five different ranges of values. The results were represented using the digital terrain model in the background. Toxics 2022, 10, x FOR PEER REVIEW 14 of 18 Figure 6. Spatial distribution of nitrates in the study area.
The obtained results were consistent with other studies that emphasized the spline method was one of the most suitable techniques for mapping nitrate concentrations in spring water. Indeed, the use of the spline method for water quality spatial analysis has been verified by numerous scientists. Guo et al. [62] investigated the five typical spatial interpolation methods, namely kriging, natural neighborhood, IDW, spline and trend surface, and found that the IDW and spline methods are appropriate for platform and small undulating areas. The national neighbor method and spline method provided the most accurate estimates of nitrates in the aquifer in the Qazvin plain [63]. The results of Zabihi et al. [64] indicate that the spline method is a good estimator of groundwater spring potential in the Iran area. The performance of the spline method was excellent in a case study in the Khalkhal region (Iran) [65].
The spatial distribution of nitrates provides some insight into the geochemical processes during the evolution of groundwater in the study area. The maximum level of nitrate concentrations was located in springs at lower altitudes in the north and western parts of the study area, in which lower precipitation and more intense agricultural activities occur compared to the eastern area. This factor could induce a greater exposure to anthropogenic contaminants and result in higher nitrate concentration.
The minimum level of nitrate, instead, can be attributed to the high-altitude area, which has no human activity and in which the lithological distribution may give an explanation for these values, as mentioned above.
In addition, in the Sila Massif, it is necessary to underline the subsequent reforestations since the postwar period and the related cleaning and maintenance activities of the woods that have deprived the subsoil of the organic substrate necessary for denitrifying bacterial colonies.
The study area was selected for its characteristics, along with the type of land use and human pressure, which are common to other wide Mediterranean areas, and therefore, despite its relatively small extension, it can well represent wider territories.
Moreover, the investigated area represents an extraordinary case of geological homogeneity, with the occurrence of purely acid metamorphic rocks. In it, there are few and The obtained results were consistent with other studies that emphasized the spline method was one of the most suitable techniques for mapping nitrate concentrations in spring water. Indeed, the use of the spline method for water quality spatial analysis has been verified by numerous scientists. Guo et al. [62] investigated the five typical spatial interpolation methods, namely kriging, natural neighborhood, IDW, spline and trend surface, and found that the IDW and spline methods are appropriate for platform and small undulating areas. The national neighbor method and spline method provided the most accurate estimates of nitrates in the aquifer in the Qazvin plain [63]. The results of Zabihi et al. [64] indicate that the spline method is a good estimator of groundwater spring potential in the Iran area. The performance of the spline method was excellent in a case study in the Khalkhal region (Iran) [65].
The spatial distribution of nitrates provides some insight into the geochemical processes during the evolution of groundwater in the study area. The maximum level of nitrate concentrations was located in springs at lower altitudes in the north and western parts of the study area, in which lower precipitation and more intense agricultural activities occur compared to the eastern area. This factor could induce a greater exposure to anthropogenic contaminants and result in higher nitrate concentration.
The minimum level of nitrate, instead, can be attributed to the high-altitude area, which has no human activity and in which the lithological distribution may give an explanation for these values, as mentioned above.
In addition, in the Sila Massif, it is necessary to underline the subsequent reforestations since the postwar period and the related cleaning and maintenance activities of the woods that have deprived the subsoil of the organic substrate necessary for denitrifying bacterial colonies.
The study area was selected for its characteristics, along with the type of land use and human pressure, which are common to other wide Mediterranean areas, and therefore, despite its relatively small extension, it can well represent wider territories.
Moreover, the investigated area represents an extraordinary case of geological homogeneity, with the occurrence of purely acid metamorphic rocks. In it, there are few and no significant anthropogenic presences, and the only inputs that contribute to the nitrate vulnerability in spring waters are attributable to the agricultural uses of soil and especially to forestry as well as to atmospheric depositions.
In addition, in the study area, investigations for vulnerability to nitrate occurrence in spring waters are scarce. The present survey confirmed the importance of assessing nitrate pollution as an environmental issue with multiple implications in terms of water quality deterioration, biodiversity loss, human health problems, global carbon-cycle alterations and climate change.

Conclusions
Springs perform numerous tasks, including ensuring sources of potable water and providing recreational, ecological and cultural value; moreover, they provide a way to assess groundwater quality since they are the connection between groundwater and surface water and give direct information about the state of groundwater in the aquifers that feed them.
The main focuses of this survey were to classify the spring water types of the Sila Massif (Calabria, southern Italy) according to their hydrogeochemical features and thus identify the factors controlling mineral waters using spatial variables focused on lithological settings and finally determine the vulnerability to nitrate occurrence in the spring waters.
The characterization of the waters of 59 springs in the study area has shown that they are strongly affected by the geological nature of the rocks of the aquifers (granites, magmatites, acid granulites, biotic gneiss, etc.). Indeed, the waters are generally poor in mineral salts (very low mineral content waters) with low calcium and magnesium content (soft).
The hydrochemical facies of the aquifers from which the water springs draw, according to the classification proposed by Langelier-Ludwig, can be classified as chloride-sulfate alkaline waters, those arising mostly from the mountainous slopes that crown the Sila Massif, and bicarbonate Ca-Mg waters that flow mostly in the central part of the Sila.
According to these considerations, most of the sampled points, having lower nitrate values, demonstrate the overall good quality of the spring waters in the Sila Massif. Therefore, the monitoring of springs in the Sila Massif revealed the scarce vulnerability to a potential alteration of their groundwaters by nitrates.