Hydrogeochemical Characters of Karst Aquifers in Central Italy and Relationship with Neotectonics

: Groundwater from karst circulation systems of Central Italy were sampled and analyzed, in 2018, for delineating a preliminary, general geochemical framework of their relationship with neotectonics, in an area characterized by a frequent and often destructive seismicity. We determined ﬁeld physical-chemical parameters, concentrations of main dissolved ions and gases and isotopic composition of water ( δ 18 O, δ D) and total dissolved inorganic carbon ( δ 13 C TDIC). We discriminated between “normal” hydro-karst systems and multi-component aquifers, composed of meteoric groundwater that have also interacted with rocks of di ﬀ erent lithological nature and / or deep ﬂuids. These multicomponent aquifers are of potential interest in the monitoring of neotectonics activity, because changes in the stress ﬁeld associated with the preparatory phase of an earthquake may a ﬀ ect the permeability of rocks, in turn causing variation of their chemical-isotopic character. The geographical distribution of these aquifers seems to be controlled by tectonics. In fact, the Olevano–Antrodoco–Sibillini thrust separates the more anomalous sites, located westwards of it, from the groundwater bodies at its eastern side, showing a more typical karst character.


Introduction
One of the most attractive challenges in geoscience is earthquake forecasting. Different approaches to this problem have been used in the last decades, based on the interpretation of variations and anomalies in geodetic, geophysical and geochemical parameters [1][2][3][4][5][6][7][8][9]. In particular, the relationship between seismogenic processes and fluid circulation has attracted the attention of several researchers, in term of the theoretical modelling of the physical processes responsible for chemical anomalies in geofluids induced by variations of crustal stress and strain [5,[10][11][12][13][14][15][16][17][18]. Both chemical characters and flow rates are sensitive to the seismic cycle and depend on the tectonic style [9,[19][20][21][22].
Karst areas are of relevant interest for the study and the monitoring of these processes because caves are geo-features whose genesis and development is inextricably linked to tectonic discontinuities that: (a) if of large dimensions and hydraulically active, concentrate groundwater drained from large rock volumes into single points (springs, pools and karst lakes and resurgences), which can offer hydro-geochemical information about large territories; (b) if accessible to humans, are good candidates for the installation of scientific instruments, relatively poorly disturbed by surface noise-generating processes [23].
Central Italy is an ideal place for studying the relationship between neotectonics and hydrogeochemical features of karst aquifers. In fact, this area ( Figure 1) is affected by a frequent and  [24], traces of main faults able to generate seismicity [25] and boundaries of the main karst areas [26].  [24], traces of main faults able to generate seismicity [25] and boundaries of the main karst areas [26].
Various studies carried out in Central Italy indagated the relationship between karst circulation and tectonic activity during recent years. Barberio et al. [9] monitored a set of springs in the area hit by the 2016Amatrice-Norcia seismic sequence, finding variations of pH, As, V and Fe concentrations months before the onset of seismic activity, and of Cr immediately after it. Barbieri et al. [27] and, previously, Petitta et al. [28], studied some springs and wells in the San Vittorino Plain, concluding that the variable geochemical character of groundwater is determined by the mixing of fresh water and deeper raising fluids, comprehending CO 2 released from a deep reservoir controlled by recent and/or active tectonics. A previous study by Favara et al. [18] reported on wide temporal chemical-physical variations in groundwaters with hydrothermal components during the 1997 Colfiorito seismic sequence, inferring that the observed anomalies had been induced by permeability variations related to crustal deformation, even in absence of seismic ruptures.
The relationship between geochemical transients in groundwater and neotectonics can be framed in the vaster conceptual scheme of the "crustal transients" proposed by Bernard [29]. According to this author, crustal transients are anomalous changes of physical-chemical parameters driven by unstable processes acting in the Earth's crust. Under this light and in our specific case, anomalous variations in geochemistry of groundwater can be considered 'fluid transients'. These do not necessarily precede or follow an earthquake, but both fluid transients and earthquakes are the effects of a more general instability, which can be individuated in the evolution of the stress field that leads to the generation of an earthquake (or not, in case of aseismic deformation). Changes in the stress field may affect permeability of rocks, causing in turn variation of parameters as flow rate and temperature of groundwater or its chemical composition [5,18,30].
Such anomalies could be particularly relevant if the affected aquifers are structured as multi-component circulation systems, and these components are characterized by significant differences in terms of permeability, chemical and/or isotopic composition, temperature, etc. If we consider the simplest case, e.g., a bi-component aquifer, the ideal simultaneous conditions for having potential efficient proxies of stress-related permeability changes are: where P 1 i and P 2 i are the values of the parameter i in the aquifer components 1 and 2, ϕ 1 and ϕ 2 the flow rates of the two components, and K 1 and K 2 their permeabilities. Of course, the necessary condition for making this system work is that the combined effect of inequalities (1), (2) and (3) leaves Pi over its instrumental detection limit. Such a system can be efficiently described using its radioelectric equivalent: a radio signal with a high-power carrier wave (the component with high flow rate and low concentration), modulated in amplitude by a low-power wave (the component with low flow rate and high concentration). Moreover, under condition (3), permeability changes driven by stress field variations should be negligible in component 2, as discussed by Madonia et al. [31], and consequently expected changes in flow rates should affect only component 1. Following these premises, the "fluid transient" ∆P 1 i /dt at the time tj is given by: As a practical example, we can consider a temperature transient in a bi-component aquifer, where a low-enthalpy hydrothermal member with T = 30 • C and ϕ = 1 Ls −1 mixes with a main meteoric member with T = 10 • C and ϕ = 1000 Ls −1 . If permeability variations in the hydrothermal member, induced by changes in the stress field, cause a ∆ϕ = 1 × Ls −1 , applying Equation (4) the resulting transient ∆T will be equal to 0.03 • C, e.g., a signal about tenfold the typical resolution of a commercial temperature logger.
These theoretical premises guided our work; we aimed to individuate aquifers where a main karst component is modulated by minor contributions, in liquid and/or gaseous phase, making these groundwater bodies potentially able to record, through fluid transients, variations in the local stress field.
Following these considerations, and with the aim of depicting a general geochemical framework of the possible relationship between karst groundwater circulation and neotectonics in the Central Apennines, we collected and analyzed water samples from karst areas with close spatial relationships with faults able to generate seismicity. We interpreted these data within the conceptual framework as described in the following chapter, with the aim of individuating sites of potential interest for the monitoring of neotectonic activity.

Study Area Settings
The central sector of Apennines, oriented NW-SE, is a part of the Apennines fold and thrust belt, built between Oligocene and Quaternary with a main eastward vergence. It is composed of Meso-Cenozoic sedimentary successions deposited over a Permian-Triassic phyllitic-quarzitic metamorphic basement. These successions include (from bottom to top) Triassic evaporites with anhydrites and dolomites and Meso-Cenozoic limestones and dolostones, topped by Plio-siliciclastic marine and continental deposits [9].
The intensely seismogenic extensional tectonic regime is controlled by NW-striking fault systems that has generated earthquakes up to Mw 7.0 (Avezzano earthquake, 1915), the last of which was the 2016 Amatrice-Norcia seismic sequence (highest Mw 6.5).
Huge fractured and karst aquifers are hosted by the Meso-Cenozoic carbonate sequences, often separated and sealed by clayely low-permeability layers, recharged by up to 1000 mm of precipitations (mainly rain and snow from the highest elevations of the mountain chain) and drained by springs often with impressive discharge rates (up to 18 m 3 ·s −1 ), characterized by quite steady regimes, with a depletion phase in late summer-early autumn times following a mild springtime increment [9].
Karst areas occupy a large part of the Central Apennines, and most of the main springs are related to karst circulation systems [26] (Figure 1).

Anaytical Methodologies
We collected 24 water samples from natural and artificial caves, wells, springs and pools, located in karst areas of Central Italy, in late spring 2018. The sampling campaign was carried out far (2 weeks) from the last rainy days, for avoiding the collection of samples not representative of the basal discharge but affected by shallower, fast infiltration. For the selection of the sampling sites we used a two-steps blind criterion. With the first step we selected the most important sites with respect their karst character, but completely ignoring their relevance with respect neotectonics. With the second step, we extracted from this list a subset of locations obeying to the criterion of close proximity to faults able to generate seismicity, but ignoring if these had been considered in previous geochemical studies aimed towards neotectonic investigations. In this way we ensured the statistical neutrality of the investigated data population.
Physicochemical parameters, such as temperature (±0.1 • C), pH (±0.1 U) and electrical conductivity (≤0.5%), were measured directly in the field, using Thermo Fischer portable instruments. We collected six different water samples at each location; four of these, for chemical and isotopic analyses, were stored in double-capped polyethylene bottles. Samples for determination of major ions were filtered using 0.45 µm Millipore MF filters and stored in separated bottles, acidifying that ones for cation determination with 100 µL of HNO 3 Suprapur. The other two bottles were used for isotopic analyses and HCO 3 content determination. Samples for analyses of dissolved gases and δ 13 C of total dissolved inorganic carbon (TDIC) were collected in glass bottles, hermetically capped with a crimper.
All the analyses were carried out in the laboratory facilities of INGV, Sezione di Palermo. Major ions were determined by ion chromatography, using a Thermo Dionex (Waltham, MA, USA) ICS 1100 instrument, with an accuracy better than 5%. Total alkalinity was determined by volumetric titration with 0.1 N HCl. The chemical composition of the dissolved gasses (He, H 2 , O 2 , N 2 , CH 4 , CO and CO 2 ) was determined according to the procedure by Capasso and Inguaggiato [32]; this method is based on the equilibrium between the liquid and the gas phases in the headspace of the bottles. We used an Agilent (Santa Clara, CA, USA) 78909 GC instrument, with Ar as gas carrier and TCD-FID as detectors; analytical accuracy was better than 5%.
The determination of δD and δ 18 O of water was performed by using a Thermo (Waltham, MA, USA) Delta Plus XP IRMS, coupled with an online pyrolysis system (TC/EA Thermal Conversion Elemental Analyzer) for the D/H ratio, and by CO 2 -water equilibration technique using a Thermo Delta V Plus instrument equipped with a Gas Bench II, for the 18 O/ 16 O ratio. The results were reported in δ% versus the V-SMOW standard, with a precision better than ±1% for δD and ±0.1% for δ 18 O. The isotopic composition of the Total Dissolved Inorganic Carbon (δ 13 C TDIC ) was determined with the method proposed by Capasso et al. [33], using a Thermo Delta V Plus IRMS coupled with an on-line gas preparation system Gas Bench II. The results were reported in δ% versus the V-PDB standard, with a precision better than ±0.1% . The correspondent values of δ 13 C of CO 2 in equilibrium with TDIC, used in graphical representations, were calculated using the formula proposed by Zhang et al. [34].
Field and laboratory data are summarized in Table 1; locations of sampling points are represented in Figure 1.

Results
A first, overall view on the geochemical facies characterizing the sampled groundwater is given by the Langelier-Ludwig diagram [35] represented in Figure 2, whose classification is based on four duplets of major ions, reporting to 50% the sum of concentrations (expressed as meq L −1 ) of anions and cations, respectively, obeying to the electroneutrality condition. As expected for karst groundwater, most of the points fall in the lower right corner of the lower right quadrant (bicarbonate-earth-alkaline water), dominated by carbonate species. Besides, two trends of mixing with different endmembers are evident. The first is a shift toward the composition typical of seawater, dominated by the duplets Cl-SO 4 as anions and Na-K as cations, which affect the samples STI (Stifone, in the cited work of Favara et al.) [18], ANB, ANV, VES and FOR. The other is an evolution toward the chlorine-sulphate-earth-alkaline composition (lower left quadrant), typical of evaporitic deposits or interaction with sulfur gas species, shown by samples SUS, CLI and SCI.
Further details about the mixing of karst water and other geochemical end-members are given by the ternary anion plot of Figure 3. Here, the relative abundances (based on concentrations in mg L −1 ) of Cl, SO 4 and HCO 3 are shown, allowing to investigate separately possible sources of Cl and SO 4 , not discernible in the Langelier-Ludwig diagram. Samples SCI, CLI and SUS follow a trend compatible with a mixing with a pure SO 4 source, while the points that in the previous plot of Figure 2 seemed to follow the mixing line with sea water occupy an intermediate region between the two evolutive trends.
the related groundwaters are oversaturated in relation to a certain mineral species, i.e., there are precipitation conditions and, vice versa, negative values indicate that dissolution occurs. As expected for karst waters, most of the samples are slightly oversaturated with respect to calcite and all undersaturated with respect to gypsum (more soluble), with a sequence from higher to lower undersaturation conditions following the concentration increments shown in the previous Figure 4; the sole exception is PFB, strongly undersaturated with respect to both phases.    Table 1. Summary of field and laboratory parameters. Longitude and latitude in decimal degrees WGS84, elevation in m a.s.l., type of site (artificial cave, spring, cave, well, pool), temperature in • C, pH, electric conductivity in µS cm −1 , major ions in mEq L −1 , alkalinity in mEq L −1 , δ 18 O and δD in % V-SMOW, δ 13 C of TDIC in % V-PDB, dissolved N 2 and CO 2 in cc L −1 STP, dissolved CH 4 in ppmv; nd is not determined, bdl is below detection limit. The possible role of the dissolution of gypsum/anhydrite as a source for SO 4 can be highlighted by the calculation of the parameter Ca*, that is obtained by:

ID
Once extracted, the residual Ca not linked to the dissolution of carbonates (calcite/dolomite) we plotted it in the binary diagram Ca*/SO 4 of Figure 4, where two different group of samples are evidenced. The first one, characterized by very low values of Ca* and SO 4 , includes all the samples from natural caves; the other one, related with higher values of both parameters, shows a good alignment following the line at slope 1, representative of the stoichiometric dissolution ratio of gypsum/anhydrite, especially for samples with concentrations higher than 2 meq L −1 .   We used the PHREEQC software [36] for calculating the saturation indexes of the main carbonate species (calcite) and gypsum, reporting the results in Figure 5. Positive values indicate that the related groundwaters are oversaturated in relation to a certain mineral species, i.e., there are precipitation conditions and, vice versa, negative values indicate that dissolution occurs. As expected for karst waters, most of the samples are slightly oversaturated with respect to calcite and all undersaturated with respect to gypsum (more soluble), with a sequence from higher to lower undersaturation conditions following the concentration increments shown in the previous Figure 4; the sole exception is PFB, strongly undersaturated with respect to both phases.
Composition of dissolved gases is reported in Figure 6, where the relative abundances (based on concentrations in mM L −1 ) of N 2 , CO 2 and CH 4 are shown; the composition of water in equilibrium with the atmosphere is also represented. Most of the points fall along the N 2 -CO 2 line, indicating a typical condition of groundwater, variably enriched in carbon dioxide during the infiltration through the soil or thanks to contributions of deep origin. Reducing conditions in the aquifers, mostly linked to degradation of organic matter or hydrothermal contributions, are evidenced by the presence of methane, which was measured in concentrations of up to several tens of ppmv (see Table 1) in fewer than a half of the samples.    Composition of dissolved gases is reported in Figure 6, where the relative abundances (based on concentrations in mM L −1 ) of N2, CO2 and CH4 are shown; the composition of water in equilibrium with the atmosphere is also represented. Most of the points fall along the N2-CO2 line, indicating a typical condition of groundwater, variably enriched in carbon dioxide during the infiltration through the soil or thanks to contributions of deep origin. Reducing conditions in the aquifers, mostly linked to degradation of organic matter or hydrothermal contributions, are evidenced by the presence of methane, which was measured in concentrations of up to several tens of ppmv (see Table 1) in fewer than a half of the samples. Isotopic composition of water is reported in the D versus δ 18 O diagram shown in Figure 7; 18 O values span from about −10 to −6 δ‰ (vs V-SMOW), which are typical for meteoric waters at these latitudes. The points fall in the area delimited by the lines representative of groundwater [37] and rainwater [38] of Central Italy. Isotopic composition of water is reported in the D versus δ 18 O diagram shown in Figure 7; 18 O values span from about −10 to −6 δ% (vs V-SMOW), which are typical for meteoric waters at these latitudes. The points fall in the area delimited by the lines representative of groundwater [37] and rainwater [38] of Central Italy.  [37] for groundwater and Longinelli and Selmo [38] for rainwater of Central Italy are also reported.
The most probable source of sulphate is the dissolution of the Triassic evaporites over which the carbonate sequence lay; this is demonstrated by the binary diagram SO4-Ca* (Figure 4), which indicates that calcium not linked to the dissolution of calcite and dolomite is in a stoichiometric relationship (1:1) with sulphate as in gypsum/anhydrite. The higher the concentration of both SO4 and Ca*, the better their correlation is, and at the opposite sides of the line there are water samples from caves (low SO4 and Ca*) and pools (high SO4 and Ca*), respectively. These differences reflect two distinct hydrogeological settings. Karst pools are located in the lower altitudinal belt (from 77 to 384 m a.s.l., see Table 1) and generated by the intersection between the topographic and the basal piezometric surfaces, involving deep aquifers seemingly hosted by both the older (and geometrically underlying) evaporites and the younger carbonate sequences. On the contrary, caves are in the higher altitudinal belt (from 727 to 1058 m a.s.l.); thus, samples here collected are representative of groundwater circulating in suspended karst aquifers (MAL, OMM, SFE) or related to direct infiltration through the rock cap over the caves (drip-water, as in BEA and CER).
Among the samples with a complex geochemical character, sites ANV and ANB deserve particular attention. They are both artificial tunnels carved for the exploitation of suspended aquifers in the subsoil of an old school (ANV) and of the famous Palace of the Pope Bonifacio 8th (ANB) in the historical city of Anagni (formerly known as Anagnia, in the Roman Age), located 65 km southeast of Rome (Figure 1). Anagni is located on a topographically isolated hill at about 400 m a.s.l., constituted of arenaceous turbidites [39] and elevating from the Sacco river plain, where travertines deposited after the quaternary volcanic cycle of the volcanic district of Central Italy diffusely outcrop. The particularity of these samples is given by the fact that remains of ancient Roman thermae (1st-2nd centuries) and groundwater catchment systems (as that sampled in ANB) have been found in Anagni, despite the small surface of the hill and its topographic isolation seeming incompatible with the presence of a groundwater body able to provide for the needs of a large human community.
The chemical characters of this groundwater (ANB and ANV in Figure 2) is intermediate between a pure carbonate term and an aquifer hosted in flyschoid sequences, which should be more shifted toward the sea water composition. This suggests that the Anagni groundwater body could  [37] for groundwater and Longinelli and Selmo [38] for rainwater of Central Italy are also reported.

Major Ion Chemistry
The first inferences come from the analysis of major anion concentrations, illustrated in Figures 2-5. The Langelier-Ludwig diagram ( Figure 2) remarks that, although most of the samples fall in its lower right corner, as expected for karst water dominated by the dissolution of carbonate species, some points move away from this geochemical fingerprint following two different mixing lines. The first trend evolves toward the composition typical of seawater (sites FOR, VES, ANV, ANB and STI), while the second is characterized by a progressive enrichment in chlorine and sulphate (SCI, CLI and SUS). The ternary anion diagram (Figure 3) allows discrimination between these two anions, indicating that trend 2 is entirely attributable to the delivery of increasing amounts of sulphate, whereas both chlorine and sulphate drive trend 1.
The most probable source of sulphate is the dissolution of the Triassic evaporites over which the carbonate sequence lay; this is demonstrated by the binary diagram SO 4 -Ca* (Figure 4), which indicates that calcium not linked to the dissolution of calcite and dolomite is in a stoichiometric relationship (1:1) with sulphate as in gypsum/anhydrite. The higher the concentration of both SO 4 and Ca*, the better their correlation is, and at the opposite sides of the line there are water samples from caves (low SO 4 and Ca*) and pools (high SO 4 and Ca*), respectively. These differences reflect two distinct hydrogeological settings. Karst pools are located in the lower altitudinal belt (from 77 to 384 m a.s.l., see Table 1) and generated by the intersection between the topographic and the basal piezometric surfaces, involving deep aquifers seemingly hosted by both the older (and geometrically underlying) evaporites and the younger carbonate sequences. On the contrary, caves are in the higher altitudinal belt (from 727 to 1058 m a.s.l.); thus, samples here collected are representative of groundwater circulating in suspended karst aquifers (MAL, OMM, SFE) or related to direct infiltration through the rock cap over the caves (drip-water, as in BEA and CER).
Among the samples with a complex geochemical character, sites ANV and ANB deserve particular attention. They are both artificial tunnels carved for the exploitation of suspended aquifers in the subsoil of an old school (ANV) and of the famous Palace of the Pope Bonifacio 8th (ANB) in the historical city of Anagni (formerly known as Anagnia, in the Roman Age), located 65 km south-east of Rome (Figure 1). Anagni is located on a topographically isolated hill at about 400 m a.s.l., constituted of arenaceous turbidites [39] and elevating from the Sacco river plain, where travertines deposited after the quaternary volcanic cycle of the volcanic district of Central Italy diffusely outcrop. The particularity of these samples is given by the fact that remains of ancient Roman thermae (1st-2nd centuries) and groundwater catchment systems (as that sampled in ANB) have been found in Anagni, despite the small surface of the hill and its topographic isolation seeming incompatible with the presence of a groundwater body able to provide for the needs of a large human community.
The chemical characters of this groundwater (ANB and ANV in Figure 2) is intermediate between a pure carbonate term and an aquifer hosted in flyschoid sequences, which should be more shifted toward the sea water composition. This suggests that the Anagni groundwater body could receive, through tectonic discontinuities, contributions from the underlying aquifer hosted in the Sacco river travertines, represented in Figure 4 by the TUF sample; the necessary hydraulic gradient should be provided by the lateral recharge of the aquifer, coming from the main elevations of this sector of the Apennines (Ernici Mounts, north-east of Anagni).
Further information is given by the calculation of the saturation indexes of mineral phases compatible with the chemical composition of the studied waters, among which those of calcite and gypsum, representative of the carbonate and evaporite components, respectively, are shown in Figure 5. As anticipated in the previous chapter, interesting peculiarities are found in BAG, CER and VIT samples. The spring VIT, located in the karst plain of S. Vittorino, at the foot of Mt. Velino, characterized by the presence of several thermal springs as Cotilia (sample COT of this work), spilled out following an earthquake at the beginning of the 19th century from the basement of the San Vittorino church, built between the 14th and 15th centuries. The spring caused the formation of a sink hole where the church has progressively collapsed.
Sample CER is from the Cervo (Deer) cave, very rich in speleothemes, previously sampled and studied for paleo-seismological reconstructions based on the analysis of their ruptures [40]. The intense concrection is an unquestionable clue about the general oversaturation of the drip water entering this cave, contrasting with the undersaturation of our sample. The sample was collected inside a small channel formed for chemical corrosion onto a previously deposited flowstone on the cave's floor (rejuvenation), indicating that newly formed fractures, compatible with a seismic shock, were introduced into the cave groundwater in disequilibrium with calcite. Signs of similar phenomena, but linked to the neo-circulation of groundwater oversaturated in calcite, are commonly found along the cave's wall, where fractures filled by calcite interrupt the continuity of previously deposited speleothems ( Figure 8).
Finally, sample BAG-T is from the same spring named "BAG" but collected by the operators of the local water service immediately after a Mw 3.9 earthquake, occurred with an epicenter located 13 km eastward from the spring on 21 May 2018, at an 8 km depth. As shown in Figure 5, the spring shifted to undersaturation conditions at the time of the earthquake, accompanied by a significant increment of water turbidity, as reported by the operators.
Another interesting site is the Scirca spring (SCI), the main resurgence of the Mt. Cucco karst complex, whose homonymous cave, with a prevalent vertical pace, reaches a maximum depth of more than 900 m. Although SCI is the basal spring of the karst complex, it does not present a mature chemical character for a carbonate water, remarked by saturation conditions with respect to calcite being close to zero and the minor undersaturation with respect to gypsum among the samples attributable to "typical" karst groundwaters. As in the case of the Bagnara spring, SCI also manifested an increment in turbidity after near seismic activity suggesting that, thanks to its composed chemical character, it could be a site deputed for investigating geochemical transients linked to seismicity.
The wide distribution of the points in the area of the gypsum saturation index comprised between −3.5 and −0.5 and reflects the analogous behavior of Figure 4, suggesting that the evaporite component is largely diffused in basal groundwater bodies of the studied area.
Vittorino church, built between the 14th and 15th centuries. The spring caused the formation of a sink hole where the church has progressively collapsed.
Sample CER is from the Cervo (Deer) cave, very rich in speleothemes, previously sampled and studied for paleo-seismological reconstructions based on the analysis of their ruptures [40]. The intense concrection is an unquestionable clue about the general oversaturation of the drip water entering this cave, contrasting with the undersaturation of our sample. The sample was collected inside a small channel formed for chemical corrosion onto a previously deposited flowstone on the cave's floor (rejuvenation), indicating that newly formed fractures, compatible with a seismic shock, were introduced into the cave groundwater in disequilibrium with calcite. Signs of similar phenomena, but linked to the neo-circulation of groundwater oversaturated in calcite, are commonly found along the cave's wall, where fractures filled by calcite interrupt the continuity of previously deposited speleothems (Figure 8). Finally, sample BAG-T is from the same spring named "BAG" but collected by the operators of the local water service immediately after a Mw 3.9 earthquake, occurred with an epicenter located 13 km eastward from the spring on 21 May 2018, at an 8 km depth. As shown in Figure 5, the spring shifted to undersaturation conditions at the time of the earthquake, accompanied by a significant increment of water turbidity, as reported by the operators.
Another interesting site is the Scirca spring (SCI), the main resurgence of the Mt. Cucco karst complex, whose homonymous cave, with a prevalent vertical pace, reaches a maximum depth of more than 900 m. Although SCI is the basal spring of the karst complex, it does not present a mature

Dissolved Gases
The chemistry of dissolved gases is illustrated in the ternary diagram of Figure 6, which reports the relative concentrations of N 2 , CO 2 and CH 4 (×10 5 ) and where the point representative of an air saturated water (ASW), falling very close to the nitrogen vertex, is also shown. In addition to nitrogen, coming from the atmosphere, a "normal" karst groundwater is more or less enriched in carbon dioxide (with respect to the ASW) deriving from soil respiration, thus points representative of these systems should align along the N 2 -CO 2 line. A little more than a half of our samples showed this behavior, with waters from caves that occupy the segment closer to the atmospheric composition and, conversely those with thermal contributions, as COT, VIT [28] and STI [18], or related to CO 2 -rich travertine aquifers, as TUF, falling closer to the carbon dioxide vertex.
It should be remembered [41] that in Central Italy huge amounts of non-volcanic CO 2 are produced by mantle degassing and/or metamorphic decarbonation. Carbon dioxide produced at depth accumulates in carbonate aquifers, from where it migrates toward the surface along fractures, faults or following morpho-tectonic structures (horsts), as its pressure (pCO 2 ) exceeds the hydrostatic one, giving rise to a diffused anomaly that characterizes most of the aquifers of the Tyrrhenian sector of Central Italy. The third main gaseous component, methane, is originated by the degradation of organic matter contained in sedimentary rocks, or released from hydrothermal systems (as in MIC, located in the thermal area of S. Vittorino Plain), under reducing conditions.
In conclusion, the gaseous dissolved fraction allows us to discriminate groups of samples with different characteristics, resulting from the mixing of meteoric and deeper components.

O, D, C Isotopic Composition
Oxygen and deuterium isotopic composition is illustrated in the binary graph of Figure 7. The points fall in the area delimited by the local groundwater [37] and rainwater [38] lines, indicating that meteoric (rain and snow) water is the main source recharging the aquifers. Isotopic composition of groundwater normally slightly differs from that one of the local rainfalls, due to the concomitance of several factors: complex circulation circuits collecting rainfall infiltrating at different altitudes, selective infiltration due to the seasonal cycle of evapotranspiration and/or to heavy rains generating surface run off but no infiltration and isotopic fractionation during surface run-off and/or infiltration through the shallowest soil horizons [42,43]. Other differences are caused by the inter annual variability of rain isotopic composition, up to several δ% for 18 O [42,43]. The maximum offset between the two lines and/or between the rainwater line and the points representative of our samples was about 1 δ% , namely largely comprised within the typical variability of these systems. Based on the abovementioned considerations, significant contributions from other feeding sources, as deeper hydrothermal circulation, should be excluded. Moreover, the similarity between our data and those from the work of Celico et al. [37], acquired in different sites more than 35 years ago, confirm that the samples we collected can be considered as representative of the average composition of the studied aquifers.
Differently to oxygen and deuterium, carbon isotopes give clear evidence of interactions between groundwater and carbon dioxide of both organic and deep origins. In Figure 9 we show a binary plot where the δ 13 C of CO 2 in equilibrium with the corresponding TDIC (see methodological chapter for details) is plotted versus the total carbon content. The combination of these parameters fixes some important constrains for interpreting our data; it is worth noting the negligibility of the atmospheric CO 2 in controlling the isotopic composition, because of its very small amount dissolved in groundwater at these conditions (about 0.02 mmol L −1 ).
Water 2020, 12, x 13 of 19 atmospheric CO2 in controlling the isotopic composition, because of its very small amount dissolved in groundwater at these conditions (about 0.02 mmol L −1 ). This plot clearly identifies two groups of data. The first one is populated by sites showing low carbon content, enriched in the 12 C isotope with respect to the points falling in the other group, which includes the sites with higher carbon content characterized by a higher abundance of the 13 C isotope; the isotopic divide between the two groups is at about −18 δ‰. Samples of the first group show no relationship between carbon content and its isotopic composition, and include sites that exhibit a dissolved ion content typical of karst systems (see Section 4.1): for these reasons we can define them as "normal" karst groundwaters. On the contrary, points of the second group exhibit a good, positive linear relationship between δ 13 CCO2 and total carbon content, with a determination coefficient of r 2 = 0.78 and a slope of 0.25 δ‰ (mmol L −1 ) −1 . Moreover, their chemical composition (see Sections 4.1 and 4.2) generally indicate a mixing between a meteoric component and fluids of different origins. These facts suggest that the isotopic composition of CO2 in this group of groundwaters is controlled by the extent of dissolution of carbon dioxide from a deep, non-meteoric source, isotopically compatible with a mantellic or decarbonation origin, as observed by Chiodini et al. [41]. This plot clearly identifies two groups of data. The first one is populated by sites showing low carbon content, enriched in the 12 C isotope with respect to the points falling in the other group, which includes the sites with higher carbon content characterized by a higher abundance of the 13 C isotope; the isotopic divide between the two groups is at about −18 δ% . Samples of the first group show no relationship between carbon content and its isotopic composition, and include sites that exhibit a dissolved ion content typical of karst systems (see Section 4.1): for these reasons we can define them as "normal" karst groundwaters. On the contrary, points of the second group exhibit a good, positive linear relationship between δ 13 C CO2 and total carbon content, with a determination coefficient of r 2 = 0.78 and a slope of 0.25 δ% (mmol L −1 ) −1 . Moreover, their chemical composition (see Sections 4.1 and 4.2) generally indicate a mixing between a meteoric component and fluids of different origins. These facts suggest that the isotopic composition of CO 2 in this group of groundwaters is controlled by the extent of dissolution of carbon dioxide from a deep, non-meteoric source, isotopically compatible with a mantellic or decarbonation origin, as observed by Chiodini et al. [41].
Sample ANB and ANV, although not properly definable as karst waters, have been included in the graph. They both show a geochemical fingerprint compatible with water-rock interactions with carbonates and flysch deposits (Figure 2), but ANV falls in the field of the "anomalous" samples because of its higher carbon content and enrichment in the heavier carbon isotope.

Potential Suitability of the Studied sites for the Monitoring of Neotectonic Activity and Its Spatial Distribution
As written in the introduction, the main aim of this work is the evaluation of karst aquifers of Central Italy as sites of potential interest for the monitoring of seismogenic processes, following the conceptual model shown in chapter 2. This model being based on geochemical transients generated by differential permeability variations, on their own driven by stress-field changes, the data we presented have been used for evaluating the potential susceptibility of a certain karst aquifer to be affected by detectable chemical-physical anomalies during seismogenic processes, ranking in a scale from 0 to 6 the studied sites according to this criterion ( Table 2).
The least susceptible (rank 0) are the "normal" karst aquifers, which are mono-component and expected to react to permeability changes generating physical rather than chemical transients, i.e., flow rate or temperature changes. This group includes the most of the caves and, more generally, aquifers whose aquitard is located at an elevation higher than the local base level, thus not favoring mixing with fluids of deeper origin. The only cave sample that occupies a better position in the ranking (1) is CER, due to its slight anomalous undersaturation in calcite ( Figure 5), which could have been caused by a seismically-induced rock fracturing ( Figure 8) and the consequent mixing of groundwater with different CO 2 partial pressures (the saturation curve of calcite as a function of pCO 2 is convex, so the mixing of two saturated end members can generate an intermediate, unsaturated groundwater). The following position (rank 2) is occupied by karst groundwaters that exhibit clues of mixing with a different geochemical end member (Figures 2-4), which is represented by the Triassic evaporitic deposits geometrically (and stratigraphically) underlying the carbonate sequences, and that could undergo physical-chemical variations during seismogenesis. Contrary to the previous group, this group includes sites located at relatively low elevations and thus close to the local basal circulation level of groundwater. Two of these sites, SCI and BAG, showed turbidity increments in coincidence with closer seismic activity, sometimes associated to changes in their chemical state (sample BAG-T, Figure 5). The next group (rank 3) is composed of sites (PFB, ANB and ANV) showing a contrasting geochemical character, i.e., strong chemical-isotopic anomalies but limited to few parameters. Sample PFB exhibits a "normal" geochemical karst character, in terms of its dissolved solid composition (Figures 2-4), but is characterized by a dissolved CO 2 with a deep carbon isotopic signature (Figure 9). Sites ANB and ANV, although hosted in non-karst terrains (arenaceous flysch) [39] , pertain to a suspended aquifer seemingly fed via tectonic discontinuities by an underlying carbonate aquifer, and show a complex chemical character making them of interest for the detection of possible seismically-induced geochemical transients.
The highest rankings (4 and 5) have been attributed to sites with evident both chemical and isotopic clues of mixing with fluids (water and gas) both of meteoric and deep origin. Rank 4 has been attributed to the karst pools CLI and SUS, which exhibit strong evidence of mixing between carbonate and evaporitic groundwaters ( Figure 4) and a dissolved CO 2 isotopically intermediate between the "normal" karst groundwater and the signature dominated by the contribution of deep fluids (Figure 9).
At the top of the ranking (5), there are the sites with strong signs of mixing between carbonate and evaporitic endmembers ( Figure 4) and an isotopic imprinting of the dissolved CO 2 dominated by non-organic, deep contributions ( Figure 9. Most of these sites are located in the hydrothermal area of the San Vittorino plain (COT, MIC, PES and VIT) and some of these showed, in historical times (VIT) and during the 2016 Amatrice-Norcia ( Figure 10) seismic sequence (site MIC), huge flow rate variations linked to a preferential vertical fluid migration path that also gave rise to sink hole phenomena (the collapse of the San Vittorino church in the 19th century).  However, since karst aquifers are complex natural systems that require holistic approaches in order to provide mutually supported results [45], future detailed hydro-structural investigations will Figure 10. Classification of the sampling sites based on their ranking for neotectonic activity monitoring purposes (see main text for details). Faults able to generate seismicity, the Olevano-Antrodoco-Sibillini (OAS) thrust and the epicenters of the most recent destructive earthquakes are also reported.
The spatial distribution of the ranks is illustrated in the conclusive Figure 10, where the epicenters of the last destructive earthquakes occurred in this area are also reported: the 1997 Mw 6.0 Colfiorito, the 2009 Mw 6.1 l'Aquila and the 2016 Mw 6.0 Amatrice-Norcia events. As shown in the figure, most of the sites with higher ranks are located westwards of the main structural feature of the Central Apennines: the Olevano-Antrodoco-Sibillini (OAS) thrust, which separates the Umbro-Marchean slope-basinal carbonate succession to the west from the Latium-Abruzzi carbonate platform deposits to the east [44]. Large scale tectonic features play a very important role in channeling toward the surface fluids of deep origin, as testified by the position of the San Vittorino Plain thermal area (sites COT, MIC, PES and VIT), located at the intersection between the OAS thrust and one of the main faults able to generate seismicity.
However, since karst aquifers are complex natural systems that require holistic approaches in order to provide mutually supported results [45], future detailed hydro-structural investigations will be important to better clarify the functional relationship between the OAS line and the studied aquifers.

Conclusions
As stated in the introduction and in the conceptual model description, this work aimed to investigate karst aquifers of Central Italy as potential sites for the geochemical monitoring of neotectonic activity. In doing this, we followed the leading idea that geochemical transients generated during seismogenic processes could more probably occur in multi-component groundwater bodies, where relative permeability variations could drive changes in mixing proportions, generating the transients. The analysis of chemical and isotopic data focused on the discrimination between "normal" mono-component hydro-karst systems and aquifers where the carbonate component is mixed, to various extents, with end members of diverse origin. These include both meteoric groundwater that interacted with rocks of different lithological natures, as the evaporite deposits at the base of the stratigraphic succession, and deep fluids with a more marked hydrothermal geochemical fingerprint. Previous studies in Central Italy had either focused on aquifers with evidence of mixing with deep fluids or restricted the analysis to single water bodies. Our research extended the field of investigation to the whole central sector of the Apennines, individuating as objects of interest all the multi-component aquifers able to generate geochemical transients driven by stress field changes, regardless of the origin (meteoric or deep crustal/mantellic) of the fluids mixing with the main karst component.
Data analysis and interpretation led to the conclusion that 15 out of the 23 studied sites exhibited characteristics of potential interest for the monitoring of neotectonic activity. One of these sites, the spring BAG located very close the epicenter of the 1997 Colfiorito destructive earthquake, gave us the opportunity to verify the existence of a geochemical transient in coincidence with closer seismic activity.
The spatial location of these sites seems correlated to tectonic lineaments: the main features were individuated in the OAS thrust, which separates the more anomalous sites, located westwards of it, from the groundwater bodies at its east side, showing a more typical karst character.
This preliminary geochemical mapping is a single shot on a system that evolves in a 4D space (x, y, z and time) and this is its main drawback. Further investigations should be based on the implementation of a network of automated stations in the 15 sites of potential interest, with the aim of recording geochemical parameters with the character of a potential earthquake precursor, as temperature, pH, Eh, electrical conductivity and dissolved CO 2 pressure. Funding: This research was realized on the behalf of the project "Centro di studio e monitoraggio dei rischi naturali dell'Italia centrale", funded by the "Fondo integrativo speciale per la ricerca (FISR)"