Earthquake-Related Signals in Central Italy Detected by Hydrogeochemical and Satellite Techniques

Central Apennines are one of the highest seismic risk regions in Italy. A number of energetic events ( M W > 5) struck the region during the period 2004–2017, killing several hundreds of people (e.g., 294 casualties associated with the August 24th, 2016, M W 6.0 event of Amatrice). These earthquakes impacted piezometric levels, springs discharges, and groundwater chemistry across a large area, even at distances of dozens of kilometers from the epicenters. Here we present a multidisciplinary dataset based on hydrogeochemical and satellite observations associated with the seismic events that occurred in Central Italy during the period 2004–2017, which combines information derived from the application of groundwater monitoring and satellite techniques. Groundwater monitoring techniques allowed for the detection of hydrogeochemical anomalies in spring and well waters (14 water sampling points in total, with 22 variations larger than 2 σ ), while satellite techniques were applied to detect time-space variations in ground thermal emissions. We detected two significant, almost synchronous, anomalies in 2009 and 2016–2017 with both techniques, and we tentatively correlated them to crustal deformation processes. Part of the observed signals were detected before mainshocks, and they appear to be related to aseismic slip or to seismic slip eventually induced by minor fluctuations in seismicity. We argue that the combination of two factors, i.e., the shallow depth of local earthquakes and the concurrent deepening of groundwater circulation paths to several km depth, allow for the recording of variations in the stress field by geofluids released at the surface.


INTRODUCTION
Central Italy is affected by intense crustal deformation processes driven by the relative motion of the African and the Eurasian plates. In this region, active extension occurs along the axial and western zones of the Apennine mountainous belt (backarc extension area, characterized by crustal thinning and active volcanism toward the Tyrrhenian Sea). Conversely, compression mostly occurs along the frontal part of the belt, toward the Adriatic Sea (contraction area, associated with the eastward propagation of the Apennine accretionary prism). Strain values of 1.5-2 cm/year are currently estimated all over the extensive belt. An array of northwest-striking (i.e., northeast-and southwestdipping) normal and normal-oblique faults control regional seismicity, so that earthquake foci are prevalently distributed along the Apennines axis (e.g., Doglioni, 1991;Patacca and Scandone, 2001;Patacca et al., 2008;Montone and Mariucci, 2016).
Since the late 1970s, carbon dioxide emissions have been recognized to be associated on a global scale with areas of intense tectonic activity and/or extensive seismicity (e.g., Barnes et al., 1978;Irwin and Barnes, 1980;Tamburello et al., 2018). Large stresses associated with earthquakes and crustal deformations are known to cause measurable changes in the physical properties of rocks and significant hydrologic perturbations in the connected aquifers (e.g., Rojstaczer et al., 1995). Major hydrologic effects observed worldwide include spring and stream flow variations, groundwater level oscillations, triggered eruptions of mud volcanoes, gas emissions, modifications of the geochemical signature of spring and/or well waters (e.g., among many others, Wakita, 1975;King et al., 1981;Sato et al., 1986;Roeloffs, 1988;Roeloffs et al., 1989;Toutain and Baubron, 1999;King et al., 2006;Pierotti et al., 2015). Martinelli and Dadomo (2018) recently evidenced that a large part of possible geochemical precursory phenomena recorded worldwide have been found in areas characterized by i) shallow seismicity, ii) relatively high heat flux, and iii) significant crustal deformation rates. All these features are characteristics of the study area (e.g., Chiodini et al., 2004;Chiodini et al., 2011;Chiodini et al., 2020).
Moderate-magnitude earthquakes (4.0 < M W < 6.5) have been frequently recorded in Central Italy during historical and instrumental time. The hypocenters of these earthquakes were recognized in the upper crust, mainly at depths <15 km (Pace et al., 2006). The western part of the Apennine belt is also affected by intense CO 2 degassing. It has been estimated that several tons of CO 2 per day are discharged into the atmosphere via cold gas venting and by exsolution from high-P CO 2 thermal springs. An overall discharge rate of more than 20 Mt-CO 2 /year is associated to this area, roughly equivalent to one third of the total geological CO 2 degassed in Italy (e.g., Minissale, 1991;Chiodini et al., 2004;Minissale, 2004;Frezzotti et al., 2010;Chiodini et al., 2018).
Significant variations in the chemical composition of dissolved gases (e.g., Italiano et al., 2004), and in water level, temperature and chemical composition of water (Lapenna et al., 2004) were detected in thermal waters of Umbria region, Central Italy, during the 1997-1998 seismic sequence (Colfiorito M W 5.4 to 5.9 earthquakes). On April 6, 2009, a M W 6.3 earthquake rocked L'Aquila province. A seismic sequence including some 4 < M W < 5 events followed. Mainshock and aftershock effects upset deep fluids over a large area. Significant variations in the CO 2 degassing rate were detected in a gas vent located 80 km away from the epicenter (Heinicke et al., 2011). Other gas flow rate anomalies were detected on a regional scale by Bonfanti et al. (2012), and references therein) employing satellite-based techniques. Chiodini et al. (2011) reported an increase of deep originated CO 2 in local groundwater possibly connected with the L'Aquila seismic sequence. Nine more shocks characterized by 5 < M W < 6.5 occurred in Central Italy in 2016 and 2017 Festa et al., 2017, and references therein). Barberio et al. (2017) and Boschetti et al. (2019) reported water level and chemical composition variations in local groundwater, while De Luca et al. (2018) reported about groundwater pressure variations in some wells in the Abruzzo region before and during the 2016 Amatrice seismic sequence. These authors claimed that such variations "could potentially represent earthquake precursors".
Variations of the geomagnetic field (from satellites and observatories), and of atmospheric climatological data recorded during the 2015-2016 period have been retrospectively investigated to search for possible lithosphereatmosphere-ionosphere coupled effects during the preparatory phase of the Central Italy seismic sequence 2016-2017, and a few anomalies in electromagnetic parameters were identified over the year preceding the Amatrice earthquake of August 24th, 2016 (Marchetti et al., 2019). Satellite-based techniques (TIR index) captured possible intense and large scale degassing phenomena that occurred in concomitance with the 1997-1998 seismic sequence (Aliano et al., 2008). Further pieces of evidence of degassing during the 2009 seismic sequence were recorded by satellite techniques (TIR index; Genzano et al., 2009;Bonfanti et al., 2012;and references therein). Seismic sequences in Central Italy have been accompanied by crustal deformation processes (e.g., Cheloni et al., 2017). In some cases, the beginning of the crustal deformation processes was detected before the start of the sequences (Bella et al., 1995;Sgrigna and Malvezzi, 2003;Moro et al., 2017). Under these conditions, any geochemical precursor possibly detected by ground and/or satellite techniques, could be theoretically attributable to CO 2 -driven crustal permeability variations able to induce volume and shape changes in CO 2 reservoirs. These flux variations have the potential to modify the chemical composition of fluids hosted in deep and shallow aquifers. This mechanism is particularly effective in areas where natural manifestations (mostly mofettes) with high-CO 2 discharge rate occur (e.g., Chiodini et al., 2008; Figure 1).
In order to better understand the possible relations between local seismicity, groundwater chemistry, and CO 2 degassing activity, geochemical data recorded by several Regional Agencies for the Environmental Protection (ARPA) located in Central Italy were analyzed (the geographic coordinates of the water points of the ARPA network are given in Table 1 of the Supplementary Material). Furthermore, Thermal InfraRed spectral range radiations (TIR) attributable to fluctuations in CO 2 degassing rate (e.g., Tramutoli et al., 2013) were recorded by satellite techniques, and compared with other fluid-related parameters (pH, chemical composition of major constituents, P CO2 ). The multidisciplinary dataset used for this study spans the period 2004-2017, for which hydrogeochemical and TIR data were simultaneously available. The combined analysis of TIR and geochemical signals proposed here to indirectly evaluate variations in the stress field, is intrinsically innovative, and partially fills a methodological gap in the multidisciplinary approach to the study of earthquake precursors.

GEOLOGICAL AND SEISMOLOGICAL SETTING
The study area is located in the Central Apennines. The Apennine chain is the result of the opening of the Tyrrhenian sea that induced the eastward migration of a compressive front (e.g., Vai and Martini, 2001). Areas previously controlled by compressive tectonics were affected by back-arc extension due to the migration of the compressive front eastwards since the early Miocene. Crustal thinning due to extension in the Tyrrhenian area induced volcanic activities along the western margin of Italy. This region is characterized by high heat flow, positive gravity anomalies and shallow earthquakes due to the presence of a relatively thin crust (0-25 km). The eastern compressive Adriatic area is characterized by relatively low heat flow values, negative gravity anomalies and relatively high crustal thickness (30-35 km). Seismicity follows tectonic boundaries of distensive and compressive domains. All seismic sequences occurred in Central Apennines in past twenty years have been characterized by the occurrence of multiple events represented by M W > 5 mainshocks activating southwest-dipping normal fault. Each Central Apennine sequence may last 1-12 months . The most relevant seismic sequences occurred in Central Italy during the last 15 years are the destructive mainshocks that struck the L'Aquila province in 2009 (L'Aquila earthquake; M W 6.3) and the provinces of Rieti, Ascoli Piceno and Perugia in 2016 and 2017 (Amatrice earthquake; M W 6.0; Norcia earthquake, M W 6.5). Other minor seismic sequences characterized by lower magnitudes occurred in Umbria, Tuscany, Lazio and Marche regions during the same period ( Figure 1).
The main shock of the L'Aquila seismic event was preceded by a foreshock sequence that lasted four months, and was followed by an aftershock sequence that lasted eight months (Valoroso et al., 2013). The hypocenter of this earthquake was found along the Paganica-San Demetrio fault, a southwestward dipping fault segment intermediate between the eastern and the western major fault systems recognized in the Abruzzo sector of Central Apennines. The L'Aquila earthquake broke a ∼15-18 km long NW-SE trending fault, mainly confined to the upper 10 km of the crust (Serpelloni et al., 2012, and references therein). The area is characterized by a baseline heat flux of about 60 mW × m −2 , and by the absence of exploitable geothermal reservoirs.
The Amatrice earthquake started on August 24th, 2016, with a M W 6.0 mainshock apparently not preceded by foreshocks. Two months later, on October 26th, 2016, another M W 5.9 mainshock occurred 25 km to the north, and four days later, on October 30th, 2016, a M W 6.5 event (Norcia earthquake) hit the area in between the two previous quakes. Aftershocks sequence lasted about 1 year. Contrary to 2009 sequence, the 2016 seismic sequence apparently started without foreshock activity. A list of the major seismic events occurred in Central Italy during the observation period is given in Table 2 of the Supplementary Material.

METHODS
For the present study, we propose a novel approach (see Cicerone et al., 2009, for a review of previous monitoring techniques) that combines groundwater monitoring and satellite techniques to indirectly evaluate variations in the stress field, and possibly examine the effects due to the crustal deformation processes associated with such variations. Groundwater monitoring techniques were applied to detect hydrogeochemical anomalies (e.g., Gherardi and Pierotti, 2018;Martinelli and Dadomo, 2018), whereas satellite techniques were applied to detect time-space variations in thermal ground emissions (e.g., Tramutoli et al., 2013). The main advantage of the proposed approach is that it allows for i) identifying possible anomalies over large areas, thanks to the application of satellite techniques, and ii) for independently verifying these anomalies on a local scale, by ground-based inspection of the relevant geochemical parameters, with beneficial results for diagnostic accuracy. Hydrogeochemical monitoring consisted of the discrete monitoring of dissolved CO 2 on a large number of water points belonging to the ARPA (Regional Environmental Protection Agency of Italy) network (see Table 1 of the Supplementary Material). In particular, we processed a subset of 114 chemical analyses collected from springs and wells roughly distributed within and around the area of greatest damage caused by the 2016 seismic events, in the Umbria and Abruzzo regions. Finally, the MSG-SEVIRI (Meteosat Second Generation -Spinning Enhanced Visible and Infrared Imager) night-time TIR images collected over the Italian region during the 2004-2017 period were analyzed with the RST (Robust Satellite Techniques) methodology (Tramutoli, 1998;Tramutoli, 2005).

Robust Satellite Techniques Methodology and Robust Estimator of Thermal InfraRed Anomalies Index
The Earth's thermally emitted radiation measured by satellite sensors operating in the Thermal InfraRed (8-14 μm) spectral range (TIR) was considered in its apparent relation with the seismic events. The most recent literature (e.g., Tramutoli et al., 2013) now supports the hypothesis that a relationship could exist between positive TIR anomalies, and the release of volatiles originally stored at depth in rocks, as a result of stress field variations. Along with Rn, optically active gases like CO 2 and CH 4 are considered the gases most likely involved in this process. The areal extension and the geometric distribution of the thermal anomalies detected by the satellite likely also depends on the relative density of the emitted gas: the larger the amount of heavier-than-air gases present, the smaller the effect of atmospheric dispersion, and the more focused the detected thermal anomaly. On the reverse, the predominant release of gases lighter than air would produce diffusion-like patterns of the thermal anomalies. The RST (Robust Satellite Techniques) approach proposed by Tramutoli (1998), Tramutoli (2005), and Tramutoli (2007) was adopted to isolate possible TIR anomalies from other signal variations related to natural and/or observational factors that can be responsible of "false" proliferations. The methodology was already used to study the different seismic phases of earthquakes with a wide range of magnitudes (from 4.0 to 7.9), and occurred in different geotectonic settings (compressive, extensional and transcurrent). More details can be found in Tramutoli et al. (2015); Tables 3 and 4 of the Supplementary Material).
The RST methodology identifies space-time anomalies with respect to a preliminarily defined "normal" (i.e., under unperturbed conditions) signal behaviour, which is achievable by the analysis of long-term series of satellite records. TIR transients normally affect areas of several thousand square kilometers and can be detected up to several hundred kilometers far from the epicenter zone, some days to almost one month prior to a seismic event, and can last for a few days after the seismic event.
For seismic monitoring prone areas, TIR anomalies are identified by using the specific index RETIRA (Robust Estimator of TIR Anomalies; Filizzola et al., 2004;Tramutoli et al., 2005). The RETIRA index is computed on the image at hand only over locations not affected by cloudiness 1 , using the following equation: where: ΔT(x, y, t) T(x, y, t) − T(t) is the value of the difference between the punctual value of TIR brightness temperature T(x, y, t) measured at the location x, y and acquisition time t, and its spatial average T(t) is computed on the investigated area by only considering cloud-free locations, all belonging to the same, land or sea, class; μ ΔT (x, y) and Σ ΔT (x, y) are the time average and standard deviation values of ΔT(x, y, t) computed for the location x, y using only cloud free records acquired under similar observational conditions (same month of the year, same hour of the day).
The RETIRA index is intrinsically not protected by the proliferation of signal outliers, and TIR anomalies are subject to a space-time persistence analysis (e.g., Tramutoli et al., 2015). Such analysis is devoted to discriminating significant anomalies from well-known spurious effects (e.g., local warming due to cloud passages, error in image navigation/co-location process; see Aliano et al., 2008;Eleftheriou et al., 2016;Genzano et al., 2015).
Finally, it should be noted that by using the RETIRA index, the effects of natural (e.g., topography and land cover) and observational (e.g., solar illumination and satellite view angle) conditions are minimized, and the occasional warming due to daily variations and/or season time-drifts (i.e., inter-annual changes) are also reduced (see Tramutoli et al., 2005, for more details).
1 Cloudy pixels are identified by means the One channel Cloudy radiance detection Approach (OCA; Cuomo et al., 2004).

Discrete Monitoring of CO 2 in Selected Cold Aquifers of Central Italy
Two major areas of CO 2 earth degassing have been identified in Central-Southern Italy (Chiodini et al., 2004): i) the Tyrrhenian hinterland, where CO 2 is directly released into the atmosphere predominantly by two degassing structures, the so-called Tuscan-Roman and Campanian degassing structure, respectively; ii) the Adriatic foreland, where CO 2 is predominantly dissolved into groundwaters circulating through major carbonate aquifers hosted along the Apennine mountainous range. The centralsouthern part of Apennines is a highly active seismic area, and some authors (e.g., Chiodini et al., 2004;Chiodini et al., 2018) have speculated that the existence of high-pressure, CO 2 -rich fluid pockets at depth may play a major role in the generation of Apennine earthquakes. Groundwater circulation patterns in Central Italy have been described by Boni et al., (1986) and Petitta et al. (2018), and references therein), among others. From these studies, it emerged that groundwater hosted in the vadose zone is affected by fast vertical flow, while significant groundwater flow toward local springs with relatively high discharge values is generally driven by the relatively high thickness of the saturated zone. A number of springs, and wells 100 to 200 m deep, are currently monitored for environmental purposes by the Regional Environmental Protection Agencies (ARPA) of Central Italy. These manifestations are sampled for chemical analyses two to four times per year, according to the needs of local governments. All the samples have a calcium-bicarbonate chemical composition, low salinity, and outlet temperatures close to the average annual air temperature. The goal of the survey was to quantify the contents of CO 2 dissolved in the waters, and to compare any possible anomaly in CO 2 contents with seismic activity over time. Carbon dioxide fugacity values were calculated with the PHREEQC code (Parkhurst and Appelo, 2013) from the analytical concentration of total inorganic dissolved carbon and pH values measured in the field. Noteworthy, a quite high logP CO2 baseline value between -2 to -1.7 has been recognized all over the study area, and positive peaks of logP CO2 were observed roughly in the same areas during the seismic events of 2009 and 2016 (Table 5 of   Time series of recalculated logP CO2 values were compared with regional major seismic events (Figure 2). Two main groups of water points were identified: i) manifestations that sensitively responded to the M W 6.3 event of 2009 (largest in number); ii) manifestations that sensitively responded to the M W 6.0 event of 2016. Both pre-and post-seismic anomalies were observed, reflecting the intrinsic difficulty of extracting earthquakerelated signals unequivocally usable as seismic precursors from surface monitoring sites (Kumpel, 1992). Interestingly, the largest increases in CO 2 concentration were observed in the springs characterized by the highest P CO2 baseline values. Further to this, the largest anomalies associated with the seismic activity of Central Italy were observed in a well and a number of springs characterized by a large outflow of up to several hundred L/s. All these manifestations feed the Pescara river and drain the largest hydrostructure of the area, located in the Gran Sasso mountainous range. The occurrence of detectable CO 2 anomalies under such high water flow conditions was interpreted as additional, striking evidence concerning the large-scale increase in CO 2 degassing associated with the 2016 Amatrice earthquake. In the same area, similar phenomena were observed in different sampling points by Martini (2016) and by Chiodini et al. (2020).

Thermal InfraRed Anomalies
Daily TIR images collected at night-time 2 (time slot 00:00÷00:15 GMT) by MSG-SEVIRI satellite (sensor with spatial resolution of ∼3.5 × 3.5 km) over the Italian region in the period 2004-2017 were analyzed by means RST methodology and RETIRA index, with the aim to discern TIR anomalies possibly related to seismic process (Figures 3, 4).
FIGURE 3 | Example of sequence of TIR anomalies. The RETIRA index over Central Italy at the time of Amatrice earthquake (August 24, 2016, M W 6.0) (upper left box) was computed also considering low intensity anomalies (i.e., by considering pixels with ⊗ ΔT (x, y, t) > 2.5). Generally, low intensity anomalies follow high intensity anomalies, noticeably enlarging the anomaly area and filling gaps in both space (among isolated anomalous pixels) and time domains.
In order to discriminate significant anomalies from residual well-known spurious effects (outliers, geo-location errors, nighttime warm cloud passages and asymmetric distribution over the scene; see Eleftheriou et al., 2016, and references therein), in this work we considered to be significant anomalies those with pixels with ⊗ ΔT (x, y, t) > 3 that covered an appreciable spatial extension (at least 150 km 2 ), and that persistently appear in the study area (Central Italy). Two peaks of anomalous pixels with ⊗ ΔT (x, y, t) > 3 were detected over Northern-Central Italy during the 2004-2017 period (Figure 4) at the time of the Abruzzo earthquake (April 6, 2009, Mw ≈ 6.3) and of the seismic swarms that occurred in Central Italy in 2016 and 2017 (a further nine shocks with 5 < M W < 6.5 were recorded during this period).
An additional example of TIR anomaly sequence in Central Italy is reported in Figure 3. Here, the RETIRA index at the time of the Amatrice earthquake (August 24, 2009, M W 6.0) is shown based on low intensity TIR anomalies "anomalous" pixels with ⊗ ΔT (x, y, t) > 2.5). In general, low intensity anomalies follow high intensity anomalies, noticeably enlarging the anomaly area and filling the gaps in both space and time domains (e.g., Tramutoli et al., 2015, and references therein).

Cross-Checking of Independent Signals
A qualitative cross-check of the different signals considered in this study was finally done in Figure 4 by graphically superimposing geochemical (i.e., CO 2 ) patterns and TIR anomalies.
It emerged that a quite large number of springs/wells (14 water points in total, with 22 variations larger than 2σ; Table 6 of the Supplementary Material) from the ARPA network sensitively responded to both major seismic events that occurred in the Central Apennines during the 2004-2017 observation period (Figure 3). Associated with the main shock of the L'Aquila earthquake (2009), significant TIR anomalies were recorded almost concomitantly with P CO2 variations in groundwater samples from ARPA network.
Conversely, the Amatrice earthquake (2016) was largely anticipated by TIR anomalies observed up to several months before the event. In this case, the water points from the ARPA network registered an early increase in CO 2 concentration. The occurrence of possible precursory signals in Central Italy is not limited to degassing activity. Terakawa et al. (2017) reported that variations in local seismicity possibly induced by CO 2 pressure fluctuations at depth were observed before the 2009 L'Aquila mainshock, while crustal deformation processes were evidenced by GNSS before the 2016 mainshock (Panza et al., 2018). Furthermore, a possible precursory seismic pattern was observed in the same period and in the same area by Montuori et al. (2016) and by Gentili et al. (2017). Further possible precursory geophysical anomalies were also observed by Marchetti et al. (2019), Crespi et al. (2020) and Nardò et al. (2020).
The multidisciplinary approach presented in this manuscript has the potential to be applied to investigate seismic precursors outside the region considered here, in areas undergoing significant degassing, and characterized by the presence of surface water sampling points sensitive to crustal deformation processes (e.g., Martinelli and Albarello, 1997;Cioni et al., 2007;Gherardi and Pierotti, 2018). The key point to advance from a condition of mere a posteriori recognition of possible anomalies to the analysis of possible seismic precursors foresees the capability to detect the real temporal variability of the natural system under observation, and requires that relevant geochemical parameters should be strictly monitored continuously, as done elsewhere (e.g., in Tuscany, Italy; Pierotti et al., 2015;Pierotti et al., 2017). Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 584716 CONCLUSION Geochemical effects of crustal deformation processes on geofluids in Central Italy have been evidenced by ground and satellite techniques. Geochemical sampling in thermal springs and in gaseous emissions has allowed us to identify large-scale degassing phenomena in concomitance with the seismic sequence that occurred in Central Italy in 2009. The time series analysis of P CO2 values in groundwaters showed multiple degassing episodes in concomitance with seismic sequences that occurred in 2009 and 2016-2017. These large-scale CO 2 degassing phenomena have also been evidenced by satellite techniques capable of identifying transient temperature anomalies due to eventual temporary greenhouse gas emissions. Part of the anomalous signals were recorded by both techniques before the most relevant mainshocks. The contemporary observation and recording of geophysical and geochemical parameters emerged as an efficient method to gain a better understanding of seismogenic features of Central Italy. This approach has the potential to be successfully applied elsewhere, in areas of active crustal deformation, to possibly identify geophysical and geochemical signals before relevant seismic events. The continuous monitoring of relevant geochemical parameters from properly selected water manifestations is highly recommended to ensure the capability to detect signal variations occurring over short time intervals, and possibly improve the performance of the method in terms of earthquake forecasting.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
GM, FG, LP, and GF contributed to the collection and interpretation of the geochemical results. NG, ML, and VT contributed to the collection and interpretation of the geophysical results. All authors provided critical feedback and helped shape the research, analysis and manuscript.

FUNDING
Preliminary scientific communications related to present paper were done during the 14th and the 15th International Conference in Gas Geochemistry. This study benefited from funding provided by the Presidenza del Consiglio dei Ministri -Dipartimento della Protezione Civile (DPC), Project S3-2012. This paper does not necessarily represent official DPC opinion and policies. GM was partially supported by the Chinese Academy of Sciences Visiting Professorship for Senior International Scientists (2018VMA0007).