The Seismotectonic Significance of Geofluids in Italy

There is growing interest in how geofluid emissions are released in the atmosphere by the planet’s geodynamic activity, and how much they contribute to the global budget of greenhouse gases. Many workers are addressing this issue with studies conducted at global scale, so as to get the required global-scale answers. The data available at the global scale on geofluids, faults, earthquakes and volcanoes, however, are generally too coarse to provide these answers. We investigate the relationships between geofluid emissions and tectonics at a more detailed scale. Building on over a century of data on geofluid emissions and on an extensive knowledge of the region’s tectonics and seismicity, we focused on Italy, one of the areas of the globe that experience the largest release of natural CO2 and CH4. We systematically overlaid and compared data collected by a number of workers into 13 published countrywide databases concerning geofluid emissions, carbon-bearing deposits, seismogenic faults, historical and instrumentally documented earthquakes, and heat flow observations. Our results indicate that 1) thermal springs and CO2 emissions dominate in areas of mantle upwelling and crustal stretching, but also that 2) some of them occur in the extending inner Apennines, generally along major lithospheric chain-perpendicular lineaments that bound the largest normal faults. Conversely, 3) CH4 emissions and mud volcanoes dominate in areas undergoing active contraction, where no CO2 emissions are observed; in particular, we find 4) that mud volcanoes concentrate where the crests of active anticlines intersect major lithospheric chain-perpendicular lineaments. An overarching conclusion is that, in Italy, the release of geofluids is primarily controlled by deep crustal discontinuities that developed over the course of 5–10 My, and is only mildly affected by ongoing crustal strains. Geofluid emissions bring information on processes that occur primarily in the lower crust, marking the surface projection of generally hidden discontinuities that control the geometry and modes of seismic release. As such they may also provide valuable insight for improving the assessment of seismic hazard in hard-to-investigate seismically active regions, such as Italy.


INTRODUCTION: HOW DO GEOFLUID EMISSIONS RELATE TO TECTONICS AND EARTHQUAKES?
It is well established that pressure gradients existing in the Earth's crust can expel geofluids from deep geological formations toward the atmosphere. The intensity and duration of the expulsion processes depend mostly on crustal permeability (Gleeson and Ingebritsen, 2017) and on the forces that create pressure gradients. Heat from various sources in the core, in the mantle and in the crust generates fundamental tectonic forces at the Earth surface, but additional forces are supplied by topographic relief, diagenetic processes and circulating fluids (Bredehoeft et al., 1990).
Italy lies in the middle of the Mediterranean basin, one of the most active areas of the planet where all these primary forces act simultaneously, generating a large variety of tectonic environments, and hence geofluid assemblages. For many years Earth scientists have attempted to explore and understand the relationships that exist between the location and the characteristics of geofluid emissions and tectonic activity. More specifically, several workers focused on how the occurrence of significant earthquakes is related to the occurrence of anomalous fluids (Irwin and Barnes, 1980;Gold and Soter, 1984;Manga et al., 2012;Chiodini et al., 2020).
The focus was initially on time-domain relationships, and particularly on the detection of possible changes in the chemicophysical characteristics of known thermal springs or gas emissions following a large earthquake (for a worldwide review of the early works see Cicerone, et al., 2009;Wang and Manga, 2010;Martinelli and Dadomo, 2017;Girault et al., 2018;among others). Later, the availability of the first worldwide systematic georeferenced compilations of thermal springs (Waring, 1965), carbon dioxide emissions (Irwin and Barnes, 1980), methane seepages (Etiope, 2015) and other phenomena occurring in specific geological domains, shifted the attention of many workers to the space-domain. Along with the inception of computerized seismicity and active-fault catalogues, starting in the 1990s these data allowed the geofluids community to explore also the spatial relationships between the location of anomalous fluid emissions and the distribution of major faults and large earthquakes (e.g., Kuz'mina and Novopashina, 2018;Tamburello et at., 2018;Ciotoli et al., 2020).
Italy has often been the object of investigations on the origin and characteristics of geofluids in relation with ongoing geodynamic processes. This is due to the combination of its geological diversity and frequent seismicity with its multicentury tradition of field observation and data-keeping. This long tradition relied on the density of historical cities and villages, and hence on the number of potential "casual observers" of geofluid emissions.
Some early investigations are due to Barbier and Fanelli (1976), Minissale (1991) and De Rubeis et al. (1992), all of whom adopted independent research strategies. Barbier and Fanelli (1976: Figure 1) used a point-to-point or point-to-line approach to correlate qualitatively the presence of thermal springs with the location of Italian earthquakes and with "long lineations" derived from ERTS-A satellite imagery. They concluded that no clear spatial correlation could be established between thermal springs and earthquakes epicenters, and proposed that many such springs are aligned with large lineaments, possibly related to faults described generically as "active".
Fifteen years later, Minissale (1991) considered a point-to-area or point-to-line approach to correlate the presence of thermal springs with the most significant CO 2 degassing areas and Quaternary faults. He concluded that a significant part of CO 2 originates in the mantle, and that the main peri-Tyrrhenian fault systems release large amounts of CO 2 which is then dissolved into large carbonate aquifers.
De Rubeis et al. (1992) proposed an area-to-area geostatistical approach to the problem, searching for significant correlations between several geophysical parameters represented on a regular grid laid over Italy. Among other parameters, they correlated the density of Quaternary faults and the historical seismic release with the heat flow, concluding that "...underground temperature and quantity of faults do not present significant areas of correlation...."

TACKLING AN OLD QUESTION THROUGH A NEW APPROACH AND NEWLY ORGANIZED OBSERVATIONS
Exploring the relationships between geofluid emissions and active tectonics has fascinated a large number of geoscientists worldwide over at least the past 50 years, not to mention the important systematization efforts that were made around the end of the XX century. We would like to make a further step forward; to this end we propose a strategy that differs from those used in previous works for three fundamental reasons ( Figure 2): (1) we correlate seismogenic faults, seen as 3D sources projected onto the topographic surface, with point-like occurrences that include thermal springs, CO 2 emissions, gas seeps and mud volcanoes; that is to say, we use an area-to-point approach that allows for a better characterization of the processes leading to the observed emissions; (2) our approach is geographically extensive, as it encompasses the main geodynamic domains of the whole Italian peninsula, but the grain of geofluid emission data is on the order of 1 km; (3) most importantly, our approach is systematic, as it relies on carefully collected computer databases developed at national scale by independent research groups and unavailable to previous workers (Supplementary Table S1). Our paper represents the first occasion for looking at them simultaneously, taking advantage of the enormous recent progress in dealing with georeferenced information.
et al. (2018), for instance, focused on how the permeability of the crust in areas undergoing active extension allows significant degassing of CO 2 -rich fluids, whereas Ciotoli et al. (2020) explored the amount of CH 4 and other hydrocarbons that is released to the atmosphere by faulting in convergent basins. These workers used worldwide compilations of fluids, hydrocarbons and faults. We share their general conclusions on the geodynamic processes that give rise to geofluid emissions, and trust their quantification of such emissions, but the scale of their analyses makes it impossible to explore how the emissions relate to the relevant tectonic structures three-dimensionally. Achieving this goal requires that all analyses be conducted at a more detailed scale, similar to that of the tectonic domains that are responsible for the emissions themselves.
This is exactly what we do in this paper. We hope that our conclusions will provide details that may be of practical use in the investigation of regional seismotectonics. We maintain that geofluid emissions can be used to set firm constraints on the location and extent of major seismogenic sources by highlighting the existence of permanent lithospheric discontinuities; hence they may contribute valuable data for the assessment of seismic hazard in hard-to-investigate active regions, such as Italy.
Our work followed a path that is mirrored in the structure of this paper.
In the Input Data I: Available Databases of Geofluid Emissions, Input Data II: Discriminatory Data and Input Data III: Tectonic and Seismological Observations and Databases sections we describe in detail the data and databases we used, concerning known geofluid emissions, heat flow estimates, the distribution of coal, lignite and peat deposits, and the tectonic and seismological framework: all adopted databases are all freely available (Supplementary Tables S1, S2, S3). We do not discuss the composition of individual emissions, but trust their original classification by the authors of the relevant databases.
In Selecting the Terms of Comparison Section we discuss how we elaborated these databases to suit our specific needs. We focused on emissions occurring in mid-to low-heat-flow areas, that are apparently controlled by the deeper crustal circulation; as such they illuminate processes that occur in the uppermost 10-15 km of the crust, where faults exhibit fully brittle behavior and most large earthquakes nucleate.
The Thermal Springs and CO 2 Emissions Sections illustrate the criteria we used to select our thermal spring and CO 2 data to minimize the "noise" arising from high-heat-flow conditions. Similarly, the CH 4 Emissions and Mud Volcanoes Section describes how we selected our CH 4 and mud-volcano data, discarding biogenic-methane emissions that would bring information only on very shallow-seated processes. Finally, the Earthquake-Induced Geofluid Emissions Section describes how we selected relevant earthquake-induced effects from the CFTI5Med.
In the Correlations and Anticorrelations Section we discuss the correlations, anticorrelations and non-correlations that arise from the superposition of the selected fluid emissions, or earthquake-induced effects, over the available tectonic and seismological information.
Lastly, in the Discussion and Conclusion Sections we combine all data into a global overview, highlighting 1) whether and how the different types of emissions overlap, or are mutually exclusive, and 2) what can be learned from them, both in terms of their geodynamic significance and of the constraints they set on seismogenic processes.

Thermal Springs
Strictly speaking, any spring (or water well: but for simplicity, in the rest of the paper we will refer only to springs) whose average temperature is noticeably above the local mean annual air temperature may be regarded as "thermal" (Waring, 1965). A more accurate definition of "thermal spring" is given in Pentecost et al. (2003): see their Table 1 and their Recommendations  section).
Section S1 of the Supplementary Materials provides an overview of the available knowledge on Italian thermal springs. For this work we decided to refer to the Database of Thermal Springs by Martinelli and Albarello (1997), available for consultation and download as Supplementary Table S2, for two main reasons: 1) to increase the resolution of our analysis, mostly concerned with mildly but positively anomalous springs occurring in areas of moderate to low heat flow, and 2) because for these specific areas the most recent database (Trumpy and Manzella, 2017) reports only a fraction of the datapoints listed in the Database of Thermal Springs, although it does provide many new data for the peri-Tyrrhenian area.
A database of volcanic and non-volcanic CO 2 emissions in Italy was compiled in the framework of project INGV-DPCV5, funded by the Italian Department of Civil Defense (Chiodini et al., 2008). Its aim was to locate, characterize and inventory all Italian gaseous emissions (http://googas.ov.ingv.it). The database has been subsequently updated by Cardellini et al. (2014) in the framework of the MaGa Project (Supplementary Table S1). It currently holds 1,870 natural gas emission measurements from 877 locations distributed worldwide, although the majority of them fall in the Italian territory.

Gas Seeps and Mud Volcanoes
Onshore methane seepages in Italy exhibit a flow rate in the range 0.1-600 t y -1 . Available data on gas seeps in Italy were collected and organized by Martinelli et al. (2012). Most of these emissions are dominated by CH 4 . These investigators merged previous historical catalogues (Camerana et al., 1926;Zuber, 1940) and individual observations from the recent literature (e.g., Etiope et al., 2007). We will refer to gas seeps as "CH4 emissions".
It is widely accepted that CH 4 may be generated by methanogenic micro-organisms, in which case it is referred to as 'biogenic', or by the thermally-driven breakdown of larger molecules, in which case it is referred to as 'thermogenic'. It is also established that biogenic CH 4 generally originates at very shallow depth (<1 km), whereas thermogenic CH 4 generally rises from much deeper rock horizons. Hence, similarly to thermal springs, CH 4 emissions must be carefully selected to remove those dominated by biogenic gas, as they are insensitive to the deepseated tectonic processes that are the object of our investigation.
Mud volcanoes comprise a distinct sub-category of gas seeps. Martinelli and Judd (2004) investigated the Italian mud volcanoes and collected all the available evidence, later turned into a catalogue by Martinelli et al. (2012). A significant part of these data have been used by Burrato et al. (2013) to set up a web site devoted to hydrocarbon seepages in Italy (http:// hydrocarbonseeps.blogspot.com/p/hysed.html).

Input Data II: Discriminatory Data
As we stated earlier, all our geofluid emissions need to be carefully selected to minimize the masking effect of shallow-seated processes and emphasize the role of phenomena occurring at seismogenic depth. This section includes "discriminatory data", i.e., observations that are required only for achieving this goal and will not be used "per se".

Heat Flow Measurements
A number of workers discussed the distribution of heat flow in Italy, including Cataldi et al. (1995), Della Vedova et al. (2001), Pasquale et al. (1997), Pasquale et al. (2010), Pasquale et al. (2014) and Verdoya et al. (2019). We used the surface heat flux density map elaborated and published by Cataldi et al. (1995) and Della Vedova et al. (2001), and made available online as the Italian National Geothermal Database, the largest collection of Italian geothermal data through the 1980s (CNR-IGG, 2012;Trumpy and Manzella, 2017).
The Alps, the Po Plain, the Adriatic and Bradanic foretroughs are characterized by relatively low heat flow values in the range 40-70 mW/m 2 (Figure 3). A major regional anomaly affects large parts of Tuscany and Latium, where the average heat flow is 80 mW/m 2 but it may locally rise up to 200-400 mW/m 2 . Similarly high values are observed in the Campanian volcanoes (Vesuvius, Phlegrean Fields and Ischia Island). According to the literature, springs occur in areas of peninsular Italy where the heat flow is >70 mW/m 2 are representative of hot waters or steam sources originating in high-enthalpy geothermal fields (e.g., Pauselli et al., 2019). Fluid circulation in such areas is strongly affected by magma intrusions at shallow depth (3-5 km), therefore the characteristics of thermal springs cannot be used to make robust inferences about the evolution and the structure of the upper crust.

Distribution of Coal, Lignite and Peat Deposits
Biogenic CH 4 chiefly originates from peats or from lignite/coal layers. In principle, geochemical and isotopic data can distinguish peat-originated methane from thermogenic methane (e.g. Stolper et al., 2015), but, unfortunately, not all Italian methane seepages have been sampled or analyzed over the past few decades. We used a national database of occurrence of fossil carbon recently developed by Martinelli et al. (2015) as a proxy for a first-order separation of biogenic and thermogenic seepages and mud volcanoes.

Data on Active and Seismogenic Faulting
The location, geometry and size of crustal faults, and particularly their length, dip and down-dip width, necessarily affects the nature and distribution of geofluid emissions at the Earth surface. Unfortunately, data on the location and characteristics of active faults carry at least three categories of inherent ambiguities that must be identified and understood prior to attempting the interpretation of a specific geofluid emission: (1) is the fault currently active, and if so, is it also seismogenic?
Or does it slip aseismically? There is no consensus in the current literature as to whether the emission of geofluids requires a network of active faults, or whether even inactive faults may allow geofluids to rise up from shallow-to mid-crustal depth (i.e., from the depth of the brittle-ductile transition at any given site); (2) is the fault positively surface-breaking, i.e. does it reach the topographic surface? Or else, is it blind, meaning that it does not extend to the surface, but its upper tip occurs at a depth ranging from a few to 10-15 km? Notice that blind faults are themselves ambiguous, as they may still reach the surface through secondary, potentially multiple upward-propagating splays; in this case they may give rise to distributed emissions, that are inevitably more difficult to pinpoint; (3) what is the hierarchy of any given fault with respect to all the faults that lie in the same area? In any seismotectonic application it is crucial to assess whether a certain fault extends to seismogenic depth, or whether it is only a subsidiary to some larger fault beneath it, which is hence the main player in the area. Fluid emissions associated with a secondary fault are likely to be representative of very shallow-seated processes only, unless that fault belongs to a larger and deeper fault system.
Choosing an appropriate fault dataset is therefore crucial in the correct interpretation of geofluid emission and in the Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 579390 investigation of preferred pathways. Italy has two main faults datasets, both of which were conceived around the end of the 1990s: ITHACA (ITaly HAzard from CApable faults, Supplementary Figure S2, left), developed by ISPRA (Michetti et al., 2000;ITHACA Working Group, 2019), and DISS (Supplementary Figure S2, right), developed by INGV Pantosti, 2000, Valensise andPantosti, 2001;Basili et al., 2008;DISS Working Group, 2018). Their characteristics and contents are described in Section S2 of the Supplementary Materials. In addition to potential earthquake sources, the earlier versions of DISS (v. 2.0: http://diss.rm.ingv.it/diss/index.php/ component/chronoforms5/?chronoform DISS200) listed over 140 "transverse lineaments", a category of large, chainperpendicular structural features inherited from previous tectonic regimes. It must be recalled that the Cenozoic compressional architecture of the Italian peninsula largely reflects the inherited Mesozoic palaeogeographic setting of its foreland, i.e. the Adriatic plate. There is ample evidence that inherited structures control the foreland monocline panels that lie below the basal detachment of the accretionary prisms and exhibit variable dip (Mariotti and Doglioni, 2000). These are commonly interpreted as evidence for the existence of slab tears (e.g., Rosenbaum et al., 2008;Rosenbaum and Piana Agostinetti, 2015), whose crustal expression could be strike-slip or transpressional transfer zones (Vannoli et al., 2015): pseudolinear features resulting from the rectification of potentially complex elements but usually perpendicular to the thrust fronts, and hence "transverse" because they form an angle of 45-90°with respect to the dominant local structural trend. The role of such lineaments in the present geodynamic setting of Italy and their potential for substantial earthquakes are fully open to scientific debate, although there is growing evidence that they do generate both frequent background seismicity and relatively large earthquakes.
The Individual Seismogenic Sources (hereinafter ISSs) and Composite Seismogenic Sources (hereinafter CSSs) of the DISS database (introduced in Section S2) include also a number of seismogenic sources that, similarly to the transverse lineaments, form an angle with the main structural trends, but whose seismogenic capability is testified by several large preinstrumental earthquakes and confirmed by a handful of instrumentally investigated earthquakes. Most of these sources are believed to indicate the reactivation of older and deeper strikeslip faults that crisscross the foreland areas of both the Alps and the Apennines, although the principal ones concentrate between northern Apulia and eastern Sicily (Di Bucci et al., 2010). Two examples of the seismogenic potential of such deeper faults are the May 5, 1990, Potenza earthquake (M w 5.8) and the 31 October-1 November 2002, San Giuliano di Puglia earthquakes (both of M w 5.7), all of which were generated by dextral motion on E-W striking strike-slip faults at 15-25 km depth. We postulate that these deeper faults may act as pathways for the rise of thermally anomalous fluids and CO 2 , even in areas of low heat flow. Their potential role in controlling the distribution of geofluid emissions will be discussed in Correlations and Anticorrelations Section.
Inevitably, the use of these two databases in our investigation is subject to uncertainties and limitations. For instance, lacking a 3D characterization and hierarchization, only some of the faults listed in ITHACA may reach a significant depth and hence be considered as viable pathways for upper crustal geofluids. Conversely, DISS lists only primary faults, i.e. faults extending to seismogenic depth, and describes them in 3D, but the accuracy FIGURE 3 | Heat flow map of Italy from CNR-IGG (see text for details). Shaded areas indicate active tectonic domains (same as in Figure 2). Left -Complete model, showing high/medium heat flow areas in the peri-Tyrrhenian belt (contours in red/orange). Contour spacing is necessarily variable as it reflects the original elaborations. Right -Outline of the 70 mW/m 2 isoline (shown in dark grey), used in this work for separating areas were thermally anomalous fluids and CO 2 emissions are mostly related to high heat flow, to the west, from areas where they may represent the outcome of fluid circulation in upper crustal rocks, to the east.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 579390 of their location may not always be sufficient to explore their relationships with fluid emissions. In Nevertheless, in the following discussion we will rely primarily on the faults described in the DISS database, as they are more likely to be relevant for the scopes of this work.

Earthquake Catalogues
We rely on two major, regularly updated earthquake compilations available for Italy and its surroundings: the Catalogue of Strong Italian Earthquakes (Guidoboni et al., 2018(Guidoboni et al., , 2019, a large compilation that summarizes the knowledge acquired through nearly 4 decades of modern research in historical seismology, and the Database Macrosismico Italiano, or DBMI , whose current version (2.0) provides macroseismic intensity datasets relating to Italian earthquakes in the time window 1,000-2017.
In some instances we plotted also instrumental seismicity recorded by the INGV national network and by other local networks over the past 35+ years (from ISIDe-Italian Seismological Instrumental and Parametric Database: ISIDE Working Group, 2007).
The characteristics and contents of these compilations are described in Section S3 of the Supplementary Materials.

SELECTING THE TERMS OF COMPARISON
The tectonic and seismological significance of the databases upon which our analysis is based may differ with the geographical area concerned, the local predominance of subcrustal geodynamic processes, and fracturing that may have been inherited from older tectonic regimes. In this section we discuss how the geofluids observations have been selected from their respective databases to remove potential noise and emphasize their relationships with deep-seated seismotectonic processes. But first we develop a simplified scheme of the dominant strain regime acting in the mid-upper crust over the Italian peninsula and its surrounding areas and seas.
In conjunction with seismicity (Guidoboni et al., 2018;Guidoboni et al., 2019;Rovida et al., 2019;Rovida et al., 2020) and tectonic stress data (Montone and Mariucci, 2016;Mariucci and Montone, 2020), the availability of orderly fault databases allows the most active portions of the Italian peninsula to be separated into four distinct tectonic domains. The resulting model is shown in Figure 2 and is replicated in Figures 3-8 to serve as context.
We first subdivided the study area into a dominantly extensional and a dominantly contractional domains. The former was further subdivided into a peri-Tyrrhenian belt of Plio-Quaternary rifting and mantle upwelling (e.g. Peccerillo, 2005: shown in light blue in Figure 2), and an elongated region of post-orogenic extension running along the axis of the Apennines (DISS Working Group, 2018: shown in blue), where large shallow normal faults cause most of the largest earthquakes of the Italian peninsula. This latter region has also undergone continuous regional uplift over the past 1-2 My; a very effective process in causing erosional denudation and promoting the visibility of the geofluid emissions. The purely contractional domain extends along the frontal portions of the Alps and of the Apennines (frontal with respect to the direction of orogenic transport), i.e. in the Po Plain, in the peri-Adriatic and peri-Ionian regions, and in much of the corresponding offshore area (DISS Working Group, 2018: shown in red).
The Italian peninsula is also surrounded by foreland areas that extend along the Adriatic microplate, i.e. from the central and eastern Po Plain to the Adriatic and Ionian offshore, in southern Sicily and in the Sicily Channel. These foreland areas are subjected to NW-SE compression exerted by the northwestward motion of the African plate: the resulting contraction is primarily accommodated by deeper strike-slip faults that are known to exist in the central and southern Apennines (already mentions in Data on Active and Seismogenic Faulting Section), where they strike nearly E-W, and in southern Sicily, where they strike NNE-SSW (e.g., Di Bucci et al., 2010: all shown in green in Figure 2). Notice that the "post orogenic extension" and "deeper strike-slip" domains overlap in limited areas of the southern Apennines.
In the following, for all databases, we will refer to the selected and removed data as the "selected dataset" and "removed dataset", respectively. In our figures each emission type is shown with a specific color, but all removed data are shown in grey, whereas red is used to indicate specific emissions that are referenced in the discussion.

Thermal Springs
The average annual temperature in mountain and plain areas of Italy is in the range 5-20°C, thus for this work we decided to focus only on springs whose mean temperature falls in the range 20-150°C (see Figure 4). The most common temperature range for the springs reported in the Database of Thermal Springs (see Thermal Springs Section) is 25-40°C. Springs in this temperature range are widespread all throughout the Italian peninsula, providing evidence for ongoing tectonic activity: as recently discussed in a worldwide overview by Bense et al. (2013), a high density of thermal springs is always associated with ongoing tectonic activity, whereas thermal springs tend to be very sparse or absent in tectonically quiescent areas.
The recent tectonic history of Italy, and the goals of our study, required that all thermal springs falling in the geothermal or volcanic areas along the Italian peninsula be removed prior to the beginning of our analysis. As mentioned earlier, springs associated with high heat flow carry information only on shallow-seated processes, hence they are useless or even misleading if one is to explore processes driven by forces acting on the full thickness of the seismogenic layer.
It is known that the "excess heat" observed in high heat flow areas (e.g., Tuscany) is large enough to mask more localized thermal contrasts between warm and cold springs. For instance, Celati et al. (1990) demonstrated that when the heat flow reaches 100-200 mW/m 2 , no faulting due to tectonic activity-that is to say, no contribution from deep-seated fluid circulation-is needed to justify the occurrence of thermally anomalous springs. Pauselli et al. (2019) stated that "Five crustal geotherms (extending to the Moho) and the corresponding rheological Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 579390 profiles confirm that the 70 mW/m 2 isoline corresponds to a major tectonic boundary, across which the thermal, structural, and seismic properties of the lithosphere go through a significant change". In the framework of a large worldwide compilation of heat flow data, Lucazeau (2019) found that 70 mW/m 2 is a global average. We hence postulate that warm or hot springs falling in areas where the heat flow is <70 mW/m 2 are likely representative of thermally anomalous fluids reaching the ground surface by means of geological structures due to sustained tectonic activity at scales of ∼100,000 years. In summary, prior to analyzing the data we removed all thermal springs falling in areas where the heat flow is >70±5 mW/m 2 (to account for uncertainties in the heat flow estimates, particularly in northern Tuscany and Campania: see Figure 3). Table 1 shows the distribution of data into the four seismotectonic domains described earlier, before and after the selection process.

CO 2 Emissions
CO 2 emissions listed by Cardellini et al. (2014) consist of mofete characterized by flow rates in the range 1-2,000 t d −1 , and can be considered unaffected by non-geological contributions. Nevertheless, and similarly to the case of the thermal springs, our interpretations require that CO 2 emissions originating in volcanic and geothermal areas be removed (see Figure 5). CO 2 emissions originate at different depths and by at least two different primary processes: 1) the anatexis of carbonate rocks at relatively shallow depth (3-5 km), which dominates in areas of mantle upwelling due to high heat flow, and 2) degassing from the mantle lithosphere, which is virtually ubiquitous in geodynamically active areas such as the Italian peninsula. According to Chiodini et al. (2004), the discharge of CO 2 in Italy reflects a combination of these two sources, with the upper crustal component dominating in high-heat-flow areas and the mantle component prevailing elsewhere. The "shallow" and FIGURE 4 | Distribution of thermal springs (see Thermal Springs, Thermal Springs and Thermal Springs: Potential Evidence for Lithospheric Discontinuities). Tectonic domains, seismogenic sources and transverse lineaments are shown as in Figure 2. Removed datapoints are shown in grey. Selected datapoints are shown in red if they fall within 5 km of the edge of the surface projection of an Individual Seismogenic Source, otherwise in white. Red stars indicate additional critical datapoints discussed in the text. Two rectangles indicate the area discussed in Test Case 1: The Garfagnana Region, Northern Apennines and Test Case 2: The Olevano-Antrodoco-Sibillini Line, Central Apennines Sections. The pie chart shows the distribution of thermal springs in the different active tectonic domains (see Table 1). Notice that 42% of the springs fall in the area of rifting and mantle upwelling.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 579390 "deep" CO 2 emissions have a distinct δ 13 C signature, but analytical results are available only for a minority of cases. We used the heat flow as a proxy to discriminate between "possibly shallow" and "positively deep" CO 2 ( Figure 5), using the same threshold adopted for discriminating "shallow" and "deep" thermal springs (70 mW/m 2 ).

CH 4 Emissions and Mud Volcanoes
As we discussed in Distribution of Coal, Lignite and Peat Deposits Section, CH 4 emissions must be selected to remove those dominated by biogenic gas, whose shallow origin is incompatible with the goals of this paper. To this end we removed all CH 4 seepages falling within 5 km of any occurrence of fossil carbon as shown in the database compiled by Martinelli et al. (2015: see Distribution of Coal, Lignite and Peat Deposits Section and Figure 6). In principle, all seepages left should be of thermogenic, i.e. of deeper origin.
In consideration of their genetic similarities, in parallel with the CH 4 seepage dataset we also examined all mud volcanoes listed in Martinelli et al. (2012). These data were also interpreted after removal of occurrences that fall close to a known fossil carbon reservoir (Figure 7).

Earthquake-Induced Geofluid Emissions
Earthquake-induced effects on the environment were one of the main innovations of CFTI5Med (see Earthquake Catalogues Section). Its authors systematically analyzed, geo-referenced observations and filed testimonies describing 2,337 such effects. Notwithstanding the difficulties involved in translating historical testimonies into corresponding natural phenomena as we know them today, these effects were classified into 32 categories belonging to five macro-groups. One of these groups is named "Gas emissions/Others" and contains two subcategories that are relevant to the scopes of this paper: Selected datapoints are shown in red if they are individually discussed in the text, otherwise in white. The areas marked with A, B and C are discussed in Test Case 1: Central-Northern Apennines, Test Case 2: Central-Southern Apennines and Test Case 3: Southern Apennines Sections, respectively. The pie chart shows the distribution of CO 2 emissions in the different active tectonic domains (see Table 1). Notice that as many as 61% of the CO 2 emissions fall in the area of rifting and mantle upwelling.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 579390 1) Gas emissions and mud-cone eruptions, containing 79 testimonies, the oldest from the great December 5, 1456 southern Italy earthquake sequence; and 2) H 2 S emissions, containing 60 testimonies, the oldest from the October 3, 1624 earthquake in southeastern Sicily. These testimonies are plotted in Figure 8 and will be discussed in subsequent sections.

CORRELATIONS AND ANTICORRELATIONS
In this section we discuss what can be learned about the seismotectonics of Italy by overlaying all geofluid observations selected in Selecting the Terms of Comparison Section over the countrywide distribution of large earthquakes, seismogenic faults and other major tectonic lineaments. In the process, we plan to learn also how the distribution of thermal springs, CO 2 emissions and CH 4 emissions relates to the recent geodynamic evolution of the Italian peninsula.
Thermal Springs: Potential Evidence for Lithospheric Discontinuities Figure 4 shows the distribution of our "selected" and "removed" datasets of thermal springs (see Thermal Springs Section). A first notable feature is that they tend to occur rather regularly along the axis of the Apennines. To investigate this circumstance more quantitatively we calculated a variogram, i.e. a description of the spatial continuity of the data, using the classical methodology described and supplied by Pannatier (1996). Figure 9 shows variograms calculated for any azimuth and for a 135°N trend; this corresponds to the average strike of the Apennines, from their northern boundary, near the Tuscany-Liguria border (∼44.0°N), to the Pollino Massif (∼40.0°N), near the Basilicata-Calabria border, over a distance of about 700 km. Comparison between the two trends shows that the most significant peak is found along the Apennines trend at a distance of about 30 km. In other words, the dominant feature of the spatial distribution of all springs is a ∼30 km spacing in the NW-SE direction. This distance is slightly larger than the average length of the largest  Table 1). Notice that 54% of these emissions fall in areas of active contraction.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 579390 Individual Seismogenic Sources (see Data on Active and Seismogenic Faulting Section) that are encountered along the crest of the Apennines. Current empirical relationships show that in a typical normal-faulting environment, such as that dominating the current tectonics of the Apennines, a 30 kmlong fault typically generates a M w 6.8 earthquake: this estimate matches the range of the largest observed M w of Apennines earthquakes (6.5-7.0). The spacing between thermal springs tends to decrease in geothermal areas (e.g., McNitt, 1978;Minissale, 2018), and to increase in cooler, tectonically stable areas. These features were qualitatively described over 50 years ago by Waring (1965), who compiled a rather complete worldwide catalogue of thermally anomalous springs.
We then overlaid the "selected dataset" onto the seismogenic sources (ISSs and CSSs) and the "transverse lineaments" introduced in Data on Active and Seismogenic Faulting Section. This procedure reveals a rather systematic anticorrelation between earthquake sources and the distribution of thermally anomalous springs all along the Apennines. A closer inspection shows that (1) thermal springs appear to be anticorrelated with the surface projection of major seismogenic sources, implying that the relevant crustal volumes are less permeable than adjacent areas. This phenomenon has also been inferred by Micklethwaite et al. (2010) from the distribution of ore deposits, which occur mainly along transverse features located at fault tips, rather than along large, welldeveloped faults; (2) conversely, many thermal springs cluster at the edges of the best defined sources, as suggested over 2 decades ago by Curewitz and Karson (1997). More precisely, in many instances where well-identified sources define a clear segmentation pattern, the springs seem to occur along tectonic lineaments that mark the transition from one segment to the other. The evidence from the DISS database suggests that these discontinuities act as segment  Table 1). Notice that as many as 85% of these peculiar phenomena fall in areas of active contraction.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 579390 boundaries for earthquake strain accumulation, whereas the evidence from our selected thermal springs is that they also allow deep fluid circulation to reach the surface. Test Case 1: The Garfagnana Region, Northern Apennines Section discusses a significant example from the Garfagnana region (northern Tuscany);and (3) some thermal springs (and CO 2 emissions) line up along major lithospheric lines of weakness that cross the Apennines; these lines may hold significant seismogenic potential, particularly in the southern Apennines, or be only mildly seismogenic as a result of passive reactivation by regional tectonic forces. Test Case 2: The Olevano-Antrodoco-Sibillini Line, Central Apennines Section discusses the example of the Olevano-Antrodoco-Sibillini Line.
The postulated correlation between thermally anomalous springs and transverse lineaments is supported by recent experience during the search for exploitation of natural geofluids (hydrocarbons, steam etc.). Important physical constraints were reviewed by Hickman et al. (1995), Caine et al. (1996), Holness (1997), Manning and Ingebritsen (1999), Sibson (2000), Kennedy andvan Soest (2007), Faulkner et al. (2010), Ingebritsen and Manning (2010), Micklethwaite et al. (2010), Bense et al. (2013), andJang et al. (2018), among several others. These investigators have shown that fluid flow through fractured rock masses depends on their permeability, and specifically 1) on the hydraulic conductivity of fractures, and 2) on the degree to which fractures form continuous pathways or barriers through the rock volume.
In the next two subsections we discuss at a higher spatial resolution the interplay of tectonics and thermal springs in two selected test areas of the northern and central Apennines, respectively.

Test Case 1: The Garfagnana Region, Northern Apennines
The area referred to as Garfagnana ( Figure 10) coincides with the upper and middle reach of the Serchio River valley at the northwestern edge of Tuscany. Garfagnana is located between the axis of the Apennines and the smaller Apuan Alps, which separate it from the Tyrrhenian coastline. According to the DISS database, the Upper Serchio Valley was the locus of the September 7, 1920 catastrophic earthquake (M w 6.5). The The X axis reports the distance between pairs of springs, whereas the Y axis shows the value of the variogram function. The most significant peak is found along the function calculated along the Apennines trend at slightly less than 30 km distance.
FIGURE 10 | Seismogenic faults, transverse lineaments and thermal springs in the Garfagnana area, northern Tuscany (see Figure 4 for location). The black boxes are the surface projections of Individual Seismogenic Sources from the DISS database: source ITIS050 is assumed to have caused the September 9, 1920, M w 6.5 earthquake, while source ITIS051 is assumed to be quiescent and capable of generating a M w 6.1 earthquake. The transverse lineaments are shown with yellow dashed lines and identified by a number (see text). The springs are shown by diamonds and are labeled with their ID from the Database of Thermal Springs. The 3D spiral symbol locates a H 2 S emission that was activated following the 1920 earthquake (Casoni site, from CFTI5Med; see Earthquake-Induced Geofluid Emissions Section). Notice that the most recent earthquakes fall near the edges of the largest seismogenic normal faults, that their faulting mechanism is consistently strike-slip.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 579390 13 same database shows that the Lower Serchio Valley hosts a seismogenic source similar to that responsible for the 1920 earthquake, but which has not been active historically over its entire length of 15-18 km. These two Individual Seismogenic Sources are labeled in Figure 10 (ITIS050-Garfagnana North and ITIS051-Garfagnana South, respectively).
The entire region is crisscrossed by a number of pre-existing tectonic structures referred to as "tectonic lineaments" (see Data on Active and Seismogenic Faulting Section), which seem to control the location and length of the two mentioned seismogenic sources, effectively acting as segment boundaries.
Several significant thermal springs and earthquakes fall along all three of these lineaments, again suggesting that they are not shallow-rooted features generated during the latest compressional tectonic phase, but rather mark the surface expression of important lithospheric discontinuities, as highlighted in the geological literature (e.g., Fazzuoli et al., 1985;Puccinelli et al., 2015). In detail (each spring is identified with its name, ID and reported temperature): Sarzana-Equi Terme Line (#1)-Spring: Equi Terme, ID 267 (24°C). This lineament bounds the Apuan Alps to the north and coincides with a seismogenic relay ramp discussed in Stramondo et al. (2014), following the June 21, 2013, M w 5.1, Lunigiana earthquake. This event was caused by oblique slip on a dominantly normal fault striking perpendicular to the main Apennines trend, and may represent a smaller repeat of the April 11, 1837, M w 5.9, Equi Terme earthquake. The lineament is supposed to have hosted the October 10, 1995, M w 5.0 earthquake, that is assigned a pure strike-slip mechanism (Frepoli and Amato, 1997).
Massa-Monte Cervarola Line (#2)-Springs: Turrite Secca, ID 347 (34°C) and Pieve Fosciana, ID 811 (39°C), both reported as falling on recent surface faults in the geological literature (Puccinelli et al., 2015). This lineament marks the Monte Perpoli threshold, a physiographic and structural discontinuity the separates the Garfagnana basin into two halves. It was the locus of the January 25, 2013, M w 4.7, strike-slip earthquake. According to CFTI5Med, the September 7, 1920 event triggered a H 2 S emission near Casoni, about 15 km east of the causative fault.  Figure 4 for location). The black boxes are the surface projections of Individual Seismogenic Sources from the DISS database: they are assumed to have caused-from north to South-the September 26, 1997 (M w 6.0), the October 14, 1997 (M w 5.6), the August 24, 2016, M w 6.2, Amatrice (sources from Valensise et al., 2016), the February 2, 1703, M w 6.7, and the April 6, 2009, M w 6.3, L'Aquila earthquakes. Transverse lineaments are shown by yellow dashed lines, and the OASL is highlighted in red. Thermal springs are shown by blue diamonds along with their ID from the Database of Thermal Springs. CO 2 emissions are shown by green diamonds and labeled with their name from the database available at http://googas.ov.ingv.it. The 3D spiral symbols locate gas emissions, mud cone eruptions and H 2 S emissions that were activated following the 1703  Viareggio-Val di Lima-Bologna Line (#3)-Springs: Bagni di Lucca cluster, ID 111-129, formed by 19 independent springs known and exploited for at least two millennia (34-55°C). This lineament transects the southern Apuan Alps, mirroring Lineament #1. It was the locus of two pure strike-slip earthquakes on June 7, 1980, M w 4.6 (Eva and Solarino, 1992) and September 7, 2014, M w 4.1 (http://terremoti.ingv.it/tdmt; Scognamiglio et al., 2006), and of the October 27, 1914, M w 5.6 earthquake, which has recently been assigned a focal depth of 24 km (Sbarra et al., 2019).
Additional springs are reported in the official geological map of Italy (Puccinelli et al., 2015) as falling on the surface projection of the Garfagnana South source near its intersection with the Massa-Monte Cervarola Line. The most important is the Gallicano spring, ID 487 (23°C).

Test Case 2: The Olevano-Antrodoco-Sibillini Line, Central Apennines
The NNE-SSW Olevano-Antrodoco-Sibillini Line (OASL, formerly known as Ancona-Anzio Line: Parotto and Praturlon, 1975;Castellarin et al., 1978, Castellarin et al., 1982Patacca et al., 1990;Calamita et al., 2012;Di Domenica et al., 2014, with references) marks the contact between units overlain onto the Latium-Abruzzi carbonate platforms and units belonging to the Umbro-Marchean Apennines (see Figure 11). It is interpreted as the surface trace of a regional paleogeographic feature that may also have acted as a lateral ramp or transpressional fault during the buildup of the Apennines, and its deeper portion may coincide with a large inherited Mesozoic fault.
The OASL is easily recognized only in the portion of the Apennines chain that is currently undergoing uplift and extension (in dark blue in Figures 4-8): nevertheless, according to the literature (e.g., Castellarin et al., 1978), and as suggested by thermal springs and CO 2 emissions (see Figure 11), and even by an align-ment of travertine deposits (Minissale, 2004: see his Figure 4), the OASL may well extend beyond the boundaries of this domain, beneath the volcanic deposits of the peri-Tyrrhenian area and the young terrigenous deposits of the peri-Adriatic area.
Furthermore, the OASL is the locus of several minor earthquakes, all seemingly shallow-seated (focal depth <10 km), that occur spatially and temporally clustered along its entire length of about 130 km. For instance, 13 clearly felt earthquakes are reported over a period of 33 years (Apr 7, 1930, M w 4.5;Nov 9, 1930, M w Figure 11).
The OASL acted as a major segment boundary during the mentioned 14 January-2 February 1703 earthquake sequence, when two M w 6.9 and 6.7 shocks occurred just to the NW and to the SE of the lineament and a M w 6.0 occurred on 16 January very close to the OASL itself (from CFTI5Med). The OASL is suspected to have played a similar role during other earthquake sequences, such as the August-October 2016 events in the Amatrice-Norcia-Visso area Bonini et al., 2019). The OASL falls between the two fault patches that experienced the largest slip during the August 24, 2016, M w 6.0, Amatrice earthquake, although this is commonly considered a single event (shown by two black boxes and red labels in Figure 11).
In summary, a distinct alignment of thermal springs and CO 2 emissions in a low-heat-flow inner portion of the peninsula highlights the presence of the OASL, a major tectonic lineament that played-and still plays-an important role 1) in the geological evolution of the central Apennines, and 2) in controlling the segmentation of seismogenic faults in central Italy, while 3) releasing its own low-level but continuous seismicity.

CO 2 Emissions: Combined Evidence for Carbonate Dissolution and Mantle Degassing
Natural CO 2 emissions are characteristic of areas undergoing extension. A recent paper by Tamburello et al. (2018) explored the relationships between CO 2 degassing and active faulting on a global scale, concluding that this process is greatly favored by the presence of surface-breaking normal faults. This is one of the first papers that focused explicitly on the relationships between CO 2 emissions and active faulting: some of its conclusions, however, are disproved by our observations on central Italy, one of the areas of the globe where this process has been best investigated. Figures 2, 5 show that Italian CO 2 emissions occur dominantly in the peri-Tyrrhenian domain (light blue band in both figures; see pie chart in Figure 5 and Table 1). Very few are seen along the main axis of the central and southern Apennines, and none are reported for the northern Apennines, except for the Pieve di Santo Stefano-Caprese Michelangelo area (discussed in Test Case 1: The Garfagnana Region, Northern Apennines Section), and for some CO 2 associated with gas seeps, which were demonstrated by Capozzi and Picotti (2010) to be not of deep origin. Furthermore, no CO 2 is generally reported in the composition of fluids rising from several wells drilled along the axis of the Apennines, mostly in the early phases of the hydrocarbon exploitation in Italy (see VIDEPI Project: https:// www.videpi.com/videpi/videpi.asp). Finally, CO 2 emissions are virtually absent in areas where large reverse and thrust faults accommodate ongoing contraction (red bands in Figure 5). Most importantly, CO 2 emissions are largely anticorrelated with the occurrence of large seismogenic normal faults (blue band in Figure 5), and in particular with the Individual Seismogenic Sources. In fact, the Italian case suggests that the control by extensional tectonics invoked by Tamburello et al. (2018) is restricted to areas of ongoing crustal rifting and mantle upwelling, such as the peri-Tyrrhenian volcanic belt of Tuscany, Latium and northern Campania. Therefore, their conclusion that there exists a "...positive spatial correlation between gas discharges and extensional tectonic regimes..." must be partially reconsidered, at least for Italy and for areas that share a similar setting worldwide. As discussed earlier, there is a useful distinction between extension generated by crustal thinning, that dominates in the peri-Tyrrhenian domain, and post-orogenic extension, the engine of nearly all seismicity in central Italy (and of most large earthquakes countrywide: e.g., Basili et al., 2008;Meletti et al., 2008;Carafa et al., 2015, all with references: see Figure 2, right). Similarly, the "...concordance between CO 2 discharge, seismic activity, and major faults..." is a primary characteristic of the rift area the authors use to support this statement (the East African Rift), but is not observed in Italy, as shown by our data.
Our observations have two main implications. The first is that generic tectonic extension alone is not sufficient to promote a significant discharge of CO 2 . The second is that, at least in Italy, CO 2 discharge is abundant almost exclusively where the heat flow is rather high (>100 mW/m 2 ), implying that most emissions may originate from carbon dissolution at relatively shallow depth. The ascent of CO 2 may be favored by a network of small-scale normal faults and extensional joints, but in high-heat-flow areas the rheology of the relevant crustal volumes prevents the generation of larger faults, and hence of significant earthquakes. In fact, the occurrence of significant earthquakes is largely anticorrelated with the existence of abundant CO 2 emissions, as clearly seen over the entire peri-Tyrrhenian extensional domain and specifically in central and western Tuscany (see Figures 2, 5). Tamburello et al. (2018) also contend that CO 2 degassing is favored by the subvertical attitude of normal faults, whereas in contractional domains the rise of mantle fluids would be prevented by the low-dip of thrust faults and by the stacking of multiple thrust sheets. The evidence available for the Italian peninsula, however, does not support this statement. For a combination of reasons that we will not discuss here [the reader may refer to Doglioni et al. (1999), andTiberti et al. (2017), for a discussion of this topic], the attitude of normal faults generated by post-orogenic extension is generally rather low. One of the principal extensional structures of central Italy, commonly referred to as the Etrurian Fault System (e.g., Boncio et al., 2000, among several others), dips consistently less than 45°, the opposite of what would be expected for a typical Andersonian normal fault, and often less than 30°. In fact, many well investigated Italian seismogenic normal faults are positively assigned a low dip, including the causative faults of the September 26, 1997, April 6, 2009 (see Table 5 in Vannoli et al., 2012), August 24, 2016  and October 30, 2016 (see Tables 1, 2 in Bonini et al., 2019) earthquakes.
Nevertheless, there are significant departures from the above rules, and substantial CO 2 degassing may take place even in the innermost portions of the Apennines. In Data on Active and Seismogenic Faulting Section we discussed the existence and role of the "transverse lineaments", most of which straddle the crest of the Apennines. In Test Case 1: The Garfagnana Region, Northern Apennines and Test Case 2: The Olevano-Antrodoco-Sibillini Line, Central Apennines Sections we showed examples of how these discontinuities affect both the location of geofluid emissions, and particularly of thermal springs, and the geometry and extent of the seismogenic sources. In the next three subsections we discuss at a higher spatial resolution the interplay of tectonics, CO 2 emissions and thermal springs in three selected test areas of the central-northern, central-southern and southern Apennines, respectively.

Test Case 1: Central-Northern Apennines
A spectacular case of CO 2 degassing is found in the area of Pieve di Santo Stefano-Caprese Michelangelo (shown with A in Figure 5), near the boundary of Tuscany, Umbria, Marche and Romagna. Degassing is so abundant that the CO 2 has long been exploited by the bottled beverage industry (e.g., Chiodini et al., 2004). Here a 20 km-wide, NE-SW trending discontinuity appears to break and displace both the main extensional trend running along the crest of the Apennines, and the trend of compressional structures running along their eastern (Adriatic) flank (see Figure 2, right). It is generally interpreted as a long-lived lithospheric tectonic lineament cross-cutting the chain from the Tyrrhenian to the Adriatic Sea (e.g., Nirta et al., 2007 and reference therein). As such, it experienced different kinematics through time, accommodating the differential motion between different portions of the Apennines during the Miocene rotation of the chain and the subsequent opening of the Tyrrhenian Sea (Nirta et al., 2007). Bonini (2009) investigated the structural control of the abundant release of CO 2 -rich fluids at Pieve Santo Stefano. Based on a detailed fault-pattern analysis, he recognized that CO 2 emissions locate along steep NE-trending faults associated with the Val Marecchia Line, perhaps the most important of all the transverse lineaments crossing this portion of peninsular Italy. According to Nirta et al. (2007) and Bonini (2009), the Val Marecchia Line corresponds to a set of overstepping, NE-SW major transverse lineaments described in the geological literature for over 50 years (starting with Ghelardoni, 1965), and also with a distinct promontory of the heat flow pattern (Pauselli et al., 2019: see Figure 3). Surprisingly, this corridor marks also a substantial gap in seismic moment release along the crest of the Apennines, in terms of both major earthquakes and background seismicity; Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 579390 perhaps the most significant gap of the entire northern and central Apennines.

Test Case 2: Central-Southern Apennines
A roughly E-W, 55 km-long alignment of six CO 2 -rich emissions is found west of the Apennines regional divide at an average latitude of 41.65°N, near the Lazio-Molise administrative boundary (from west to east: Fontana Liri, Lago Fibreno, Polla Monticchio, Casa Ianni, Capo Volturno, and Isernia-Castel Romano; shown with B in Figure 5). It falls along the westward extension of a well-known, ∼180 km-long, ∼20 kmwide, E-W-trending major shear zone that affects (from west to east) portions of the Adriatic crust beneath the external Apennines thrust belt, the Apulian foreland and the Molise-Mattinata-Gondola shear zone (see Di Bucci et al., 2010;Kastelic et al., 2013 and references therein). This large discontinuity is the locus of recent, instrumentally investigated earthquakes, such as the Oct 31-Nov 1 2002 Molise events (both with M w 5.7), and of a number of large historical earthquakes, such as the July 30, 1627, M w 6.7 event.
The DISS database represents this important seismogenic area with a large CSS and eight ISSs, one of which is assumed to be the source of one of the main shocks of the M w 7.2, December 5, 1456 earthquake, the largest-or second largest-Italian earthquake of all times (Fracassi and Valensise, 2007). Available focal mechanisms consistently show right-lateral strike-slip over the entire structure, which is commonly interpreted as an inverted, left-lateral Mesozoic shear zone (Di Bucci et al., 2010).
The identified alignment of CO 2 emissions appears to continue this large shear zone westward, suggesting its possible extension to the west beneath the Apennines foothills of southern Latium, and may shed light on the origin of damaging historical earthquakes that occurred further to the west, between Sora and Frosinone.

Test Case 3: Southern Apennines
A further E-W alignment of at least 11 CO 2 emissions and 24 thermal springs is seen west of the Apennines drainage divide around the latitude of 40.70°N. From west to east the CO 2 emissions are found at San Benedetto, near Polla and Contursi (cluster of nine emissions), and at Terme di San Cataldo (shown with C in Figure 5 Similarly, to the previous case, this discontinuity follows a well-known, though much shorter, E-W-trending, deep shear zone that has been the locus of instrumentally-investigated strikeslip earthquakes, such as the May 5, 1990, Potenza event (M w 5.8).
The same shear zone could be the source of the July 31, 1561, M w 6.3, and August 19, 1561, M w 6.7 earthquakes, the two largest shocks of a sequence that affected the Basilicata-Campania border in historical times.

CH 4 Emissions and Mud Volcanoes: The Signature of Contractional Domains
The occurrence of CH 4 emissions and mud volcanoes exhibits a spatial pattern that is complementary to that of the thermal springs and CO 2 emissions. Ciotoli et al. (2020) recently explored the modes and rates of CH 4 seepage, concluding that this process is chiefly associated with convergent, hydrocarbon-bearing sedimentary basins. They added that the details of the interplay between seepage and tectonics are generally unknown except for a few local studies. We maintain that in this respect Italy makes a perfect case study, due to the combination of good-quality data on tectonics and active faulting, an extensive dataset of CH 4 seepages, data from carbon-bearing deposits, and data from several hundred wells drilled for hydrocarbon exploration starting in the 1950s (see the Gas Seeps and Mud Volcanoes and CH 4 Emissions and Mud Volcanoes Sections).

CH 4 Emissions and Mud Volcanoes: Distribution
A simple inspection of our dataset (Figures 6, 7, Table 1) shows that the vast majority of mud volcanoes (85%) occur in regions currently undergoing compression and shortening, whereas gas seeps are more sparse: 54% fall in areas currently undergoing contraction, but 20% of them are 'undefined' and 15% occur in areas of post-orogenic extension. Most likely this different behavior has to do with the construction of the Apennines range, which proceeded as a propagating wave causing compression at the front and extension in the back. Under these circumstances, a significant number of gas seeps occur in areas that were undergoing contraction until recently, but that are no longer recording any shortening. For instance, in the northern Apennines (Bonini, 2013: see his Figure 3) documented the existence of number of fossil gas seeps aligned near the boundary between currently active extension, to the SW, and currently active contraction, to the NE.
The actively contracting areas correspond with the hydrocarbon-bearing basins of the northern and central Apennines, of eastern Calabria and southern Sicily, along with their associated Adriatic and Ionian foredeeps; an association that agrees with observations in other productive basins worldwide (Etiope and Klusman, 2002;Van der Meer et al., 2002). It is also expected that many gas seeps and mud volcanoes exist in the Adriatic and Ionian offshore; in fact, some have already been reported in the recent literature (e.g. Donda, et al., 2019), but they are still too sparse to allow statistically significant inferences.
Perhaps surprisingly, the gas seeps and mud volcanoes are not uniformly distributed, but over 80% of them concentrate in two specific regions that are-or used to be, prior to heavy industrial exploitation over the past 70 years-the two largest onshore CH 4 fields countrywide (see red symbols in Figure 6 and Figure 7): a relatively small (∼18,000 km 2 , about 6% of the entire country) rectangular area encompassing the northeastern flank of the northern Apennines and part of the adjacent southern Po Plain, and an even smaller area (∼8,000 km 2 ) of centralsouthern Sicily.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 579390 More specifically, gas seeps and mud volcanoes occur almost exclusively where the pre-Pliocene bedrock-a set comprising variably deformed rocks-is exposed, and are very sparse or absent in the main present-day alluvial basins, e.g., the Po Plain, where the CH 4 reservoir rocks are buried beneath a thick blanket (100 m to over 8 km) of Plio-Pleistocene and Holocene deposits (e.g., Vannoli et al., 2015, with references). This circumstance causes truncation and sealing of potential ascending pathways, but may not prevent CH 4 from being dispersed into the atmosphere through pervasive microseepage. In fact, slow degassing from reservoirs located beneath the southern side of the Po Plain is very likely, as these reservoirs are the largest of the entire hydrocarbon province (e.g., Etiope et al., 2007).
A significant aspect of the distribution of gas seeps is their apparent anticorrelation with the largest earthquakes. Following a statistical study of the distribution of CH 4 fields in the southeastern Po Plain, Mucciarelli et al. (2015) proposed that their productivity is anti-correlated with the presence of seismogenic faults capable of M 5.5+ earthquakes. A visual inspection of Figure 6 shows that the concentration of gas seeps is very low along a ∼110 km stretch of coastal Romagna and Marche, between Cesena and Ancona. The same area shows a concentration of seismogenic faults and of significant earthquakes (Figure 2), and public data from the hydrocarbon industry show that it is devoid of deep wells-and presumably of productive CH 4 fields-whereas they are numerous to the northwest and southeast (https://www.arcgis.com/home/user. html?user DGSUNMIG). The absence of gas seeps can hence be interpreted in terms of lack of CH 4 in the deep reservoirs, and both lines of evidence would be indicative of a larger seismogenic potential of the area under investigation.

CH 4 Emissions and Mud Volcanoes: Structural Control
According to the literature (e.g., Jakubov et al., 1971;Etiope et al., 2007;Bonini, 2012;Etiope et al., 2014;Mazzini and Etiope, 2017;Maestrelli et al., 2019), the reservoir of a mud volcano is generally located over the crest of a growing anticline, provided that the fluids are properly trapped by a sealing rock unit. This configuration is in agreement with the setting of the main fields of mud volcanoes in Italy-Regnano, Montegibbio, Nirano, and Ospitaletto-all of which occur in the hangingwall of the Pedeapenninic Thrust Front, a seismogenic system running along the foothills of the northern Apennines (PTF; DISS Working Group, 2018; Figure 12), and specifically near the crest of the anticlines associated with it. It is also well known that fault intersections act as preferential pathways for deep fluids, ultimately allowing them to easily reach the surface (see Mazzini and Etiope, 2017 and references therein). With specific reference to Italy we remark that: (1) The reservoir-anticlines are usually generated by relatively shallow thrusts. Figure 12 shows that these are concentrated near the nearly emergent PTF, or where the Pliocene deposits are thin or absent, or where the thrusts are very shallow; and (2) A seal-bypass system is required to allow fluids to reach the surface. An efficient superficial path may be provided by extrados fracture arrays developed at the fold crest, i.e. parallel to the anticline axis (e.g., Bonini, 2012;Maestrelli et al., 2019). Nevertheless, we maintain that most mud volcanoes ultimately owe their existence to a deeper path that may be provided by some of the transverse structures described in Data on Active and Seismogenic Faulting Section: well-established, long-lived subvertical lineaments whose activity is testified by frequent background seismicity (see Figure12A) and that may represent the main pathways for fluid expulsion.
In fact, many Italian mud volcano fields are geographically coincident with zones where the dip of the regional monocline changes more drastically (Mariotti and Doglioni, 2000), likely indicating the presence of deep slab tears (Rosenbaum et al., 2008;Rosenbaum and Piana Agostinetti, 2015), which in turn correspond to major chain-orthogonal transverse faults in the upper crust (e.g., Vannoli et al., 2015 for the Po Plain). The adjacent fields of Nirano, Montegibbio, and Ospitaletto are located where the foreland monocline plunges beneath the northern Apennines thrust belt, with a dip that varies laterally between 8°and 23°(see cluster #1 in Figure 12B). This area appears to be a crucial crossroad that controls the extent of forward propagation of the PTF, that advanced further to the east of it and less to the west. Furthermore, this area of deformation locates near a major transverse structure generally known as the Secchia Line (Castaldini et al., 1979;Castellarin et al., 1985), whose role is testified by the sharp western termination of aftershocks following the 2012 Emilia earthquake sequence ( Figure 12A). In turn, the Nirano mud volcanoes fall near a chain-orthogonal fault cutting an anticline that has deformed Plio-Pleistocene clayey sediments (Gasperi et al., 1987;Bonini, 2012).
The mud volcanoes located south of Bologna (San Martino in Pedriolo, San Clemente, Casa Bubano etc.: see cluster #2 in Figure 12B) also concentrate where the PTF intersects the Prato-Sillaro Line ( Figure 12A), another transverse lineament identified as the possible surface of a slab tear (Rosenbaum and Piana Agostinetti, 2015). The Prato-Sillaro Line plays both an active and a passive role, because it is also the locus of earthquakes swarms and occasional moderate-size earthquakes (see Vannoli et al., 2015: Figure 12A).
In summary, the data suggest that most of the largest Italian mud volcanoes locate where an actively growing anticline intersects a primary transverse lineament. Therefore, we propose a systematic, genetic relationship between the existence and location of mud volcanoes and changes in dip of the foreland monocline (see Figure 12C).

DISCUSSION
In this work we investigated in detail the relationships among geofluid emissions, tectonics and seismogenic sources in Italy Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 579390 (refer to Figures 13, 14 for a summary of methods, approaches and results). Our work has the threefold goal of improving the understanding of the relationships between fluid emissions and tectonics in Italy, improving seismic hazard estimates by constraining the location of poorly visible faults, and contributing to the ongoing efforts to calculate the global budget of greenhouse gases. Some of our findings and interpretations challenge those arising from global compilations of geofluids and tectonics, as they stem from relatively high-resolution evidence gathered on Italy: a rather small country, but one endowed with a wealth of geophysical, geochemical and geological data. Over the course of our investigation we have come to realize that the high resolution of Italian data is in fact crucial for discriminating between alternative processes leading to apparently similar outcomes, and hence to possible misinterpretations. It is also crucial for describing the fine geographic grain of the observed phenomena, down to a resolution of less than 1 km: for example, the data show a consistently different behavior between the nature, number and distribution of geofluid emissions in peninsular Italy, which is globally affected by younger and faster tectonic activity, with respect to the more slowly deforming Alps, the Apulian foreland and Sardinia (Figure 2). For this reason, this discussion will refer mostly to the Apennines chain, including its foredeep and foreland areas beneath the Po Plain and the Adriatic and Ionian Seas.
Given the good spatial resolution of Italian data, the nature of the relationships between the different types of geofluid emissions and the ongoing tectonic strains can be analyzed in great detail (refer to Figure 14). First of all, our data and analyses confirm the substantial distinction between the location, behavior and motivations of thermal springs ( Figure 4) and CO 2 emissions ( Figure 5) on the one hand, and gas seeps ( Figure 6) and mud volcanoes (Figure 7) on the other hand. With few exceptions, all geofluid emissions may be assigned to one of these two large geofluid families, and occur in areas currently undergoing extension and shortening, respectively.  Bigi et al., 1992), gas seeps (yellow diamonds) and mud volcanoes (white diamonds) in the northern Apennines (see Figure 7 for location). Active tectonic domains as in Figure 2. Earthquakes of M w ≥ 5.5 are (Continued ) , and main structural elements. The 3D spiral symbols locate gas emissions, mud cone eruptions and H 2 S emissions that were activated following the earthquakes whose year of occurrence is indicated in italics (from CFTI5Med). Sections of the northern Apennines foreland monocline are shown by white lines; numbers indicate its inferred dip (from Mariotti and Doglioni, 2000). The green lines are isopachs of the Pliocene-Pleistocene marine and continental deposits (from Bigi et al., 1992). Panel C -Sketch illustrating the structural position of the mud volcanoes discussed in this work. More specifically, our results indicate that thermal springs and CO 2 emissions dominate in areas of mantle upwelling and crustal stretching, where only minor seismicity is expected owing to the limited thickness of the seismogenic layer; they are found also in the inner Apennines, a region currently undergoing uplift and extension, and specifically along major lithospheric, chain-perpendicular lineaments bounding large normal faults capable of M w 7+ earthquakes. In contrast, CH 4 emissions and mud volcanoes dominate in areas undergoing active or recently active contraction, where no CO 2 emissions are observed and where earthquakes generally do not exceed M w 6.0; in particular, the main mud volcano fields appear to fall where the crests of active anticlines intersect major chain-perpendicular lineaments.   Figure 2 for location), summarizing our main findings. The upper panel summarizes the geodynamic, tectonic and seismological characteristics of the three domains crossed by the section. The central panel summarizes the distribution of all geofluid emissions (i.e., prior to the selection described in the Selecting the Terms of Comparison Section) across the same domains. The lower panel show the topography encountered along the section, and an idealized trend of its Early Pleistocene-to-Present uplift (the southwestern and northeastern maxima reflect the uplift observed in the area of the peri-Tyrrhenian volcanism and along the crestal portion of the Apennines respectively).
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 579390 20 An overarching-and largely unexpected-conclusion arising from our elaborations is that the geographic pattern and the nature of most geofluid emissions does not depend on the current geodynamic regime and on its spatial variability, but rather on the deeper architecture of the crust that results from mostly pre-Pliocene tectonic events. There is no doubt that the large-scale geographic pattern of most geofluid emissions is controlled by Italy's current tectonic setting; its fine-grain, however, exhibits little concordance with ongoing surface processes, being rather controlled by paleogeography and by the existence of buried weaknesses extending to lower crustal or even lithospheric depth. In other words, what is crucial for determining the geographic pattern and the nature of geofluid emissions is not only the nature of the current tectonic strains, but rather their integral over a 5-10 My-long geodynamic history. The only exception to this rule is the peri-Tyrrhenian volcanic and geothermal district (Figure 2), which includes most of coastal and central Tuscany and Latium, plus selected portions of coastal Campania: a region that has experienced crustal stretching and mantle upwelling in Pliocene and Quaternary times.
More specifically, most of our selected geofluid emissions (see Selecting the Terms of Comparison Section for a discussion of selection criteria) take advantage of the existence of transverse lineaments discussed earlier (e.g., Data on Active and Seismogenic Faulting Section). All of these lineaments predate the onset of the current extensional regime active over most of the Apennines, but are coeval with the initial phases of the Miocene-Pliocene buildup of the chain. Many of them are parallel to the current trend of extension, or form a small angle to it, and therefore should be only mildly active or totally inactive under the extensional tectonic regime that currently dominates in the core of the Apennines. But the existence of thermal springs and CO 2 emissions in the inner chain suggests that the transverse lineaments are continuously kept alive by one-or more likely, by a combination-of the following ongoing processes: (1) the response of the upper crust to the activity of the slab tears, which dominates in the northern Apennines and is reflected in occurrences such as those described in the test cases of Garfagnana (Test Case 1: The Garfagnana Region, Northern Apennines Section: Figure 10), centralnorthern Apennines (Test Case 1: Central-Northern Apennines Section: area marked with A in Figure 5), and northern Apennines (CH 4 Emissions and Mud Volcanoes: Structural Control Section: Figure 12C); (2) the progressive regional uplift of the Apennines chain, which takes place with a bell-shaped pattern peaking at nearly 2.0 mm/y near the crest and slowly tapering to zero toward the Tyrrhenian and Adriatic coasts (e.g., Bordoni and Valensise, 1998;Mancini et al., 2007: Figure 14). It is generally accepted that the regionalscale uplift of the Apennines is closely related with the ongoing extension, and hence with the occurrence of large earthquakes near the chain crest, but also that the rate of uplift may be different in adjacent portions of the chain, possibly as a result of changes in the geometry of the Apennines subduction caused by slab tears. The Olevano-Antrodoco-Sibillini Line (OASL) test case (Test Case 2: the Olevano-Antrodoco-Sibillini Line, central Apennines Section and Figure 11) provides an example of a regional structure that forms an angle of 45-60°with the current direction of extension in the central Apennines, and consequently should not be active in the current stress regime. Yet, the frequent low-to mid-level seismicity associated with the OASL and other similar lineaments suggests that differential regional uplift is likely responsible for their progressive rejuvenation-albeit in an indirect fashion-and hence for the rejuvenation of the pathways feeding geofluid emissions. This process, which certainly dominates in the central Apennines, most likely contributes to rejuvenating transverse lineaments throughout the chain as far as Liguria, to the north, and Sicily, to the south, and may also account for the limited number of thermal springs and CO 2 emissions found in the core of the Alps; (3) the ongoing activity of transverse lineaments segmenting the Adriatic side of the northern Apennines accompanies the residual activity of both the thrust fronts expressed at the surface and of those buried beneath the Po Plain (Panels B and C in Figure 12); (4) the right-lateral reactivation of deep E-W lineaments beneath the central and southern Apennines and in the Apulian foreland, and of left-lateral NNW-SSE lineaments in the Hyblean foreland; a primary effect of the residual NNW-oriented push of the African Plate with respect to stable Eurasia (areas marked with B, C in Figure 5).
In summary, tectonic stretching is not the primary reason for the existence of emissions belonging to the 'extensional geofluid family' in the core of the Apennines, although it is the indeed the reason for the generation of large normal-faulting earthquakes. We suggest that the ongoing extension, the seismicity and the release of geofluids are simultaneous but independent outcomes of the ongoing regional uplift of the chain: a process ultimately caused by the buoyancy of the orogen following the retreat of the underlying subducting slab and laterally modulated by the existence of slab tears. This finding may require revision of previous hypotheses (e.g., Chiodini et al., 2004;Chiodini et al., 2020) regarding the relationships between CO 2 degassing and seismicity in the central Apennines.
The close spatial correlation between CO 2 emissions belonging to the 'extensional geofluid family' and deep discontinuities suggests processes originating at depth, possibly below the brittle-ductile transition (e.g., Frezzotti et al., 2009, and references therein), and hence that the CO 2 is at least partly mantle-derived. This finding is further supported by the spatial anticorrelation between these emissions and the footprint (i.e., the projection onto the topographic surface) of the largest extensional faults that exist in peninsular Italy, which generate large earthquakes along the main axis of the Apennines. Extensional geofluid emissions almost never occur within the surface projection of the largest active extensional faults: rather, they concentrate at their edges, hence in the gaps between one fault and the next one in the same system, generally referred to as "segment boundaries" in the seismotectonic literature. Such boundaries are assumed to be permanent properties of the Earth's crust, at Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 579390 least at the scale of one to a few My, ultimately leading to the characteristic earthquake concept (e.g., Schwartz and Coppersmith, 1984) and to its implications for the location and size of future earthquakes. A further finding is that most gas seeps-even those fed by thermogenic methane-originate at shallow (1-10 km) or even very shallow depth (1-3 km). Our observations confirm that the presence of gas seeps is a direct indication of underlying hydrocarbon reservoirs, combined with limited permeability of the overlying crust. Conversely, their absence may indicate that they gas reservoirs never developed, or that the gas has already been lost to the atmosphere due to enhanced faulting and fracturing.

CONCLUSION
Earthquakes causing significant shaking and permanents strains in a large rock volume allow the generation of new, episodic fluidflow regimes, which may become permanent when rocks having limited permeability are pervasively fractured by the upward propagation of seismogenic faults. Such geophysical conditions are easily attained near the edges of seismogenic faults, where transverse lineaments are known to occur from independent lines of evidence (e.g., Curewitz and karson, 1997).
Our Italy-based hypotheses have important implications for the understanding of deep fluid circulation in active regions and for improving the quantification of natural CO 2 and CH 4 emissions. We have shown that the location of CO 2 emissions in the inner portions of the Apennines is not controlled by ongoing continental extension, but rather by other large-scale processes that are responsible "also" for rather localized extension. In contrast, the extrados faults that allow the development of gas seeps and mud volcanoes should imply ongoing activity of the underlying thrusts; alternatively, fluids might accumulate also at the core of inactive anticlines, where they may escape to the atmosphere through cross-fold joints.
The demonstrated interplay between large Apennines normal faults and transverse lineaments may also help constraining the location and geometry of blind and hidden seismogenic sources that remain undetected in the Italian territory. The relatively small earthquakes (generally ≤ M w 5.0) that are spatially coincident with geofluid emissions in the core of the Apennines may be interpreted as generated by passive reactivation of the transverse lineaments under the current stress/train regime: a circumstance that acts as an upper boundary to their size, and allows for a more confident identification of the actual extent of the main chain-parallel seismogenic sources.
Our work provides also a detailed framework for exploring the preparation of large earthquakes based on geofluid emissions, potentially contributing to the identification of sites that are most promising for earthquake prediction investigations (e.g., Martinelli, 2020).
All the observations and hypotheses made in this paper were made possible by the systematization of knowledge concerning widely separated manifestations of the activity of the Earth. Although the databases we relied on are continuously improved, we maintain that they already form and exceptional body of information, probably unique worldwide, making our inferences reasonably robust.

DATA AVAILABILITY STATEMENT
This article was made possible by a great deal of published data, that actually represent the paper's main strength. The vast majority of these data are made available through Web-GIS databases, listed in Supplementary Table S1 along with the corresponding URL. All these databases were fully functioning on the date of submission of the second version of this article. Since the most relevant of these databases (database of Thermal Springs, by Martinelli and Albarello, 1997) is also the oldest and the most difficult to access, we reported it in the Supplementary Material section for the convenience of all readers, after correcting a few typos the existed in the original version. 1 | Summary of the distribution of geofluid emissions collected for this paper in the four tectonic districts discussed in Selecting the Terms of Comparison Section. For each emission type the first line reports all available data, whereas the second line (in italics) reports only the data left following the selection discussed in the Selecting the Terms of Comparison Section. Numbers in parentheses indicate the % of each type of emissions falling in each tectonic district, relative to the total (last column to the right). Notice that some datapoints were deliberately double-counted as they fall in overlapping regions. a Includes one datapoint falling into the region of overlap between thrusting and strike-slip, and two datapoints falling into the region of overlap between post-orogenic extension and strike-slip. b Includes one datapoint falling into the region of overlap between post-orogenic extension and strike-slip.