Active degassing of crustal CO2 in areas of tectonic collision: A case study from the Pollino and Calabria sectors (Southern Italy)

Randazzo, Paolo; Caracausi, Antonio; Aiuppa, Alessandro; Cardellini, Carlo; Chiodini, Giovanni; Apollaro, Carmine; Paternoster, Michele; Rosiello, Angelo; Vespasiano, Giovanni

doi:10.3389/feart.2022.946707

ORIGINAL RESEARCH article

Front. Earth Sci., 30 September 2022

Sec. Geochemistry

Volume 10 - 2022 | https://doi.org/10.3389/feart.2022.946707

Active degassing of crustal CO₂ in areas of tectonic collision: A case study from the Pollino and Calabria sectors (Southern Italy)

Carmine Apollaro⁵

Michele Paternoster^2,6

Angelo Rosiello³

Giovanni Vespasiano⁵

¹Department of Earth and Marine Sciences (DiSTeM), Università di Palermo, Palermo, Italy
²Section of Palermo, National Institute of Geophysics and Volcanology (INGV), Palermo, Italy
³Department of Physics and Geology, Università degli Studi di Perugia, Perugia, Italy
⁴Section of Bologna, National Institute of Geophysics and Volcanology (INGV), Bologna, Italy
⁵Department of Biology, Ecology and Earth Sciences (DiBEST), University of Calabria (UniCAL), Arcavacata di Rende, Calabria, Italy
⁶Department of Sciences, University of Basilicata, Potenza, Italy

Carbon dioxide (CO₂) is released from the Earth’s interior into the atmosphere through both volcanic and non-volcanic sources in a variety of tectonic settings. A quantitative understanding of CO₂ outgassing fluxes in different geological settings is thus critical for decoding the link between the global carbon budget and different natural processes (e.g., volcanic eruption and earthquake nucleation) and the effects on the climate evolution over geological time. It has recently been proposed that CO₂ degassing from non-volcanic areas is a major component of the natural CO₂ emission budget, but available data are still sparse and incomplete. Here, we report the results of a geochemical survey aimed at quantifying CO₂ emissions through cold and thermal springs of the tectonically active Pollino Massif and Calabrian arc (Southern Italy). The chemical ad isotopic (He and C) composition of fifty-five dissolved gas samples allows to identify two different domains: 1) a shallow system dominated by gas components of atmospheric signature (helium, hereafter He) and biogenic origin (C), and 2) a deeper system in which crustal/deep fluids (CO₂ and He) are dominant. The measured He isotope ratios range from 0.03 to 1.1 Ra (where Ra is the He isotopic ratio in the atmosphere) revealing a variable atmospheric contamination. Furthermore, the He isotopic data indicate the presence of traces of mantle He contributions (2%–3%) in the thermal groundwater. The prevailing low R/Ra values reflect the addition of crustal radiogenic ⁴He during groundwater circulation. Using helium and carbon isotope data, we explore the possible sources of fluids and the secondary processes (dissolution/precipitation) that act to modify the chemistry of pristine volatiles. For the thermal springs, we estimate a deep C output of 2.3 x 10⁷ to 6.1 x 10⁸ mol year⁻¹. These values correspond to deep CO₂ fluxes per square km comparable with those estimated in several active and inactive volcanic areas and in continental regions affected by metamorphic CO₂ degassing (e.g., the southern margin of the Tibetan Plateau).

1 Introduction

The current rise in atmospheric CO₂, and its link with the global climate change, provides a strong motivation to understand the natural processes that control the nature and magnitude of geological CO₂ cycling (Evans, 2011). The release of carbon dioxide into the atmosphere via Earth degassing has played a crucial role in controlling global planetary temperature over geological time via the greenhouse effect (Foster et al., 2017). The modes and rates of geological CO₂ release are thus crucial to understand the compositional evolution of the atmosphere through geological time, life on Earth, and climate changes (Berner and Lasaga, 1989; Kerrick, 2001; Dasgupta, 2013; Aiuppa et al., 2019; Fischer and Aiuppa, 2020; Guo et al., 2021). Despite continuous improvements via direct measurements, models, and global extrapolations, the CO₂ Earth degassing output remains poorly constrained, hampering full understanding of the geological carbon cycle (Berner and Lagasa, 1989; Burton et al., 2013; Fischer, 2013; Fischer et al., 2019; Fischer and Aiuppa, 2020). The release of CO₂ from the Earth’s interior into the atmosphere occurs in different tectonic settings (Lee et al., 2019), through volcanic and non-volcanic sources, and on a global scale, it is known that CO₂ discharges are associated with tectonically/seismically active zones (Barnes et al., 1978; Chiodini et al., 2004; Tamburello et al., 2018). Quantitative estimates of CO₂ outgassing fluxes in different tectonic settings are thus critical for decoding the link between the global carbon budget and climate evolution from a whole-Earth carbon cycling perspective (Zhang et al., 2021). In the last decades, the number of studies on CO₂ degassing in non-volcanic areas has risen exponentially, emphasizing the important contribution of these areas to the Earth carbon budget (e.g., Chiodini et al., 2020, 2004; Minissale, 2004; Becker et al., 2008; Italiano et al., 2008; Groppo et al., 2022, 2017, 2013; Rolfo et al., 2015; Lee et al., 2016; Tamburello et al., 2018; Caracausi and Sulli, 2019; Frondini et al., 2019). The first regional-scale CO₂ Earth degassing studies led to the catalogue of Italian CO₂-rich gas emissions (googas.ov.ingv.it and www.magadb.net) and to the regional map of deeply derived CO₂ degassing in central Italy that uses the quantification of carbon dissolved in regional groundwater systems (Chiodini et al.,. 2000; 2004; 2011). Some studies (Chiodini et al., 2004, 2020; Miller et al., 2004) also demonstrated a relation between CO₂ degassing and seismogenesis in the Italian Apennines, pointing to the presence of gas triggering earthquakes. The Mt. Pollino region, at the southern end of the Apennines (southern Italy), has been historically recognized as one of the most hazardous seismic gaps in the intra-Appenine seismogenic belt (Napolitano et al., 2021), but it has recently been affected by seismic sequence occurred between 2010 and 2014 and characterized by about 10,000 earthquakes with highly variable magnitude (strongest events ML 4.3 and ML 5.0; De Matteis et al., 2021; Pastori et al., 2021). Moreover, recent studies identified fluid-related dynamics responsible for historical and recent seismicity of the area (Sketsiou et al., 2021). The Calabrian arc, further to the south, is one of the most active seismogenetic areas in Italy (Italiano et al., 2010; Neri et al., 2020), which has been repeatedly affected by catastrophic seismic events with 5.9 < M < 7.2 during the last centuries (18 times from 1626 to 1908; Gruppo di Lavoro CPTI, 2004; Boschi et al., 2000). The two areas are characterized by the presence of several springs, some representing low-enthalpy geothermal resources (Zarlenga, 2011; Vespasiano et al., 2014, 2015a, 2015b, 2015c, 2016, 2021; Apollaro et al., 2015, 2016, 2020). The geochemical and isotopic compositions of Calabrian and Pollino waters have previously been investigated to define their geochemical features and geothermal potential (Bencini and Ciracò, 1982; Duchi et al., 1991), to investigate a link with seismicity and implications for a fluid–fault relationship (Gurrieri et al., 1984; Calcara and Quattrocchi, 1993; Italiano et al., 2010; Apollaro et al., 2020) and to evaluate potential natural metal contamination of spring waters (Margiotta et al., 2012, 2014; Paternoster et al., 2021). However, no attempt has been made so far to model the water–gas interaction processes and to quantify the regional-scale budget of CO₂ sequestrated/transported by aquifers at depths and released into the atmosphere upon spring discharge.

In this study, we present the results of a geochemical study of cold and thermal springs from both the Calabrian arc and the Pollino region. Our goals are to 1) investigate the relationships between Earth degassing and geological features in the two areas; 2) assess the presence and eventual origin of deep volatiles released in the hydrothermal basins and the surrounding areas; 3) model the processes at depths that can modify the pristine chemistry of deeply rising volatiles, potentially affecting the deep carbon budget; and 4) estimate the total deeply derived CO₂ output. For this aim, we combine helium isotopes (³He/⁴He), total dissolved inorganic carbon (TDIC), and dissolved carbon isotopes (δ¹³C_DIC) of groundwaters to explore the origin of carbon and to develop a model of water–gas–rock interaction. The results are then compared with the CO₂ output from some active tectonic regions and volcanic areas worldwide.

2 Geological and hydrogeochemical background

The Calabrian–Peloritan Orogen (CPO) is a well-developed, arc-shaped segment of the circum-Mediterranean orogenic belt between the southern Apennines and the Sicilian Maghrebides, bounded by two main tectonic lineaments: the Sangineto line to the north and the Taormina line to the south (Cirrincione et al., 2015; Tortorici, 1982a, Tortorici 1982b; Figure 1). Incorporation of the Calabria terranes into the Apennine–Maghrebian chain is related to the processes responsible for the formation of the Tyrrhenian Basin since the late Miocene (Alterberger et al., 2011 and reference therein). In this context, the Calabrian arc represents an accretionary wedge, caused by the collision of the Eurasian and African plates (Amodio Morelli et al., 1976; Tortorici, 1981), consisting of a series of ophiolite-bearing tectonic units (Liguride Complex; Ogniben, 1969) and overlying basement nappes (Calabride Complex; Ogniben, 1969) with Paleozoic metamorphic and plutonic terranes that represent the remnants of Caledonian, Hercynian, and Alpine orogens (e.g., Amodio Morelli et al., 1976; Schenk, 1981; Zanettin Lorenzoni, 1982; Atzori et al., 1984; Del Moro et al., 1986; Zeck, 1990; Messina et al., 1994). The CPO is classically subdivided into a northern and southern sector, separated in correspondence to the Catanzaro Strait Basin, a Neogene–Quaternary basin connecting the Ionian and Tyrrhenian seas (Tortorici, 1982a; Chiarella et al., 2012, 2016; Longhitano et al., 2014; Brutto et al., 2016). The two sectors differ for the structural style and assemblage of the chain. The northern block exhibits overthrust of alpine and pre-alpine crystalline units on carbonate tectono-stratigraphic units. In the southern block, the chain is made up of alpine and pre-alpine crystalline units, while the Apennine carbonate rocks are not present beneath the crystalline–metamorphic units (Apolllaro et al., 2019a and reference therein). In particular, the Calabria arc terrane consists of three main groups of stacked tectonic units (Tursi et al., 2021) that can be summarized, from bottom to top, as follows: 1) the Lower Complex, characterized by Apennine units with Meso-Cenozoic phyllites and partly metamorphosed carbonate rocks exhibiting high pressure (̴1.4 GPa) and low temperature (̴390°C) metamorphic imprint (Iannace et al., 2007); 2) the Intermediate Complex, composed of ophiolite units of the Ligurian Tethys’ oceanic lithosphere (Liberi et al., 2006), which records HP/LT Eocene metamorphism with peak conditions at 2.0–2.1 GPa and 470–490°C (Tursi et al., 2020); and 3) the Upper Complex, which consists of the Hercynian continental crust, showing a local Alpine metamorphic overprint at 0.3–0.7 GPa and 200–450°C in the Sila Massif and Catena Costiera (Piccarreta, 1981; Acquafredda et al., 1994; Graessner and Schenk, 2001; Liberi et al., 2011; Ortolano et al., 2020) and up to 1.1–1.2 GPa and 540–570°C in the Aspromonte Massif (Cirrincione et al., 2008). According to the current geodynamic models, the evolution of the Calabrian arc was driven by the south-eastward retreat of the Ionian slab (Malinverno and Ryan, 1986; Jolivet and Faccenna, 2000; Faccenna et al., 2001). During the Eocene, subduction of the Ligurian Tethys oceanic crust underneath the continental margin, represented by the Calabria terrane, (Stampfli and Borel, 2002; Rossetti et al., 2004; Vitale et al., 2019), is thought to have occurred at 47–20 Ma (Borsi and Dubois, 1968; Schenk, 1980; Beccaluva et al., 1981; Thomson, 1994, 1998; Rossetti et al., 2001, 2004; Shimabukuro et al., 2012). Currently, the active subduction residue of the ancient, 200-km wide, subducting slab dipping 70° towards NE is found beneath the Calabrian arc with the presence of deep seismicity (150–300 km) (Lucente et al., 1999; Spakman and Wortel, 2004; Chiarabba et al., 2008; Neri et al., 2009; Neri et al., 2012). Different studies show a rapid deepening of the Ionian Moho beneath Calabria, illustrating the geometry of the subduction zone (Piana Agostinetti et al., 2009; Scarfì et al., 2018). The estimated current plate convergence velocity between the two plates is 3–5 mm/year (Neri et al., 2020; Mattei et al., 2007), and the rollback of the subducting slab occurs at about 2 mm/year rate (Hollenstein et al., 2003; Devoti et al., 2008; Nocquet, 2012).

FIGURE 1

FIGURE 1. (A) Simplified geological map of the Calabrian arc and the surrounding region (modified after Amodio-Morelli et al., 1976 and Bonardi et al., 1988b) with (B) the geological section on the bottom (after Van Dijk and Scheepers, 1995; and Van Dijk et al., 2000; modified) and on top (C), the location of the study area and tectonic simplified sketch of the Calabrian arc (after Tansi et al., 2005, modified). TF, Taormina Fault; CTF, Catanzaro Trough Fault; SF, Sangineto Fault. For a detailed overview of the individual springs, see Paternoster et al. (2021), Apollaro et al. (2020), Apollaro et al. (2019c), Apollaro et al. (2012), Vespasiano et al. (2021).

The Calabrian arc is one of the strongest seismic areas in Italy (Neri et al., 2020) has experienced several destructive earthquakes with estimated magnitudes of about 7 or higher (Scarfì et al., 2018). After the destructive 1908 earthquake (Rovida et al., 2016), a few events with M > 4 and about 200 shocks with a magnitude between 3 and 4 (out of a total of 3,800 events) have occurred between 1980 and 2005 (Gruppo di Lavoro CPTI, 2004; Castello et al., 2006). Crustal thickness reaches about 35–38 km in correspondence of the highest portion of the chain (Di Stefano et al., 2009), and the recorded seismicity is marked by focal depths <30 km (i.e., crustal depths; Neri et al., 2020; Boschi et al., 2000). Since the Middle Pleistocene, an intense WNW–ESE-oriented regional extensional phase occurred, resulting in a longitudinal faults system with NNE–SSW strikes and parallel to the mountain system (Figure 1A), consisting of a 10–50 km-long distinct normal fault segment running along the western side of the Calabrian arc (Tansi et al., 2005; Catalano et al., 2008). The development of the rift-zone, coupled with contrasting vertical movements, such as mountain chain uplifting of 0.5–1.2 mm/year in the last 1–0.7 Myr (Monaco et al., 1996; Ferranti et al., 2008; Faccenna et al., 2011), is still an active process is (Dumas and Raffy, 2004) and probably represents the response to the isostatic rebound due to the detachment of the Ionian-subducted slab (Tortorici et al., 2003; Tortorici et al., 1995 Westaway, 1993). GPS data show differential motion of the Calabrian arc relative to both the Nubia and Eurasia plates, which causes active extension on the region with the developments of the aforementioned extensive faults (Mattei et al., 2007 and reference therein). These normal faults are considered to be major seismogenetic faults (Monaco et al., 1996; Neri et al., 2006), with the NE trending fault systems of the Messina Straits, Gioia Basin, and Mesima Valley believed to have generated the major earthquakes of the area (Rovida et al., 2019, 2020; Neri et al., 2020).

The first comprehensive geochemical data-set for fluids circulating over the Calabrian arc has been presented by Italiano et al. (2010). The authors used the chemical and isotopic (C and He) compositions of groundwater and dissolved gas to show the aquifers contain deeply derived CO₂-rich gas with a radiogenic He signature, a consequence of long residence in the crust. Clear fault–fluid relationships have been found in some of the investigated sites, with the thermal character of the investigated waters being linked with deeper hydrological circuits and normal geothermal gradients (30°C/km). Recent studies have reconstructed the conceptual geothermal model of some Calabrian sites using a “site-specific” (Vespasiano et al., 2021), multidisciplinary approach, involving geological, hydrogeological, and geochemical data. These studies highlighted different features between the thermal waters from north to south, related to complex geologic and tectonic settings of the region (e.g., Apollaro et al., 2012, 2016, 2019a, 2019b, 2019c, 2020, 2021; Vespasiano et al., 2012, 2014, 2015a, b, 2015c, 2016). Moreover, a deep component for dissolved gases of Pollino water have been identified from Apollaro et al. (2020) on the base of their C and He isotopic compositions.

3 Sampling and analytical methods

In total, 55 water samples were collected (see Figure 1) during two field campaigns in February and July 2019 (Table 1). Water temperature, pH, Eh, and electrical conductivity (EC) were measured in situ by means of high-resolution multi-parametric probes (Hanna Instruments HI-9828). Total alkalinity was measured in situ by acidimetric titration with 0.05 N HCl using methyl-orange as the indicator. The water samples were filtered in situ through a 0.45-μm pore-size membrane and acidified with supra-pure HNO₃. Different sample aliquots (one filtered and two filtered and acidified) were collected. All samples were stored in high-density polyethylene bottles for laboratory analysis. Major elements were determined by high-performance liquid chromatography (HPLC) by using a Thermo Scientific Dionex™ ICS-1100 equipped with Dionex IOnPac AS23 and Dionex IonPac CS12A columns for the determination of anionic (F⁻, Cl⁻, SO₄²⁻, and NO₃⁻) and cationic (K⁺, Na⁺, Ca²⁺, and Mg²⁺) species, respectively. The computed charge balance resulted <5% in all water analyses. All measurements were performed at the Department of Biology, Ecology and Earth Sciences laboratories of the University of Calabria (Cosenza, Italy). TDIC and saturation indexes (SI) with respect to the mineral phases (calcite, dolomite, and gypsum) were calculated using the PHREEQC Interactive computer code (Parkhurst and Appelo, 1999), considering the measured pH, the temperature, and the water composition.

TABLE 1

TABLE 1. Location, physico-chemical parameters, flow rates, chemical and isotopic compositions of Calabria and Pollino springs. T—thermal; C—cold; < d.—below detection limits; - not measured.

The water samples used for the analysis of dissolved gases, δ¹³C of TDIC, He and Ne isotopes (³He, ⁴He, and ²⁰Ne) were sampled in glass bottles according to Capasso and Inguaggiato (1998) and analyzed in a few days from their collection to prevent any contamination and/or loss of volatiles. The chemical composition of the dissolved gases was analyzed by using the method described by Capasso and Inguaggiato (1998), which is based on the equilibrium partition of gas species between a liquid and a gas phase. The analysis was performed utilizing a Perkin Elmer Clarus 500 gas chromatograph equipped with a 3 m packed column (100/120 Shincarbon, Ar gas carrier) and two detectors, a thermal conductivity detector (TCD), and a flame ionization detector(FID), and using Ar as the carrier gas. H₂, O₂, N₂, and CO₂ were measured by means of the TCD, while CH₄ and CO were determined through a FID coupled with a methanizer. Analytical errors for CO₂, N₂, H₂, CO, CH₄, and O₂ are within 3%. Analyses of the dissolved noble gases (He and Ne) and He isotopic composition (³He/⁴He) were performed by using the methodology proposed by Inguaggiato and Rizzo (2004), which is based on the isotope equilibrium between the liquid and host gas phases (e.g., N₂). The extracted gases from waters are purified in a high-vacuum purification line that is directly connected to the mass spectrometers (Rizzo et al., 2019 and references therein). He and Ne isotopes are analyzed using a static vacuum mass spectrometer (GVI Helix SFT) with a double collector to detect ³He and ⁴He ion beams simultaneously with a multi-collector Thermo-Helix MC Plus mass spectrometer (isotopic ratio precision within ±0.5%.). The ³He/⁴He ratio was determined by measuring ³He in an electron multiplier detector and ⁴He in an axial Faraday detector. The isotopic composition of TDIC (δ¹³C_TDIC) is measured by using the method proposed by Capasso et al. (2005), using a Thermo Scientific Delta V Advantage continuous flow isotope ratio mass spectrometer. All δ¹³C_TDIC values were reported relative to the Vienna Pee Dee Belemnite (VPDB) international reference standard, and the analytical precision is ±0.15‰. All sampling and analytical devices are provided by the Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Palermo.

4 Results

The physico-chemical parameters of the collected waters are reported in Table 1 together with the computed TDIC values, chemical and isotopic compositions (He and C) of dissolved gases. Major ions water chemistry is presented in Supplementary Table S1. The 55 samples have been subdivided in two categories on the basis of their discharge temperature following the classification of Apollaro et al. (2020) and Italiano et al. (2010): cold (T<20°C) and thermal (T≥20°C). The cold waters show compositions typical of shallow air-saturated waters (ASW) with N₂ concentrations ranging between 10.3 and 17.6 ccSTP/l and CO₂ from 1.5 to 24.4 ccSTP/l. O₂ concentrations are between 2.6 and 8.6 ccSTP, while CH₄, CO, and H₂ concentrations are very low (CH₄≤5.5×10^–04 ccSTP/l; CO≤3.5×10^–4 ccSTP/l; and H₂≤1.6×10^–3 ccSTP/l). Both He and Ne are present in trace amounts (up to 9.4×10^–5 ccSTP/l and 2.9×10^–4 ccSTP/l). In thermal waters, the measured He amounts are much higher (up to 0.11 ccSTP/l) than in cold samples, while Ne ranges between 1.4×10^–4 and 4.3×10^–4 ccSTP/l (Table 1).

In the ternary diagram CO₂−N₂−O₂ of Figure 2, thermal samples cluster along the CO₂−N₂ axis, as implied by their low to negligible O₂ contents (≤1.4 ccSTP/l). These samples have CO₂ concentrations from 0.1 to 36.6 ccSTP/l and N₂ concentrations from 0.004 to 22.6 ccSTP/l. Sample RC is the unique thermal sample falling close to the ASW field (Figure 2). This is because the sample was taken from a small stream exposed to the atmosphere and thus differs from other thermal samples from underground circuits not exposed to the atmosphere.

FIGURE 2

FIGURE 2. CO₂–O₂–N₂ ternary diagram for the dissolved gas from cold and thermal springs. Thermal samples (red) fall along the axis between CO₂ and N₂, as implied by their being the two dominant gas species, showing low O₂ values. The cold waters (blue circles) show compositions typical of shallow water falling on mixing line between the air-saturated water (ASW) and CO₂-rich end member. The label Latronico Group refers to a group of samples from the Basilicata–Calabria border (Lucane thermal area). The cold group includes samples RMC, SP, SF, SC1, and SC2. The thermal group includes samples FS, SC, SXC, RC, and SCE. More information about these samples can be found in Apollaro et al. (2020).

The δ¹³C_TDIC ranges from −14.8‰ to −5‰ (vs. V-PDB) in cold samples and from −9.9‰ to −4.4‰ in thermal samples (Table 1).

The ³He/⁴He isotopic ratio (R) of each sample is normalized to the same ratio in air (R_a=1.386∙10^–6; Ozima and Podosek, 2002). The R/R_a ratios range from 0.8 to 1.1 and from 0.03 to 0.5 R_a, for cold and thermal waters, respectively. The ⁴He/²⁰Ne ratios are up to 659 for thermal samples, well above the ASW ratio (0.295 at 25°C; Ozima and Podosek, 2002), indicating a negligible atmospheric contamination. On the contrary, the ⁴He/²⁰Ne ratios of the cold waters, ranging between 0.26 and 0.35, indicate a dominant atmospheric derivation (Figure 3).

FIGURE 3

FIGURE 3. ³He/⁴He ratios (expressed as R/R_a) versus ⁴He/²⁰Ne ratios. All samples fall along mixing lines between three possible end-members characterized by distinct He isotopic signatures: 1 Ra, for air-saturated water (ASW; Ozima and Podosek, 2002); 0.01–0.02 Ra, for pure crustal fluids dominated by radiogenic ⁴He produced by U and Th decay (Ballentine and Burnard, 2002); and 6.1 ± 0.9 Ra, for the European subcontinental lithospheric mantle, ESCLM (Gautheron & Moreira, 2002) and ⁴He/²⁰Ne values (⁴He/²⁰Ne ratios >1,000 for crust and mantle and 0.295 for ASW, respectively; Sano et al., 1985; Ozima and Posek, 2002). Thermal samples show helium isotopic composition near the crustal value of 0.02 R_a with small percentages in mantle contribution (2%–3%) and negligible atmospheric contamination, while cold waters and only one thermal sample show ASW-like composition. Data for comparison (small diamonds) from Italiano et al. (2010).

5 Discussion

In the following sections, we discuss the sources of fluids (both CO₂ and He) in the studied areas, and the secondary processes that affect the groundwater He–CO₂ signature during circulation and storage in the aquifers.

5.1 Helium

In natural fluids, He is typically fed by three distinct sources: the mantle, the crust, and the atmosphere (Sano et al., 1997). Each of these sources has a distinct He isotopic signature and ⁴He/²⁰Ne ratio (>1,000 for crust and mantle, 0.295 for air-saturated water at 25°C; Sano and Marty, 1985; Ozima and Podosek, 2002). Therefore, the contributions of these three different sources can be solved by using binary mixing equations. Applying the approach of Sano et al. (1997), and assuming that all ²⁰Ne is of atmospheric origin, we estimated low atmospheric contributions (<10%) for the thermal samples, along with small percentages (up to 2%–3%) of mantle He contribution, in agreement with the data reported by Italiano et al. (2010) (Figure 3). Dissolved gases from the cold waters and one thermal sample (Cotr.S) have an ASW-like composition, which indicates high atmospheric contamination probably due to a shallow hydrological circuit. We find a statistically significant (R²=0.86) negative correlation between the He isotopic composition and its concentration (Figure 4). The lower R/Ra values of the He-rich thermal waters with respect to that of the cold waters is thus explained by the addition of crustal He, rich in radiogenic ⁴He produced by U and Th decay (Ballentine and Burnard, 2002), during deep/prolonged circulation in the crust. It is noteworthy that the He-enriched samples have been collected in areas geologically dominated by metamorphic rocks characterized by high U and Th concentrations (e.g., 3.3 ppm of U and 19.4 ppm of Th for the “Sila” gneiss; Micheletti et al., 2007 and 297 ppm of U and 155 ppm of Th for zircons from the “Catena costiera” gabbros; Liberi et al., 2011). Hence, it is reasonable that these low He isotopic ratios reflect, in addition to possible long residence times, also the high radiogenic ⁴He production in such U-Th-rich lithologies.

FIGURE 4

FIGURE 4. ³He/⁴He vs. He concentrations. A good correlation between the helium isotopic composition and helium concentrations (R²=0.86) is found. The lower R/R_a values can be explained by an addition of crustal radiogenic ⁴He to the thermal waters that also shifts the He isotopic ratio from the ASW value toward the crustal radiogenic end member (0.01–0.02Ra). In this scenario, the samples which show higher He values and lower ³He/⁴He could be interested by longer residence periods in the crust. Data for comparison (small diamonds) from Italiano et al. (2010). Blue = cold waters; red = thermal waters.

5.2 Carbon

The relationship between the TDIC and δ¹³C_TDIC (Figure 5) can provide additional constraints on the sources of fluids. Indeed, deeply rising fluids ascending through the crust interact with rocks and groundwaters that cause changes in carbon abundance and its isotopic composition (Randazzo et al., 2021). From the carbon mass balance approach developed by Chiodini et al. (2000, 2020), we estimate, for each sample, the external carbon contribution, C_ext, (i.e., the C fractions not resulting from carbonate rock dissolution) and its isotopic composition δ¹³C_ext (see Supplementary Appendix A).

FIGURE 5

FIGURE 5. C_ext vs. δ¹³C_ext diagram. Cold waters show the lowest C_ext and δ¹³C_ext values with a negative correlation between the two variables due to the dissolution of isotopically “light” biogenic CO₂ falling along the mixing lines (green) between the infiltration water and the biogenic CO₂ end members (δ¹³C of −20‰, −24‰ and −28‰), while thermal waters are characterized by higher δ¹³C_ext values linked to the presence of deep CO₂ source end member. The theoretical curves have been computed considering C_inf (carbon linked to infiltration processes) contents of 1, 2, and 3 mmol/L and carbon isotopic composition for deep end member of −1‰ (black lines).

The relationship between C_ext and δ¹³C_ext suggests the presence of three distinct carbon sources in the groundwaters (Figure 5). The cold waters fall along the mixing lines (green lines in Figure 5) between a meteoric component (C_ext= 0.03 mmol/L; δ¹³C_ext= −7‰ to −5.5‰; Apollaro et al., 2020; Chiodini et al., 2011) and a set of end members whose isotopic compositions are in the range of biogenic CO₂ (i.e., δ¹³C from −20‰ to −28‰; Deines et al., 1974; Hoefs, 2018; Valley and Cole, 2019). We ascribe this component to biogenic (soil) CO₂ dissolving into groundwaters during their infiltration. Conversely, the samples characterized by more positive δ¹³C_ext values (thermal waters) imply the addition to the shallow CO₂ component with a heavier isotopic signature (δ¹³C of −1‰). This heavy C signature well matches that of deep (crustal/mantle) CO₂ released along the Apennine (Chiodini et al., 2020, 2004). However, also in light of the He isotope evidence previously sjhown (Figure 3), a major mantle C contribution is very unlikely, and we conclude, therefore, that the thermal waters are dominated by a crustal CO₂ component.

In collisional contexts, large crustal CO₂ fluxes (Becker et al., 2008; Gaillardet and Galy, 2008; Perrier et al., 2009; Evans, 2011; Skelton, 2011; Girault et al., 2014; Menzies et al., 2018; Guo et al., 2021) can be sustained by either regional metamorphism (Groppo et al., 2013; Eberhard and Pettke, 2021) or mechano-chemical CO₂ production (Italiano et al., 2009). Metamorphic processes can operate via either 1) decarbonation reactions at relatively high temperatures within calc-silicate rocks (Groppo et al., 2013, 2017, 2020) or 2) dehydration reactions of mineral phases (Eberhard and Pettke, 2021), in which CO₂ degassing is triggered by prograde heating arising from conductive heating triggered by slab breakoff (von Blanckenburg and Davies, 1995), slab rollback (Sizova et al., 2019), or by thermal relaxation of the crust following tectonic thickening upon continent–continent collision. These processes may operate in combination and have certainly interested the past evolution of the Calabrian arc. Notably, metamorphic reactions can have a large CO₂ yield not only where calc-silicate minerals in high-grade rocks and/or limestones are abundant, but also in contexts with relatively few carbonate rocks and/or where siliciclastic metasediments with low carbon contents (<2 wt% C, Pitcairn et al., 2006) prevail. As aforementioned, the entire study area is made up of important metamorphic complexes (Tursi et al., 2021) that include metabasic rocks, felsic granulites, metapelites, and metacarbonate rocks (Schenk, 1984). The conditions for metamorphic CO₂ production are thus certainly met in the study area. Because such metamorphism occurs at very slow rates, amagmatic CO₂ mobilization along convergent plate boundaries can endure over millions of years (Eberhard and Pettke, 2021).

Considering that the sampled springs fall on major active tectonic discontinuities, responsible for the regional crustal seismicity (Rovida et al., 2019; Rovida et al., 2020; Neri et al., 2006, Neri et al., 2020) and for the circulation and discharge of the thermal waters themselves (Italiano et al., 2010; Apollaro et al., 2012; Vespasiano et al., 2012; Vespasiano et al., 2015a; Tiberti et al., 2017; Apollaro et al., 2019b; Apollaro et al., 2020; Vespasiano et al., 2021), it is likewise possible that mechano-chemical CO₂ production is an additional source for crustal CO₂, as already proposed for other active seismic areas as Central Apennines (Italiano et al., 2008), eastern Alps (Italiano et al., 2009), and Japan (Nojima fault; Famin et al., 2008).

5.3 C/³He relationship

During their migration and storage in the crust, fluids can undergo different processes that modify their chemical and isotopic compositions. Insights into these processes, and into volatile sources and sinks, can be derived from a joint analysis and interpretation of He and C isotopic signatures (Holland and Gilfillan, 2013; Barry et al., 2020; Randazzo et al., 2021). In order to reconstruct the original signature of deeply sourced fluids, the samples that are dominated by the atmosphere-sourced volatiles (e.g., the cold water) are initially filtered out. Biogenic carbon can be derived from soil or from a deep source as thermal decarboxylation and pyrolysis of organic matter into metapeliti (δ¹³C from −30 to −20‰, Evans et al., 2008). In light of the He isotope evidence, it is reasonable to think that the carbon present in cold waters come from shallow environments (i.e., soil). In fact, carbon from thermal decarboxylation and pyrolysis should be linked to high He concentrations and crustal isotopic signature, while the cold waters are characterized by atmosphere-derived He. Then, we analyze the remaining sample in a δ¹³C_ext vs. C_ext/³He ratio space (Figure 6), in which the potential C–He sources typically plot in distinct compositional fields: mantle and two crustal sources (biogenic vs. crustal-metamorphic) (modified by Sano and Marty, 1995). We find that our samples have C_ext/³He ratios of 4.4×10⁹ to 9.1×10¹¹ that, coupled with the δ¹³C_ext values, would be consistent with a mixing between crustal-metamorphic (C_ext/³He= 1×10¹³ and δ¹³Cext=−3‰ to +3‰; Evans et al., 2008; Becker et al., 2008; Dai et al., 1996; Sano and Marty, 1995) and mantle (CO₂/³He= 2–4x10⁹ and δ¹³C=−4‰; Marty et al., 2020) fluids in proportions of 88% and 12% (average values; biogenic component would account for 0.8%). However, mantle component fractions (up to 91%) are much higher than calculated above from helium isotopes (2%–3%, Figure 4). This discrepancy can be reconciled taking into account the impact of secondary processes on both C_ext/³He ratios and δ¹³C_ext.

FIGURE 6

FIGURE 6. Plot of δ¹³C_ext versus C_ext/³He. Changes in δ¹³C are calculated following the method from Gilfillan et al. (2009) using the Rayleigh fractionation equation either for precipitation or for dissolution. In the case of precipitation, there is zero ³He loss from the CO₂ phase and CO₂/³He changes in proportion to the fraction of the remaining CO₂ phase, while for CO₂ dissolution, the change in the CO₂/³He ratio is calculated following the Rayleigh equation. The gradual loss of CO₂, with a decrease in the C_ext/³He ratio and the δ¹³C according with Rayleigh-type gas dissolution at different pHs, is shown from broken lines and slim solid lines, while the predicted trend for carbonate mineral precipitation is from the black solid line. Deep end member with C_ext/³He=1×10¹³ (crustal range; Sano and Marty, 1995; O'Nions and Oxburgh, 1988) and δ¹³C_ext from −3 to 3‰ (the mean value of δ¹³C for metamorphic CO₂ is 0.3‰; Dai et al., 1996; Hunt et al., 1996; Clark and Fritz, 1997; Evans et al., 2008). The computed model fit nicely the entire dataset with the samples most affected by the secondary processes that also have the highest He concentrations and the lowest values of R/Ra.

5.4 Secondary processes

CO₂ and He have contrasting solubility in water (Ellis and Golding, 1963; Vogel et al., 1970). As such, the two elements undergo selective gas/water partitioning as deeply rising fluids interact with aquifer(s), ultimately altering the C_ext/³He ratio. Carbon isotopes are likewise fractionated during gas–water interactions (Randazzo et al., 2021 and references therein). Thus, both elemental ratios and δ¹³C can be severely modified by secondary processes such as gas dissolution in water and solid-phase (carbonate) precipitation (Gilfillan et al., 2009; Barry and Bekaer, 2021; Randazzo et al., 2021).

We investigated the possible role of secondary processes (e.g., partial gas dissolution in water and calcite precipitation) during fluid transfer trough the crust by modeling (see Gillfillan et al., 2009) their impact on C_ext/³He ratios and δ¹³C_ext (Figure 6). The process can be modeled as 1) an open-system degassing (Rayleigh-type) at isotopic equilibrium (between phases) and 2) calcite precipitation (Gilfillan et al., 2009). We want to clarify that for a thick crustal sector as Calabrian orogen (thickness up to 38 km; Di Stefano et al., 2009), our model is evidently a simplified approach.

He isotopes indicate a negligible mantle component in thermal waters (up to 2%); hence, the mantle component computed by the C–He relationships (up to 91%) in a simple approach based on mixing between the mantle and crustal end members could be an artefact. Assuming a crustal-metamorphic deep end-member (C_ext/³He of 1×10¹³; Sano and Marty, 1995; O'Nions and Oxburgh, 1988 and δ¹³C_ext of 0.3‰, i.e., the mean value for metamorphic CO₂; Becker et al., 2008; Evans et al., 2008; Dai et al., 1996; Figure 6) as pristine gas composition, in order to explain the variability of the C_ext/³He ratio and δ¹³C_ext in the samples, we used a two-step model: 1) the partial dissolution of He and CO₂ in groundwater and the progressive variation of the C_ext/³He ratio and δ¹³C_ext in the residual gas (Figure 6) and 2) the total dissolution of the residual gas (step 1) into a shallow groundwater. The computed model curves show increasing extents of gas dissolution, over a range of pH values at a fixed temperature (30°C mean sample temperature). It is noteworthy that even using different temperatures (from 10 to 40°C), the two models do not show significant differences.

The results of the modeling fit well with the C_ext/³He ratio and δ¹³C_ext of the thermal waters, suggesting that processes of partial gas dissolution occur at depth. Despite this, we cannot exclude that the lowest δ¹³C values are not at least partially reflecting a biogenic origin (e.g., sample GdN), and carbonate precipitation (together with CO₂ dissolution at a lower pH than 5.7–7; Gilfillan et al., 2009) has not taken a role. In light of this, it is plausible to think that the amount of carbon (i.e., CO₂) present below the study area is much more than what we can measure on the surface and that it can be distributed along a multilayer aquifer of which we can only sample the final member. Our model additionally highlights the role played by the gas–water interaction in determining the composition of fluid released in the studied area and also identifies a metamorphic process as a potential source of CO₂.

5.5 Carbon fluxes

Here, we use the estimated external carbon contributions (C_ext) from thermal and cold waters (Figure 5) to constrain the external carbon outflow through the investigated Calabrian groundwater systems. For each site, we combined the spring flow rates with the computed C_ext. The so-calculated total C_ext flux is 3.63×10⁸ mol year⁻¹ of which 4.19×10⁷ mol year⁻¹ and 3.21×10⁸ mol year⁻¹ from thermal and cold waters, respectively.

For the cold water, we assume that all the C_ext is due only to carbon from shallow biogenic sources and atmospheric CO₂ (C_inf; i.e., C_ext = C_inf). This assumption based on the δ¹³C_ext values is also supported by the presence of atmospheric He dissolved in the cold water, in contrast with the thermal waters that are crustal-He rich (Figures 3, 5). For the thermal waters, the C_inf was deconvolved from total C_ext (see Supplementary Appendix A). The average C_inf of both cold and thermal waters results is 2.5 ± 1 mmol/L, a value very similar to that estimated for recharge waters of the central-southern Apennines (2.31 ± 0.61 mmol/L; Chiodini et al., 2004; Frondini et al., 2019; 2020).

As the deeply derived dissolved carbon (C_deep) is given as C_ext − C_inf, the deep carbon budget is obtained for each spring by multiplying the dissolved C_deep content by the relative water flow rate (see Supplementary Appendix A). The total deep-derived CO₂ output associated with the investigated waters is 2.6×10⁷ mol year⁻¹. However, this value is not representative of the flux from the entire Calabrian orogen and is closely related to the flow rates of the investigated springs as shown by Figure 7.

FIGURE 7

FIGURE 7. CO₂ budget vs. total flow rate for different volcanic and non-volcanic aquifers. The CO₂ budgets are closely related to the flow rates of the investigated springs. The values for CO₂ lost for secondary processes are not included in these budgets. Blue diamonds= non-volcanic area; red diamonds= volcanic area. Data from Caracausi et al. (2015), Becker et al. (2008), Evans et al. (2008), Newell et al. (2008), and Chiodini et al. (2004).

Therefore, in order to better compare the contribution of deep carbon released from different areas, we also computed the specific flux for each spring, that is, the deep carbon outputs normalized for the catchment areas (see Supplementary Appendix A).

We estimate a value of 2.6x10⁷ mol year⁻¹ of deep carbon for the thermal water. Including the percentage of carbon lost due to secondary processes (Section 5.4), we estimate that a value of 6.1×10⁸ mol year⁻¹ of deeply derived carbon for an area of 2,880 km² (Supplementary Table S2A) produces a specific flux of 2.1×10⁵ mol km⁻² year⁻¹.

This value ranges from two orders of magnitude lower to the same order of CO₂ fluxes defined for some volcanic aquifers (dissolved CO₂ in Figure 8; e.g., Etna 5.7×10⁶ mol km⁻²year⁻¹, Vulcano 3.4×10⁶ mol km⁻² year⁻¹, Mammoth Mountain 1.7×10⁷ mol km⁻² year⁻¹, Mt. Amiata 1.1×10⁶ mol km⁻² year⁻¹, Roccamonfina 4.8×10⁵ mol km⁻² year⁻¹; Caracausi et al., 2015 and reference therein). Compared with the values from the active tectonic region and collisional orogen (violet gradient background in Figure 8), our result is one order of magnitude lower than that of the deep CO₂ flux estimated for the central Appenine (2.4×10⁶ mol km⁻² year⁻¹; Chiodini et al., 2000) and well fitted the range of values estimated for some Himalayan areas. In the latter case, our values are higher than those calculated for the southern margin of the Tibetan Plateau (9.9×10³ mol km⁻² year⁻¹; Newell et al., 2008), of the same order for the Narayani Basin (4.4×10⁵ mol km⁻² year⁻¹; Evans et al., 2008) and one order of magnitude lower than flux values estimated by Becker et al. (2008) for the Marsyandi Valley area (1.1×10⁶ mol km⁻² year⁻¹). For the Marsyandi Valley and Narayani basin, the values estimated included the CO₂ lost for degassing processes, but no correction related to the carbon content link to carbonate dissolution or to the biogenic source was made. Therefore, it is probable that the values related to the deep CO₂ are slightly lower. In any case, the values calculated for the Calabria region are comparable to those of areas in which a strong outgassing of CO₂ from a metamorphic source has been identified. We want to clarify that our estimate has many assumptions, and thus an unknown uncertainty, but our intent here is to present a possible limit on CO₂ fluxes for the Calabria region. However, if accurate, our values suggest that the Calabrian orogen is an important contributor to the global carbon budget today.

FIGURE 8

FIGURE 8. Deep CO₂ flux from volcanic and non-volcanic areas. Calabria Cor includes the percentage of deep CO₂ lost for secondary processes (Section 5.4). The Marsyandi Valley and Narayani Basin values refer to total C without extrapolation from biogenic and deep components and included CO₂ lost for degassing processes (Becker et al., 2008; Evans et al., 2008). Orange-yellow background for volcanic areas; purple background for non-volcanic areas. Southern Tibetan Plateau from Newell et al. (2008). Modified after Caracausi et al. (2015).

6 Conclusion

Understanding carbon isotopic signatures and the processes that affect them are critical to make accurate CO₂ flux estimates and identify the origin of carbon. In this study, we investigate the sources and sinks of fluid dissolved in groundwaters from the Pollino and Calabria regions. Chemical and isotopic compositions (He and C) of sampled waters allowed us to identify two different domains: 1) a shallow system dominated by gas components of the atmospheric signature (He) and biogenic origin (C), and 2) a deeper system in which crustal/deep fluids (CO₂ and He) are dominant.

The external carbon contributions have been calculated following the mass balance approach and coupled with the helium data, which allowed us to identify a deep CO₂ (i.e., crustal/metamorphic) associated with the fluids released in the hydrothermal basins and to detect the secondary process (dissolution/precipitation), which modify the pristine chemical and isotopic compositions of fluids, affecting the deep carbon budget. The samples with the highest He concentrations and the lowest ³He/⁴He ratios also are the most affected by carbon removal processes. This could indicate longer residence time in the crust and/or a more complex circulation system (multilayer) of which the sampled waters are only the last member. We also proposed metamorphic processes as a source of CO₂, but on this, way more studies are required.

For the investigated springs, a total deep-derived CO₂ output of 2.6×10⁷ mol year⁻¹ was computed, but the value strongly depends on the flow rates and does not represent the flux from the entire Calabrian orogen. Extrapolating from our model, for the percentage of CO₂ lost due to secondary processes, we estimated a maximum value of 6.1×10⁸ mol year⁻¹ for the deeply derived carbon. Scaling our estimate of deep CO₂ flux to the whole study area (2,880 km²), a value of 2.1×10⁵ mol km⁻² year⁻¹ is obtained. The value compares to other globally significant carbon fluxes as that defined for the Himalaya orogen and Central Appenine in Italy and despite the considerable uncertainties, represents the first estimate of CO₂ flux for the Calabria region based on the field sampling and modeling of secondary processes facilitating a comparison with estimates from other collisional orogens. Given the extent of flux values extrapolated, we emphasize that more studies should be conducted to implement knowledge on possible sources, circulation systems, and deep CO₂ release by defining the contribution that a collision area such as the Calabrian arc can provide to the global carbon budget.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Author contributions

PR collected all samples, carried out the analyses, processed the data, and wrote the manuscript. AC helped in sampling and analyses supervision. CA, GV, and MP helped in sampling and analyses. AR helped in sampling. AA and CC helped in organization and funds. All the authors contributed to the discussion and co-wrote/corrected the manuscript.

Funding

This work was supported by the MIUR project PRIN 2017-2017LMNLAW “Connect4Carbon” and DCO Grant 10881-TDB “Improving the estimation of tectonic carbon flux.”

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2022.946707/full#supplementary-material

References

Acquafredda, P., Lorenzoni, S., and Zanettin Lorenzoni, E. (1994). Palaeozoic sequences and evolution of the calabrian-peloritan arc (southern Italy). Terra nova. 6, 582–594. doi:10.1111/j.1365-3121.1994.tb00525.x