Towards the Understanding of Hydrogeochemical Seismic Responses in Karst Aquifers: A Retrospective Meta-Analysis Focused on the Apennines (Italy)

Earthquakes are known to affect groundwater properties, yet the mechanisms causing chemical and physical aquifer changes are still unclear. The Apennines mountain belt in Italy presents a rich literature of case studies documenting hydrogeochemical response to seismicity, due to the high frequency of seismic events and the presence of different regional aquifers in the area. In this study, we synthesize published data from the last 30 years in the Apennine region in order to shed light on the main mechanisms causing earthquake induced water changes. The results suggest the geologic and hydrologic setting specific to a given spring play an important role in spring response, as well as the timing of the observed response. In contrast to setting, the main focal mechanisms of earthquake and the distance between epicenter and the analyzed springs seems to present a minor role in defining the response. The analysis of different response variables, moreover, indicates that an important driver of change is the degassing of CO2, especially in thermal springs, whereas a rapid increase in solute concentration due to permeability enhancement is observable in different cold and shallow springs. These findings also leave open the debate regarding whether earthquake precursors can be recognized beyond site-specific responses. Such responses can be understood more comprehensively through the establishment of a regional long-term monitoring system and continuous harmonization of data and sampling strategies, achievable in the Apennine region through the set-up of a monitoring network.


Introduction
Major earthquakes (Mw 6.0 and above) are catastrophic events that can dramatically impact human beings, human infrastructure, and the environment [1]. Earthquake ground motion and ground rupture can affect the water chemistry and flow rates of aquifers in the epicentral areas and cause problems for emergency management, as well as longer term issues for water district managers [2,3]. The number and quality of case studies worldwide that have reported seismically induced groundwater anomalies (i.e., changes relative to background levels), discovered by sampling following major shocks [2,[4][5][6][7][8], has increased exponentially in the past decades. The role of fluids in the seismogenic processes is attracting the attention of the scientific community. Likewise, groundwater resources are vulnerable

Main Hydrogeochemical Features
In the Italian peninsula, the most important regional aquifers are located in the Apennines, along the Tyrrhenian margin of the mountain belt [39]. Central Italy, in particular, hosts many regional karst aquifers with high flow rates, reaching values up to 20 m 3 /s [38,39]. These aquifers are important sources of drinking water, feeding major metropolitan areas and bottling plants (Italy is the top country in Europe for production and per-capita consumption of bottled natural mineral waters [54]). Moreover, a high number of thermal springs, including thermal spas for health treatments (well the Apennines is about 500-1000 mm along the Tyrrhenian and the Adriatic coasts, while in the highest mountains along the central carbonate ridge it reaches values of 2000-2500 mm [37]. In the whole area, precipitation easily infiltrates into the Mesozoic limestone outcrops, which characterize large provinces of the Apennines landscape, generating an extensive and well-developed karstic circulation system. This results in a large number of karst springs, which are present throughout the mountain belt [36,38,39].  [40]. Numbers in legend indicate, respectively: (1) Adria plate foreland, mostly Mesozoic limestone areas; (2) foreland, Pliocene and Pleistocene deposits (delimited by the −1000 m isobath); (3) Po Plain and Adriatic foredeep domains characterized by a compressional tectonic regime; (4) thrust belt units accreted during the Alpine orogenesis in the Alps and in Corsica; (5) areas affected by extensional tectonics: these areas can be considered as a back-arc basin system developed in response to the eastward roll-back of the west-directed Apenninic subduction; (6) outcrops of crystalline basement (including metamorphic alpine units); (7) regions characterized by oceanic crust: newly formed  [40]. Numbers in legend indicate, respectively: (1) Adria plate foreland, mostly Mesozoic limestone areas; (2) foreland, Pliocene and Pleistocene deposits (delimited by the −1000 m isobath); (3) Po Plain and Adriatic foredeep domains characterized by a compressional tectonic regime; (4) thrust belt units accreted during the Alpine orogenesis in the Alps and in Corsica; (5) areas affected by extensional tectonics: these areas can be considered as a back-arc basin system developed in response to the eastward roll-back of the west-directed Apenninic subduction; (6) outcrops of crystalline basement (including metamorphic alpine units); (7) regions characterized by oceanic crust: newly formed oceanic crust has been recognized in the Provençal Basin (Miocene in age) and in the Tyrrhenian Sea (Plio-Pleistocene in age) whereas old Mesozoic oceanic crust can be inferred for the Ionian Basin; (8)

Regional Geology and Deformational History
The Apennines are mainly composed of Mesozoic-Cenozoic sedimentary cover derived from the Ligurian oceanic crust and the western passive margin of the Adriatic Plate, and of the Neogene-Pleistocene piggyback basin and foredeep deposits of the active margin [33,41]. The oldest rocks of Minerals 2020, 10, 1058 5 of 28 the exposed Meso-Cenozoic succession are 1.5-2 km thick Triassic shallow water carbonate platform and euxinic gypsum-anhydrite and dolostones with minor halite basin deposits (e.g., Burano Formation Evaporites [42,43]). Afterward, the shallow-water platform domain was broken-up during the lower Jurassic, permitting the development of new platform-basin systems, characterized by downthrown blocks, facilitating the accumulation of limestone, marly limestone and marl. In the platform environment, a late Cretaceous-early Miocene hiatus was followed by deposition of early Miocene paraconformable carbonates deposited along a carbonate ramp, while sedimentation was continuous during the Paleogene in the deeper environment. In the proto-western Mediterranean area, the tectonic regime changed during Eocene-Oligocene, when two subduction systems (the late stage SSE-directed Alpine subduction system and the young NW-dipping Apennines subduction system) became active along the retrobelt of the Alpine orogeny [40]. In this way, two subduction systems, with nearly opposite polarity, were present in a relatively narrow area for a short time (i.e., about 20 Ma [44]).
During the late Miocene, the southern Neotethyan passive margin became involved in the evolution of the Apennines, accreting passive margin sedimentary cover during the "eastward" roll-back of the NW-dipping Apennines subduction system [28]. The Apennines slab roll-back induced subsidence and boudinage of large portions of the Alps, which have been dismembered into the Apennines-related backarc basin: the Provencal and the Tyrrhenian basins ( Figure 2) [45].
Within the foreland basin setting, tectonically controlled sedimentary basins were developed. Hemi-pelagic marls, followed by turbiditic siliciclastic sandstones were deposited in the foredeep, ahead of the propagating deformational front [46,47]. North-eastward migration of thrust fronts developed different tectonic units that brought carbonate ridges onto turbiditic basins with NW-SE orientation [48]. Normal faulting appears to have been ongoing during the whole Quaternary and is still active today, producing one of the most seismically active provinces of the Mediterranean region ( Figure 1). The external (eastern) sector (Adriatic Domain) is characterized by normal heat flow, negative Bouguer gravity anomaly, a crustal thickness of about 35 km and a compressional seismicity characterized by hypocentral depths of 8-14 km. In contrast, the internal (western) sector (Tyrrhenian Domain) of the Apennines is characterized by thinned crust (20-25 km), high heat flow and CO 2 degassing, positive magnetic anomalies, shallow crustal earthquakes (6 to 10 km of focal depth), and positive gravity anomaly [39,49,50]. In this sector, in fact, the Apennine chain is presently affected by crustal delamination, with melts from the mantle wedge feeding active and quiescent volcanic systems; extensive alkaline K-rich volcanic systems are located in the Roman, the Neapolitan areas, as well as in Sicily [51,52], whereas silicic volcanic systems are present in Tuscan Magmatic Province [36,53].

Main Hydrogeochemical Features
In the Italian peninsula, the most important regional aquifers are located in the Apennines, along the Tyrrhenian margin of the mountain belt [39]. Central Italy, in particular, hosts many regional karst aquifers with high flow rates, reaching values up to 20 m 3 /s [38,39]. These aquifers are important sources of drinking water, feeding major metropolitan areas and bottling plants (Italy is the top country in Europe for production and per-capita consumption of bottled natural mineral waters [54]). Moreover, a high number of thermal springs, including thermal spas for health treatments (well documented since the Roman times) are present along the Apennine chain [55]. For these reasons, the area is well investigated and monitored regarding the chemical quality of the main regional aquifers, including both cold and thermal waters, which has allowed for documentation of transient and long-term hydrogeochemical changes during seismic crises.
Four main types of aquifers are present in the area: • karst carbonate litho-structural units, which can be considered the main reservoir of freshwater, characterized by a dense fault system and a widespread fracture-grid that supports infiltration of rainfall and circulation of water. The aquifers are generally hosted by Mesozoic formations made up of limestones, marly limestones and dolostones; • alluvial aquifers hosted by Upper Jurassic to Early Oligocene basin deposits of pelagic limestones, marly limestones, and siliciclastic deposits; • volcanic deposits, hosted in the volcanic complexes belonging to potassic series (e.g., Roman and Neapolitan areas) and silicic lavas (e.g., Tuscany); often geochemical anomalies regarding trace elements are reported in these aquifers (e.g., As, Li, U [53,55]); • evaporitic sequences aquifers, hosted in the Triassic dolostones and evaporites [56,57].
The tectonic setting of the area affects the spatial distribution and depth of the aquifers. There is, in fact, a discontinuity in aquifers from east to west. The eastern sector preserves the compressional structural features of the Adriatic Domain and the hydrogeology is characterized by the presence of vast recharge areas, with only a few important springs discharging most of the water flowing through the aquifer, causing similar temperatures of water in springs and in the infiltration areas [58]. In contrast, in the western sector, where the compressional belt is disrupted by important extensional structures of the Tyrrhenian Domain, the lower aquifer is partly buried under a low permeability cover (clays, shales, and marls) and is characterized by discontinuous recharge areas corresponding to the outcrops of the carbonate-evaporite formations [39]. Consequently, the western sector has a widespread occurrence of thermal waters with high pCO 2 values, connected with high regional CO 2 flux derived from a deep mantle related source [39,50]. Their composition was considered consistent with a mixing of two water bodies: a shallower component (carbonate component) equilibrated with calcite and dolomite and a deeper one (selenite component) circulating in the Evaporitic Triassic basement [55]. Hydrogeochemical facies of most of the springs, in fact, belong to the Ca-Mg-HCO 3 facies, while part of the collected waters presents a higher load in SO 4 , derived mainly from evaporitic aquifers ( Figure 3).
Minerals 2020, 10, x FOR PEER REVIEW 6 of 29 vast recharge areas, with only a few important springs discharging most of the water flowing through the aquifer, causing similar temperatures of water in springs and in the infiltration areas [58]. In contrast, in the western sector, where the compressional belt is disrupted by important extensional structures of the Tyrrhenian Domain, the lower aquifer is partly buried under a low permeability cover (clays, shales, and marls) and is characterized by discontinuous recharge areas corresponding to the outcrops of the carbonate-evaporite formations [39]. Consequently, the western sector has a widespread occurrence of thermal waters with high pCO2 values, connected with high regional CO2 flux derived from a deep mantle related source [39,50]. Their composition was considered consistent with a mixing of two water bodies: a shallower component (carbonate component) equilibrated with calcite and dolomite and a deeper one (selenite component) circulating in the Evaporitic Triassic basement [55]. Hydrogeochemical facies of most of the springs, in fact, belong to the Ca-Mg-HCO3 facies, while part of the collected waters presents a higher load in SO4, derived mainly from evaporitic aquifers ( Figure 3). Moreover, shallower alluvial aquifers are present in the Apennines, generally hosted in recent detrital alluvial and lacustrine deposits. These alluvial aquifers represent a hydrogeochemical facies similar to the karst one (generally slightly enriched in K and Na) but are sourced from shallower recharge areas, which can be demonstrated by water isotopes [59,60].
In this setting, such a variety of aquifer types are located in a relatively small and seismically active area. Therefore, geochemical markers of the different aquifers can be reliable tools to understand the causes of changes in water chemistry related to seismic activity (i.e., water mixing can be easily observed in the case of aquifer breaching).

Major Recent Earthquake Sequences
A significant cluster of strong seismic sequences has impacted the Apennine Chain over the last Moreover, shallower alluvial aquifers are present in the Apennines, generally hosted in recent detrital alluvial and lacustrine deposits. These alluvial aquifers represent a hydrogeochemical facies similar to the karst one (generally slightly enriched in K and Na) but are sourced from shallower recharge areas, which can be demonstrated by water isotopes [59,60].
In this setting, such a variety of aquifer types are located in a relatively small and seismically active area. Therefore, geochemical markers of the different aquifers can be reliable tools to understand the causes of changes in water chemistry related to seismic activity (i.e., water mixing can be easily observed in the case of aquifer breaching). October 2016, respectively, ruptured the Vettore fault system; all these events generated coseismic normal fault scarps. This sequence lasted for months, causing nearly 300 casualties, significant damage in numerous towns, and leaving more than 20,000 homeless [66][67][68].

Major Recent Earthquake Sequences
Moreover, various other events of significant magnitude were reported in historical accounts, which induced important hydrogeological effects. Hydrogeochemical changes of local aquifers during these historical events are reviewed in the section that follows.

An Overview on the Published Studies
Changes in hydrochemistry correlated with seismicity are reported worldwide in different tectonic and hydrogeological settings. Literature on chemical changes in karst aquifers is rich in case studies that largely cover single seismic sequences with sporadic samplings or short-term monitoring. For example, in the Pyrenees, changes in Cl ion and Pb isotope ratios were observed before a Ml 5.2 shock [21,70]. The increase of some major ions (specifically Cl and SO4) was also observed after an earthquake in the springs system of Dead Sea Rift Valley [71]. Changes in electrical conductivity (EC) were observed after a shock in the French Alps, as an index of permeability increase and water rock interaction enhancement [72]. Also in Turkey, after a Mw 7.4 earthquake, slight changes in major ion composition as well as an increase in flow and turbidity were observed [73]. These studies exemplify some of the hydrogeochemical responses that have happened in the Apennines of central Italy, which are described in more detail in the following paragraphs.

History of Observed Hydrogeochemical Responses in the Apennines
Focusing on the setting of the Apennines in Italy, abundant evidence of widespread seismically induced hydrogeochemical responses has been captured through various snapshot and response-

An Overview on the Published Studies
Changes in hydrochemistry correlated with seismicity are reported worldwide in different tectonic and hydrogeological settings. Literature on chemical changes in karst aquifers is rich in case studies that largely cover single seismic sequences with sporadic samplings or short-term monitoring. For example, in the Pyrenees, changes in Cl ion and Pb isotope ratios were observed before a Ml 5.2 shock [21,70]. The increase of some major ions (specifically Cl and SO 4 ) was also observed after an earthquake in the springs system of Dead Sea Rift Valley [71]. Changes in electrical conductivity (EC) were observed after a shock in the French Alps, as an index of permeability increase and water rock interaction enhancement [72]. Also in Turkey, after a Mw 7.4 earthquake, slight changes in major ion composition as well as an increase in flow and turbidity were observed [73]. These studies exemplify some of the hydrogeochemical responses that have happened in the Apennines of central Italy, which are described in more detail in the following paragraphs.

History of Observed Hydrogeochemical Responses in the Apennines
Focusing on the setting of the Apennines in Italy, abundant evidence of widespread seismically induced hydrogeochemical responses has been captured through various snapshot and response-based field campaigns. Hydrochemical anomalies include dissolved gas increases, changes in carbon isotope values, and increased trace and major element concentrations. Also, while this meta-analysis is mainly focused on the chemical changes, a large number of studies reported both transient and permanent changes in water level of aquifers, as well as changes in water flow from springs, which are also important observations with which to understand the mechanisms driving the changes (e.g., [6,74]).
Ancient observations in Italy on the effects derived from seismicity in water date back at least to the 1570 earthquake in Ferrara, located a few km east of the 2012 Emilia sequence in the foredeep of the northern Apennines (Adriatic Domain), when bubbling was noticed in various wells [75].
Nearly 200 years after these early observations, a scientific study by Robert Mallet in southern Italy noted the increase in turbidity and water flow in springs after the normal fault earthquake of 1857 in the Val d'Agri Quaternary extensional basin (Tyrrhenian Domain) [75,76].
In the last century, with the development of more precise instrumentation and better knowledge about water chemistry, the first observations of chemical changes correlated with seismicity were performed following the 1968, Mw 6.4, Belice earthquake in in the compressional foredeep setting of southwest Sicily (Adriatic Domain), where some major ions as well as well pH showed changes before and after the shock [77].
Most of the attention in the first studies of the 1960s and the 1970s was focused on the search for relationships between geochemical markers and tectonic or volcanic activity, in order to decipher possible precursory phenomena [77].
Hydrochemical changes were observed during the Ms 6.9, 1980 Irpinia normal surface faulting earthquake the southern Apennines (Tyrrhenian Domain), the largest event in Italy in the more than a century [78]. An anomalous increase in dissolved Rn concentrations in 2 springs in Latium was reported [25], which was interpreted as a result of dilation during pre-earthquake processes. In this case, the hydrochemical anomalies were observed in a location far from the epicentral area (more than 250 km away). Additionally, slight increases in major ions were reported in springs closer to the epicentral area [79].
The extensional tectonic setting of the southern Apennine belt was then struck again by a moderate seismic event (Val d'Agri, Mw 4.9) in 1996, during which 3 thermal springs were continuously monitored for physicochemical features [80]. Some anomalies were discovered regarding electrical conductivity and temperature. Moreover, Italiano et al. [81] reported changes in gas composition in one spring, indicating an increase in deep degassing.
A rich literature exists concerning the hydrogeochemical responses during the 1997 Colfiorito sequence. A number of thermal springs were monitored in order to understand geochemical sources of anomalies from mantle derived fluids. In fact, the main focus was set on bubbling and dissolved gases [11,57,82]. These results highlighted CO 2 degassing on a regional scale, consistent with the extensional tectonics of the Tyrrhenian Domain, which caused a drop in the pH of the thermal waters as well as changes in the mixing proportions between shallow and deep gaseous components. Some authors also reported changes in major ion composition of the same springs [56,[83][84][85][86]. Also, some cold springs in the same setting were analyzed for their major ions, and increases of HCO 3 and Fe concentration were detected [87].
In 2000-2001, in the compressional setting of Monferrato (northern Italy, Adriatic Domain), an anomalous increase in temperature of a shallow aquifer was reported after a series of moderate shocks (Nizza Monferrato, including two events of Mw 5.1 and 4.8). This anomaly was investigated through geochemical modeling and explained as enhanced leakage from a confined hot aquifer favored by an increase in pore pressure induced during a seismic swarm [88,89].
Throughout the 2009 L'Aquila sequence, several studies report changes in major ions. Two studies compared pre-and post-seismic values of major ions in different springs of the Gran Sasso karst aquifer system [74,90] and another study reported a decrease in biodiversity of macroinvertebrates in a series of springs in the SE part of the Gran Sasso Mountains, which was correlated with changes in major ions [91]. Temperature and water level also showed anomalies during and after the seismic sequence [92]. Moreover, changes in dissolved gas composition in several thermal springs during this seismic sequence were reported, as well as Rn anomalies in proximity to an activated fault [93,94]. A detailed regional study was performed by Chiodini et al. [95], measuring major ions, isotope ratios, and dissolved gases in various springs in the epicentral area. The authors pointed out the strong increase in radiogenic crustal gases (Ar and He), suggesting deep, long-term storage of high-pressure gas in the seismic area.
Additionally, in the deep area of the Gran Sasso Mountain massif, where hydrochemistry of groundwater is constantly monitored (due to the presence of the Gran Sasso National Laboratory, where nuclear physics experiments are conducted), an increase in uranium concentration was observed before the seismic swarm. These anomalies were highlighted as a possible increase in U derived from carbonate dissolution due to CO 2 release during the dilation phase before the earthquake mainshock [96,97].
In 2012, during the Emilia Romagna sequence, a preliminary observation of changes in chemistry of dissolved gases was presented by Italiano et al. [98] and then integrated with more hydrochemical data and a longer period of gas monitoring of the same deep wells, highlighting a notable increase in dissolved CH 4 and CO 2 [65]. Two other notable changes were identified. First, some changes in major ion concentration were observed in samples collected from bottling plants in the northern part of the Apennines [99]. Secondly, an anomaly in CO 2 content was noticed and continuously monitored in a thermal spring during this seismic sequence in the Tuscany region, about 100 km from the epicentral area [100].
During the last seismic sequence that struck the central part of the Apennine belt, in 2016 and 2017, there were 4 studies published with complete analyses of springs and wells in a range of 100 km from the epicentral area. Each study reported changes in trace element composition. Some studies noted anomalies before the mainshock, and others noted anomalies that occurred during the seismic sequence [101,102]. In contrast, another study analyzing springs in the same general region observed the abrupt increase in various trace metals following each of the mainshocks, followed by a decrease to values similar to the pre-earthquake baseline [2]. A follow-up paper by Barbieri et al. [12] then tried to clarify why some springs showed an effect from deep fluid interaction, interpreting the carbon and carbon isotopic data observed during this seismic sequence. Combining carbon data with Na/Li geothermometers, these authors observed that springs that showed responses in terms of carbon isotopes have a calculated CO 2 partial pressure at depth similar to that in deep boreholes in the northern Apennines. Thus, these springs are likely to indicate the effects of deep CO 2 mixing related with seismicity.
Building on these local and regional case studies, a few works have tried to compare more data, combining responses on different earthquake sequences. One study showed a long time series of inorganic carbon data following both 2009 L'Aquila and 2016 Amatrice earthquake sequences in a system of springs and wells around Rome [103], trying to find similar responses in different sequences. However, while an increase of carbon derived from deep CO 2 was observed before the 2009 L'Aquila mainshock, this anomaly was not detected for the 2016 Amatrice earthquake. Another recent study observed the correlation of CO 2 degassing and seismicity in a group of springs in the Gran Sasso area during two recent earthquake sequences (i.e., 2009 L'Aquila and 2016-2017 Norcia-Amatrice [104]). Finally, the first proposal to compare chemical data collected from environmental agencies to observe hydrochemical anomalies was made by Martinelli et al. [105].

Mechanisms Proposed for Earthquake-Related Responses
The mechanisms proposed to explain hydrological and chemical responses observed in previous studies can principally be classified into two main groups: • mechanisms originating external to the aquifer, including the intrusion of water from a different source; • mechanisms internal to the aquifer dynamics, with changes in water-rock interaction.
Mechanisms of the first group include expulsion of water from compressed aquifers [106], aquifer breaching and fluid mixing [14,21,107], interactions with deep fluids or gases [96,101], and consolidation and liquefaction of sediment [108]. Mechanisms of the second group are opening of deep fractures, which can increase chemical release with freshly created rock surface [3] and/or changes groundwater flow and aquifer properties to increase permeability [2,90,109].
A fundamental role in the understanding of the main mechanism at play is knowing the timing and the duration of the response, which gives clear information to discriminate mechanisms, that is, mechanisms immediately correlated with the seismic momentum increase (i.e., the microfracture opening or aquifer breaching) show a fast response after the mainshock [2,73]; mechanisms correlated with increased stress buildup can show anomalies also before the mainshocks, in the preparatory phase [97]. Additionally, the duration of the anomaly could give valuable information to discriminate the main mechanism causing changes ( Table 1).

Degassing and Interaction with Deep Fluids
Crustal degassing is observed to be strictly correlated with seismic and volcanic activity [95,111]. The gases dissolved in groundwater provide information on the water-gas interactions between water containing a dissolved atmospheric component and a deep-seated gas phase that occurs in the subsurface hydrologic system [11]. Among the analyzed gases, CO 2 is the main residual component of deep sourced gases, when temperatures drop and the water vapor condenses, and is the most common gas released from colder systems over tectonically active areas. CO 2 thus plays an important role for understanding the mechanism, especially for its presence in the aquifers of the Apennines, derived either from deep sources, or from carbonate rock dissolution in the aquifer enhanced through earthquake processes [112]. These separate sources of CO 2 can be distinguished by their different carbon isotope signatures, and can play a fundamental role in identifying internal and external mechanisms at play [2].
Along with the major gas components, changes of minor and trace gases (e.g., He, Rn, Ar) as well as their isotopic fractions (e.g., 3 He/ 4 He), can be used to evaluate mantle or crustal sources of the observed hydrogeochemical anomalies [81,83]. These gases, in fact, originate from the natural radioactive decay of elements in the mantle and crust (mostly U and Th series elements [84,113]). Therefore, their abrupt change, compared to the baseline atmospheric-derived values, is interpreted as an index of the influence of deeply-derived fluids on groundwater. This suggests that the increase in permeability of rock masses before, during, or after shocks can enhance the interaction of these gases with groundwater [57]. Moreover, minor and trace gases, especially Rn, were widely tested as possible precursory signals for earthquake forecasting [17,113,114].
Degassing can also affect changes in other variables through chemical reactions. CO 2 , for example, can mobilize redox sensitive metals from oxy-hydroxide phases, causing an increase in their concentration in the preparatory phase or after the mainshocks. This mechanism was inferred for trace element increase in groundwaters during the 2009 L'Aquila and 2016-2017 Amatrice sequences [2,96,97]. The effect of interaction or uprising of deep fluids could also directly affect chemical features of waters, typically causing the increase of dissolved trace metals (i.e., Li, B, As, U) [101]. This kind of CO 2 -driven mechanism could manifest itself in differential response timing. Generally, degassing could give a fast response due to the higher mobility of gases, whereas the mechanism involving deep thermal fluids should give a slower response.

Aquifer Breaching and Mixing
Seismic waves may enhance the permeability of the affected aquifers through increasing fracture porosity, changing pore pressure, and fracturing aquitards. These seismic effects could result in the mixing of fluids from different aquifers or the upwelling of deep fluids [109,115,116]. This effect was observed in various aquifers worldwide [14,107], but was proposed in only few cases in the Apennines (e.g., Piedmont area in the north [88,89]). In the case of aquifer breaching and mixing, timing of response can be faster than other mechanisms (e.g., release of deep-seated fluids), and the duration of the anomaly can be long (on the order of years) or even permanent. Hydrochemical variables that might be expected to change include the composition and concentrations of major ions, as well as water isotopes [107,115].

Mechanisms Internal to the Aquifer
The release of seismic energy can generate an increased dilation of existing fractures in the host rock, as well as in the creation of new fractures, which can serve to enhance permeability. This is especially the case in the high flow aquifers in central Italy (which are mostly karst). Seismicity increases permeability with the opening of new micro-fractures, the interconnection of pre-existing cracks and porosity (e.g., isolated pores in travertine), and the removal of infilling from cracks [2,7,70,90]. Evidence of increased permeability is increased flow rates observed after the seismic event, especially when there is no evidence of aquifer breaching described above (no change in major ions nor water isotopes) [90,117].
Moreover, the creation of new micro-fractures can easily expose fresh unweathered rock to the water flow, which can increase the release of major ions through increased water-rock interaction [3]. Another mechanism that is proposed, especially regarding karst springs, is the effect of fracture clearing, i.e., the release of highly mineralized pore water, which had been trapped in fractures and pores. Karst aquifers characteristically have multiple circulation paths (fast flow circuits on the scale of decades travelling along fractures and discontinuities, and much slower flow paths where water travels via seepage through lower permeability pore and fracture networks) [2,24]. Shaking would serve to expel the highly mineralized pore water from the slower circuit into fast flow circuit [2,91].
These internal changes are strongly affected by aquifer hydrogeological features (e.g., recharge area, flow rates and flow circuits, and bedrock geochemistry), which mainly affect the timing and the magnitude of responses. Chemical components ( i.e., major and trace elements) derived from within the aquifer system can show anomalies in the form of transient increases in concentration, but variables that indicate water mixing or interaction with deep fluids (e.g., water isotopes, dissolved gases, trace elements associated with deep fluids) are not expected to show changes.

Meta-Analysis: Methods
While the array of case studies discussed previously report hydrogeochemical changes related to seismicity, a synthesis and analysis of all these published data are missing in literature. Through the methods described in the following paragraphs, we attempted to combine the findings of previous studies to search for similar responses among the seismic sequences reported in the study area. We then inferred the most reliable mechanism to explain the responses.

Data Collection
Data were collected from published studies that took place during the last 4 major seismic sequences in the Apennines (1997 Colfiorito, 2009 L'Aquila, 2012 Emilia and 2016-2017 Amatrice) from tables, graphs, and supplementary material. For most of the studies, data were not available in electronic format as supplementary material, especially for the oldest papers. Therefore, data were collected manually from tables, where present, and from graphs using data extraction tools (e.g., WebPlotDigitizer, [118]).
We selected the studies containing at least 5 samples collected during the seismic sequence and in the pre-or post-earthquake period, in order to have a significant baseline value for the comparison with anomalies. Data collected only once from springs, without time trends were not included in the statistical analysis, nor were data reported only as mean values ± standard deviation of preand post-seismic sequences. As well, studies monitoring a single variable for long-term trends (e.g., long term continuous monitoring of U [96]), were not included in the multivariate statistical analysis, since they only increased the noise and complicated data interpretation. All of the remaining studies fitting these criteria were evaluated to determine the most plausible mechanisms explaining chemical changes.

Data Processing
Data were first organized by the spring locality sampled and earthquake sequence monitored. A few studies analyzed the same springs on different dates (and different sequences). In these cases, the times series were merged to observe possible differences in measurements from the different authors and to extend the available time span of monitored data. After the data were organized, all the records were normalized in order to emphasize the anomalies along the time series of data collected, as well as to remove differences in signals due to the high heterogeneity of the different sampling sites (i.e., differently mineralized waters). To accomplish this, a z score normalization of data was calculated following the equation (1): where z is the calculated score, x is the observed value, m is the mean of the whole sampling period and s is the standard deviation. This normalization easily allowed for the observation of anomalous values and outliers, reducing differences derived from different measuring units and different ranges of variance [119,120]. This calculation was performed for every spring selected in the meta-analysis to highlight the anomalies observed in comparison to what were considered background hydrogeochemical values [120,121], including main seasonal trends of each analyzed spring ( Figure 5). Normalized values for all the analyzed data were reported in Table S1, supplementary material. A Principal Component Analysis (PCA) was applied to the dataset, selecting the most frequently analyzed variables in the selected studies to observe their correlations in responses. PCA is a multivariate statistical technique that, through a linear combination of the analyzed variables, generates new uncorrelated variables named principal components. This approach allows observations of the data along the maximum variance of the system. In this way, the correlation of the measured variables (i.e., loadings) could be projected on the spaces defined by the principal components. The data obtained in the different samples can be projected as scores [122][123][124]. The normalization of data through z-scores was applied for every spring analyzed and highlights the anomalous values for each variable, compared with the geochemical background of different springs. For example, as observed for Fe and Cr (Figure 5), the anomaly reported during the Colfiorito sequence is more easily observable.  [87] and Amatrice sequence [2,101]. Black lines indicate the mainshocks of the two sequences (see Section 2.4). In this case, the effects of z-score transformation are observable, especially when analyzing the iron anomaly during the Colfiorito sequence, which is not very evident when analyzing the absolute values.

Variables Showing Changes and Their Correlation
As mentioned previously, Principal Component Analysis was applied to all the data records collected for the meta-analysis, but some variables (e.g., trace metals and radon) were removed because they were present in too few samples compared to the total number of collected data. Therefore, only physicochemical features, major ions (Ca, Mg, Na, K, SO4, HCO3, Cl), dissolved gases (He, CO2, CH4, O2, N2), and water isotopes were used for this analysis. Analyzing the loading plots of the PCA ( Figure 6A,B) shows which anomalies correlate in the different seismic sequences. The first three components were analyzed in more detail, because they explain about the 70% of the total variance of the system. Component 1 of the PCA, representing one third of the variance, shows a strong separation of most of the cations from the physicochemical parameters along axis 1. The water isotopes correlate with most of the cations along axis 1, suggesting a partial relationship with aquifers containing specific source waters and all cations but Ca. Moreover, the water isotopes correlate with CO2 and He, suggesting a possible role of gases in the chemical anomalies. There is also a correlation of most of the physicochemical parameters with the anions Cl and SO4 ( Figure 6A,B), elements that  [87] and Amatrice sequence [2,101]. Black lines indicate the mainshocks of the two sequences (see Section 2.4). In this case, the effects of z-score transformation are observable, especially when analyzing the iron anomaly during the Colfiorito sequence, which is not very evident when analyzing the absolute values.

Statistical Meta-Analysis: Results
Following the selection criteria described in the methods section, the springs/wells with long time series were selected, as well as the most commonly measured variables in the different studies. Among the literature analyzed, 1771 raw data records were collected for 63 sampling points, reported in 17 published papers.
These studies reported different suites of analyzed variables, which were not well distributed among the springs nor seismic sequences; in fact, different case studies reported different sets of variables. Consequently, we present several different outputs aiming to reveal the main correlations among variables and their response timing, as well as the interpretation of local effects and of tectonic features in the different earthquake sequences.
The normalization of data through z-scores was applied for every spring analyzed and highlights the anomalous values for each variable, compared with the geochemical background of different springs. For example, as observed for Fe and Cr (Figure 5), the anomaly reported during the Colfiorito sequence is more easily observable.

Variables Showing Changes and Their Correlation
As mentioned previously, Principal Component Analysis was applied to all the data records collected for the meta-analysis, but some variables (e.g., trace metals and radon) were removed because they were present in too few samples compared to the total number of collected data. Therefore, only physicochemical features, major ions (Ca, Mg, Na, K, SO 4 , HCO 3 , Cl), dissolved gases (He, CO 2 , CH 4 , O 2 , N 2 ), and water isotopes were used for this analysis. Analyzing the loading plots of the PCA ( Figure 6A,B) shows which anomalies correlate in the different seismic sequences. The first three components were analyzed in more detail, because they explain about the 70% of the total variance of the system. Component 1 of the PCA, representing one third of the variance, shows a strong separation of most of the cations from the physicochemical parameters along axis 1. The water isotopes correlate with most of the cations along axis 1, suggesting a partial relationship with aquifers containing specific source waters and all cations but Ca. Moreover, the water isotopes correlate with CO 2 and He, suggesting a possible role of gases in the chemical anomalies. There is also a correlation of most of the physicochemical parameters with the anions Cl and SO 4 ( Figure 6A,B), elements that are found in highest concentrations within Mesozoic evaporites at depth, indicating a possible release of these elements from evaporites during a seismic event. Dissolved carbonate and dissolved CO 2 correlate well with the water isotopes with respect to components 1 and 2 ( Figure 6A), suggesting that these properties are associated with processes affecting aquifers with specific source waters. In this way, component 1 is interpreted to indicate the effect of mixing with other water (for the negative values of water isotopes and some gases, which are inversely correlated to main physicochemical parameters), whereas component 2 indicates the role of increased ion concentrations derived from shaking (with high positive values of most of the cations). A plot of the loadings for the first and third components ( Figure 6B) shows a strong correlation of Ca and sulfates and may indicate a stronger influence with evaporitic rocks or brines in some aquifer settings. Moreover, the CO 2 and alkalinity results separate differently from principal components 1 and 2. Therefore, component 3 can explain the anomalies derived from deep fluid and gases, with positive values for CO 2 , He, and temperature. Still, this interpretation is limited by the fact that only some variables have sufficient data to allow for multivariate analysis.

Effect of Timing and Distance
Timing of anomalies and the distance of the springs where they are observed from the epicenters are often used to obtain information regarding the mechanism of hydrogeochemical responses [7,125]. Since most of the studies in this meta-analysis report an anomalous increase in solute during seismic sequences, Figure 7A reports the number of variables with z values greater than 2 and 3, compared with the number of days before/after the mainshocks (in order to align time series of the different earthquake sequences analyzed). As observed, there is a higher incidence of anomalous values closer to the mainshocks: most of the anomalous values reported (with z value higher than 3) occur less than 400 days following the mainshock and about half of the anomalies observed in the studies analyzed are reported during the 6 months following the first mainshock. This is the period mostly affected by major aftershocks in the last seismic sequences. These data highlight the effect of timing in defining the possible mechanisms because most of the anomalies are observable during a short period following the main shock, whereas a few other anomalies seem to happen after a long period.
The long-term anomalies are difficult to attribute to a particular factor and may be related to other features (e.g., climate). In addition, there are some springs which showed no anomalies. These features indicate that timing can be a helpful tool in defining anomalies but still need to be validated with other markers. Notably, there are more variables showing anomalies after the mainshocks than before, but this observation may be biased by the higher number of samples collected immediately after mainshocks, compared to before. Different springs are, in fact, analyzed with high sampling frequency only after mainshocks.
Comparing the anomalous values observed and the distances of sampling points from the epicenters of the mainshocks, which is often used to try to understand the effects of seismic energy release for groundwater changes [125], there is no clearly observable trend. Roughly about half of the anomalous values were observed less than 20 km from the epicenter of the mainshock, but other anomalies were observed farther than 80 km ( Figure 7B). This lack of correlation could reflect the fact that response is more affected by the hydrogeological features of the aquifer, such as the recharge area and the bedrock geochemistry, than by distance. Therefore, in this setting the distance from epicentral area seems to show little correlation to the hydrochemical response. Nonetheless, in the case of this meta-analysis, the maximum distance of sampling points was limited to <100 km, but a few studies reported in this review (e.g., [25]) report anomalies at even higher distances. It should also be noted that when considering all the recent earthquake sequences in the Apennines, other mainshocks happened within 40 km from the epicenter of the first shock, possibly affecting these results ( Figure 7B).
The analysis of different hypocenters was not reported here because the differences in this value among the sequences were generally less than 2 kilometers (most of the mainshocks were in the range 8-10 km below ground surface [32,61,64,67]).
As a side note, there is a relatively high number of Z values >3 during all seismic sequences, highlighting the fact that the geochemical anomalies related to earthquakes mostly cause a high change of chemical variables compared with the geochemical background of the analyzed aquifer. Comparing the anomalous values observed and the distances of sampling points from the epicenters of the mainshocks, which is often used to try to understand the effects of seismic energy release for groundwater changes [125], there is no clearly observable trend. Roughly about half of the anomalous values were observed less than 20 km from the epicenter of the mainshock, but other anomalies were observed farther than 80 km ( Figure 7B). This lack of correlation could reflect the fact that response is more affected by the hydrogeological features of the aquifer, such as the recharge area and the bedrock geochemistry, than by distance. Therefore, in this setting the distance from epicentral area seems to show little correlation to the hydrochemical response. Nonetheless, in the

Understanding the Main Mechanisms Causing Anomalies
The most promising and representative mechanism to explain the observed changes in the discussed data will be proposed in the following section, focusing on the chemical variables and the hydrogeological features that are most important in the investigation of hydrogeochemical responses. In this section, the meta-analysis of available data is merged with the interpretations proposed in the literature to explain the geochemical anomalies observed in the Apennines (also in the studies which could not be included in the meta-analysis), in order to infer the most plausible mechanisms at play.

The Focal Role of Hydrogeological Setting in Defining the Mechanism
The collection and analysis of different case studies allows us to highlight the most plausible and frequently occurring mechanisms for earthquake-related anomalies in the Apennines. The observed responses vary widely depending on water flow and the hydrologic and geologic setting of the recharge areas, whereas distance from epicentral area seems to be have little influence ( Figure 7B). This evidence confirms that the main hydrogeological features of the analyzed aquifer and the geochemical background levels need to be well known in order to understand the most plausible mechanism that leads to hydrochemical anomalies.
In thermal springs, which have a higher heat flux, the water chemistry is highly influenced by its gas component. It is also expected that these springs could have a higher interaction with deep fluids or gases (i.e., CO 2 , H 2 S), which could affect redox phenomena and dissolution patterns [12,101,126]. In fact, these springs (such as the springs analyzed during 1996 Colfiorito sequence [11]) showed changes in dissolved gas composition, as well as in their isotopic ratios (e.g., 3 He/ 4 He), as an indication of higher interaction with deep fluids [11,93]. These thermal-spring settings, in fact, seem more frequently influenced by mechanisms external to the aquifer during seismic crisis. Moreover, these anomalies can cause a fast response after the earthquake shock, as reported in our meta-analysis ( Figure 7A).
As observed in various studies of the four most recent seismic sequences, as well as from the PCA loading plot of our meta-analysis, CO 2 and carbonate equilibria play an important role in the hydrochemistry of the Apennines. Analysis of the plots of PCA loadings in Figure 6A,B indicates that they are, in fact, well correlated observing principal component 2, indicating an effect internal to the aquifer ( Figure 6A). However, these variables are not correlated (see principal component 3) ( Figure 6B), which could possibly indicate a different source of CO 2 from deep fluids (these results correlate with He, another deep sourced gas). In the PCA loading plot, CO 2 concentrations correlate with isotopes of water both on components 2 and 3. This correlation can indicate a mixing of other water sources correlated with changes in permeability after the release of seismic energy [6,109] as an explanation of the increase in dissolved CO 2 . These findings indicate that CO 2 can be a driver of groundwater changes correlated with seismicity, especially for springs with high heat flux (e.g., thermal hot springs) [65,95]. CO 2 , moreover, plays an important role in the interaction with other chemicals affecting water quality; this gas, in fact, could enhance the dissolution of redox-sensitive metals in water, an effect already proposed in some case studies [2,12,96], which could affect even cold shallow springs after the degassing caused by increased permeability during seismic activity [87,102]. It seems clear that deep-seated gas release is strongly affecting the chemistry of water in central Italy and this mechanism could be a valid hypothesis to explain the observed anomalies [12].
Shallower cold springs (karst or alluvial), on the other hand, mainly showed responses in major ions, as well as in trace elements, which generally show sharp changes in their concentration caused by mixing with stagnant fluids or release of elements from poorly bound exchange sites (e.g., Figure 5). This mechanism explains the fast response of the springs after the mainshock observed in our meta-analysis ( Figure 7A). In the case of karst and alluvial aquifers, moreover, the interaction with deep fluids is less likely to happen because of the thickness of low-permeability sediment strata of shallow recharge areas. The mechanism playing an important role in the setting of karst and alluvial aquifers (typical of the studied area) is the increase in permeability following the release of seismic energy, which causes the effect of fracture cleaning, as well as new exposure of fresh unweathered rock and the release of pore water into the fast flow circuit typical of these aquifers [2,26]. In this case the observed anomalies include a transient increase in the concentration of major ions and trace elements, which was observed widely during the last historical seismic sequences in the Apennines [2,90,91]. As observed in the PCA loading plot of components 1 and 3 ( Figure 6B), the increase of ions derived from dissolution of gypsum and brines (e.g., Cl and SO 4 , well correlated with EC in Figure 6) is the main chemical change indicating this type of anomaly, especially for the karst aquifers, presenting higher thickness and a possible enrichment in SO 4 and Cl (Figure 3, [2,58]). Moreover, in more recent studies (e.g., [2]), even trace element increases are reported as a proxy of this type of mechanism.
Following these hypotheses, variables such as dissolved gases and their isotopic ratio, major ions generally dissolving in karst aquifers, and trace metals are among the most reliable variables to evaluate the mechanisms causing changes in water chemistry. Trace metals were mostly analyzed in more recent studies but revealed important trends that were used to understand the mechanisms of hydrogeochemical anomalies and are of increasing interest in this research field [2,20,127].

Tectonic Setting and Hydrochemical Response
When evaluating the role of tectonic setting on hydrochemical responses, we saw common features, regardless of the setting. The tectonic setting of the central Apennines is associated with a dip-slip focal mechanism of earthquake formation, which should result in dilation. In the northern area, the focal mechanism causing earthquakes is associated with a strike slip setting, presenting less dilation (and consequently possible lower degassing and slightly increased permeability). The major changes observed during the 2012 Emilia earthquake (in the northern Apennines) are similar to the responses of the seismic sequences in the central Apennines: anomalous increase of dissolved gases especially in deep wells [65], presumably derived from deep-seated fluid interaction, and a slight increase in major ions observed in shallower karst aquifers [99]. The results in the 2012 Emilia earthquake are similar to anomalies observed in other localities with this tectonic setting, such as the increase of dissolved gases during the Colfiorito sequence [57], and the increase in major ions observed in shallower springs observed during Colfiorito, L'Aquila and Amatrice sequences [2,74,87].
While some authors indicated that earthquake-related chemical parameters can give us information to detect the focal mechanisms causing these responses [125,128], the results of this meta-analysis do not support that view. The 2012 Emilia sequence is the only one analyzed in this meta-analysis with a stike-slip focal mechanism, offering limited evidence to ensure a lack of differences.
While the direct effect of earthquake focal mechanisms on the hydrogeochemical responses to seismicity is not elucidated in this study, the hydrogeological setting is inevitably affected by the tectonic one. In this way, the different tectonic (and hydrogeologic) settings can present different types of responses to seismicity. For example, the high variability in CO 2 degassing sources changes with the setting [95,104,129]. During the last seismic sequences, the analysis of springs that have a large recharge area and can be affected by a deep source of CO 2 (typical of the Tyrrhenian domain) showed an evident increase in metal concentrations as well as fractionation of δ 13 C isotopes, indicating gas exchange from the deep CO 2 source [2,12,93,98,101,110]. Alternatively, aquifers that are mainly fed by shallow groundwater in the karst system (mainly present in the Adriatic domain) showed chemical changes caused directly by the increase in permeability and fracturing after the shock [2,90,101].

Knowledge Gaps and Future Perspectives for Earthquake Monitoring
Earthquakes are sudden, seemingly random events that make it difficult to collect long-term data for analysis that can be used to understand a complete picture of the pre-and post-shock responses. This is the main reason for the present situation, where a multitude of reported case studies have different unsolved issues regarding their nature (e.g., the role of pre-seismic anomalies, the length of the anomaly, or the heterogeneity in response among springs). In this section, the future trends for hydrogeochemical response analysis will be proposed, and the still debated issues of this discipline will be discussed.

Geochemical Precursors: A Research Horizon or A False Claim?
Earthquake forecasting is an attractive pursuit in the geological sciences that has been of interest since the early 1960s (Weichiuan earthquake), and this topic is still attracting considerable interest and generating contrasting opinions among the experts in this field [77,113,130].
Despite efforts to locate universal and reliable precursors, these signals are, at best, inconsistent and unreliable, because they do not always occur in the same places preceding major seismic events, do not always occur consistently within strongly affected areas (e.g., [2]), and can be generated from different mechanisms, making it difficult to find a systematic and reliable marker. Therefore, skepticism is high in the scientific community that reliable geochemical precursors can be determined [131,132]. Nonetheless, precursor signals are an area of research that has received considerable interest because of the potential benefit they could bring in disaster management [130,133].
During the period analyzed in this review, evidence of possible earthquake precursors was proposed in a few studies in central Italy. Precursory effects in U concentration were observed before the 2009 L'Aquila mainshock in the preparatory phase [96], and the authors of the study suggested that this element can be used as a potential indicator of pre-earthquake processes. In another study, an anomalous increase in a few trace metals was observed the month before the 2016-2017 Amatrice earthquake mainshock [101]. Also, an anomalous peak in carbon derived from a deep source was observed in a series of springs before the 2009 L'Aquila mainshock [103].
Analysis of the data collected in this review shows that pre-earthquake anomalies reported are fairly limited when compared with the observation of post-seismic effects. There is a sampling bias that makes pre-earthquake anomalies hard to evaluate, since the majority of observations occur after the seismic event has occurred. Consequently, this area of research needs further study and analysis in order to have more application, thus requiring a more robust pre-earthquake time series, including additional long-term discrete chemical data and potentially continuous data from a greater number of wells and springs. Such data may be able to shed light on the mechanisms that drive the changes, with the potential to use them as possible proxies of the preparatory pre-earthquake phase.

The Drawback of Data Availability and Data Analysis
The sporadic sampling of hydrochemical responses to seismic activity, mainly due to the lack of consistent and coordinated data collection, is one of the main obstacles preventing improved understanding of the mechanisms of earthquake induced anomalies. The length of the time series available and the quantification of background concentrations of elements are fundamental to avoid misleading results (e.g., false anomalies). The comparison with climate data (precipitation and temperature), flow data, and the observation of normal seasonal trends in chemical variables may be helpful in understanding how other mechanisms (e.g., aquifer breaching, deep seated fluid interaction) are changing water chemistry. In this way, harmonized and complete datasets are needed for a clear comprehension of these complex mechanisms.
As observed in this meta-analysis, multivariate statistical methods give valuable information about signal anomalies observed in data analysis, and data reduction techniques permit an understanding of which variables are causing the anomalies [2,3,56,134,135], as well as showing the timing during which variables are changing, which may elucidate mechanisms at play. In this review, seasonal cyclic changes and noise were removed using statistical tools (especially for long term monitoring with probes [92,100]). However, in some studies (e.g., [90,91]), average values of pre and post seismic data were used without considering the time span used for calculating the average. This kind of approach can complicate data interpretation, especially in areas that have marked seasonal trends, because the annual variance of data could blur the signal of the hydrogeochemical anomaly.
Moreover, while not deeply analyzed in this meta-analysis, continuous monitoring of physicochemical parameters could be helpful to observe the timing of responses with a high resolution (e.g., [41,136,137]). Nonetheless, the limited number of variables that could be monitored can lead to overlooking certain anomalies that may show up in snapshot sampling of a wider array of variables. A preliminary selection of chemical variables to monitor, based on the anomalies observed in the past, would be helpful and should ideally involve periodic samples of water isotopes, dissolved gasses, and selected ion chemistry

Towards Coordinated Data Collection to Clarify Mechanisms at Play
As observed in this review, consistent data harmonization can help to understand the main drivers causing hydrogeochemical anomalies related to seismicity; in fact, correlating similarities among different sporadic case studies helped identify coherent trends explaining the main mechanisms causing anomalies. Nonetheless, a joint and coordinated long-term data collection network would improve understand these phenomena [138].
The central Apennines, with an active recent history of earthquake sequences and a great number of regional aquifers providing public water supply is a plausible study area to set up a coordinated hydrochemical monitoring network. Additionally, given that our meta-analysis involved only 1 sequence from the Adriatic Domain, additional data from this tectonic setting could be helpful in ascertaining the role of focal mechanisms and stress regimes on hydrogeochemical response. With regards to the practicality of this endeavor, the monitoring of water chemical parameters is already in practice throughout the region, due to water standards monitoring in bottling plants, as well as thermal water plants, which were already analyzed during seismic sequences (e.g., [2,11,99]), providing a head start in establishing such a sampling network. Within the historical window that includes the studies in this review, the idea of a monitoring system was proposed in several articles and feasibility studies were reported [137,139,140], but a coordinated, long-term monitoring system has yet to be established. Recently, a preliminary group of research institutes started a collaboration with regional environmental agencies to set up a national monitoring system of hydrogeochemical variables connected with geological hazard observations [141], and if established, this monitoring network will be of great benefit in this endeavor, taking advantage of the water quality monitoring infrastructure already in place for bottling plants.

Conclusions
In this review, the Apennines of central Italy was chosen as a field laboratory to understand changes in water chemistry in response to seismicity, due to the high frequency of earthquakes and the high number of case studies reporting chemical anomalies. Following a meta-analysis based on studies from the last five historical seismic sequences in this region, the main variables showing changes were revealed.
The important role of response timing and duration was highlighted to elucidate the mechanisms causing changes. A differential response was often observed between the type of aquifers analyzed. Natural background values of springs are invaluable in determining hydrochemical response to seismicity and must capture natural variation, including seasonality, to avoid misinterpreting results. Interestingly, there was a common hydrochemical response that transgressed both distance of earthquakes from monitoring sites, as well as the focal mechanism of earthquake sequences and the tectonic setting at the spring site.
Two factors appear to play an important role in karst aquifers in this region; deep CO 2 release and changes in permeability, which has led to the hypothesis that two main mechanisms show hydrogeochemical anomalies: the release of deep gases derived external to the aquifer and the enhancement in permeability and fracturing due to earthquake shaking. These two mechanisms explain the different responses observed in different springs types, correlated also with the different hydrogeological settings of the areas. Nonetheless, open questions remain unsolved after these observations. Why are the anomalies generally not well distributed, with differences in responses? Could we obtain specific responses also from pre-earthquake data? A more complete understanding of mechanisms causing anomalies may come with a more complete, long-term dataset with harmonized data, in part generated by coordinated monitoring efforts. The creation of a monitoring network for data collection and analysis for the Apennines would significantly further our understanding of the relationships between seismic events and hydrogeochemical responses in aquifers.