Tectonically active regions are often characterized by large amounts of carbon dioxide degassing, and estimation of the total CO2 discharged to the atmosphere from tectonic structures, hydrothermal systems and inactive volcanic areas is crucial for the definition of present-day global Earth degassing. The carbon balance of regional aquifers is a powerful tool to quantify the diffuse degassing of deep inorganic carbon sources because the method integrates the CO2 flux over large areas. Its application to peninsular Italy shows that the region is characterized by specific CO2 fluxes higher than the baseline determined for the geothermal regions of the world, and that the amount of endogenous CO2 discharged through diffuse regional degassing (c. 2.1 × 1011 mol a−1) is the major component of the geological CO2 budget of Italy, definitely prevailing over the CO2 discharged by Italian active volcanoes and volcanoes with hydrothermal activity. Furthermore, the positive correlation between geothermal heat and deep CO2 dissolved in the groundwater of central Italy suggests that (1) the geothermal heat is transported into the aquifers by the same hot CO2-rich fluids causing the Italian CO2 anomaly and (2) the advective heat flow is the dominant form of heat transfer of the region.

Supplementary material: The location, flow rate, extent of the hydrogeological basin, chemical and isotopic analyses of the 160 springs considered in this study, and the results of the carbon mass balance are reported in a table available at https://doi.org/10.6084/m9.figshare.c.4237025

On geological timescales, the level of atmospheric CO2 is controlled by the balance between CO2 consumed by chemical weathering and CO2 degassed by metamorphism and magmatism (Walker et al. 1981; Berner et al. 1983; Berner 1993, 2004; Kerrick & Caldeira 1993). The BLAG (Berner–LAsaga–Garrels; Berner et al. 1983) and subsequent models for the long-term global carbon cycle (e.g. Berner 2006) derive the CO2 degassing rate over geological times assuming that present-day CO2 degassing equals the flux of CO2 consumed by chemical weathering. However, the global estimates of present-day CO2 degassing flux based on this assumption are inconsistent with those based on volcanic degassing data (Burton et al. 2013, and references therein), suggesting that the fluxes of CO2 consumed by chemical weathering and those released by metamorphic–magmatic degassing should be computed separately (Fresia & Frezzotti 2015).

The number of CO2 flux measurement studies at volcanoes is continuously increasing but the CO2 emissions from the majority of volcanic sources are still unknown. Extrapolating the available data to all the active volcanoes, Burton et al. (2013) estimated a global volcanic subaerial CO2 flux of about 12 × 1012 mol a−1 and a total flux, including mid-ocean ridge degassing, of 14 × 1012 mol a−1. These values are higher than previous estimates (e.g. Marty & Tolstikhin 1998; Mörner & Etiope 2002) but are still affected by large uncertainties mostly related to the estimation of diffuse degassing from inactive volcanoes and hydrothermal structures.

The impact of metamorphic degassing on the global carbon cycle is even more uncertain (Gaillardet & Galy 2008). All the estimations derive from theoretical models and studies on the Himalayan orogen (Kerrick & Caldeira 1998; Mörner & Etiope 2002; Gorman & Kerrick 2006; Rapa et al. 2017) but there are only few estimations of present-day metamorphic CO2 degassing based on direct measurements (e.g. Evans et al. 2008; Girault et al. 2016).

The so-called non-volcanic diffuse CO2 degassing (i.e. not directly connected to active volcanoes) represents an important aspect of the global carbon cycle and its quantification is crucial to give a realistic estimate of the degassing process at a global scale (Kerrick et al. 1995; Chiodini et al. 1999, 2000, 2004b; Kerrick 2001; Lee et al. 2016; Hunt et al. 2017).

Tectonically active regions, including orogenic belts, continental rifts and, more in general, hydrothermal or inactive volcanic areas, are often characterized by widespread regional carbon dioxide degassing. In these areas carbon dioxide, derived from mantle degassing and/or metamorphism, is continuously released to the atmosphere through various processes, including emissions from focused gas vents, diffuse soil degassing, and degassing from lakes and hot and cold springs. The coincidence of carbon dioxide discharges and major zones of seismicity was first recognized by the pioneering studies of Barnes et al. (1978) and Irwin & Barnes (1980), and successive studies (Kerrick & Caldeira 1998) clearly showed that large orogenic belts, in addition to their role as an atmospheric CO2 sink through silicate weathering (Gaillardet et al. 1999), also aid in the production of CO2-rich fluids, mainly related to regional metamorphism. More recently, a very large amount of mantle CO2 is hypothesized to be released by soil diffuse degassing in the East African Rift (Lee et al. 2016; Hunt et al. 2017).

A variety of techniques have been used to estimate present-day CO2 regional degassing. In the geothermal areas of the Taupo Volcanic Zone, New Zealand, and the Salton Trough region, California, USA, the convective hydrothermal CO2 emission was computed using data on convective heat flow, temperature and CO2 concentration of reservoir fluids (Kerrick et al. 1995; Seward & Kerrick 1996). The results showed that the two areas were characterized by a similar specific CO2 flux (c. 106 mol a−1 km−2), and according to Kerrick et al. (1995) this value could be used as a baseline to compute convective hydrothermal CO2 emission from other areas of high heat flow. However, using a similar approach, Frondini et al. (2008) computed a specific CO2 flux about five times higher (4.8 × 106 mol a−1 km−2) for Tuscany and Northern Latium (Italy).

In the areas where a flux of deep CO2 directly affects regional aquifers, most of the deeply generated gas can be dissolved by groundwater because the relatively high solubility of CO2 in water (Rose & Davisson 1996; Sorey et al. 1998; Chiodini et al. 1999; Caliro et al. 2005). Hence, studying the concentration of total dissolved inorganic carbon (TDIC) in groundwater is a suitable way to infer the regional flux of endogenous CO2.

The Earth degassing process was extensively studied in Italy and numerous techniques currently used for the estimation of non-volcanic CO2 degassing at regional scale were originally developed for the study of the Italian territory.

At Colli Albani, central Italy (Chiodini & Frondini 2001), the CO2 emission was computed by summing the contributions of CO2 transported by groundwater and the CO2 released by two large zones of diffuse soil degassing. The flux of the deeply derived CO2 dissolved by groundwater was derived by the carbon balance of the volcanic aquifer, whereas that discharged from the areas of diffuse soil degassing was measured using the accumulation chamber method (Chiodini et al. 1998). The resulting area-averaged deep CO2 production rate is 2.8 × 106 mol a−1 km−2, of the same order of magnitude as the CO2 fluxes computed by Frondini et al. (2004) for Tuscany and Northern Latium. Similar values were also estimated by Gambardella et al. (2004) for the volcanic aquifers of central–southern Italy. Chiodini et al. (2000, 2004b), coupled the mass and isotopic balance of TDIC with hydrogeological data from the main springs of the Italian Apennines and integrated the groundwater data with data on vent and diffuse soil degassing, and derived a regional map of CO2 Earth degassing that showed that two large degassing structures affect the Tyrrhenian side of the Italian peninsula and that c. 40% of the inorganic carbon in the groundwater derives from magmatic–mantle sources.

In this paper we discuss some of the main methods used to measure non-volcanic degassing over large areas and their application to the significant case of the Italian Apennines. In the following sections, after a brief description of the geological, hydrogeological and geophysical features of the Apennines, we will discuss (1) the methods used to measure the regional flux of deep CO2 associated with groundwater circulation, (2) the relationships between regional flux of CO2 and areas of focused and diffuse gas emission from soil, as well as the conditions for the development of such zones of discharge, and (3) the relationships between regional CO2 degassing and advective heat flow.

The Apennines are a series of mountain ridges that run the entire length of the Italian Peninsula, extending from Liguria (northern Italy) to Sicily. They form the physical backbone of peninsular Italy and are divided into three sectors: northern, central and southern Apennines. The study area (Fig. 1), delimited to the west by the Tyrrhenian coastline, includes the central Apennines and part of the northern and southern Apennines.

The stratigraphical successions of the Apennines are essentially composed by shallow-water carbonate platform formations and pelagic formations (limestones and marly limestones), followed by foredeep deposits (arenaceous turbidites, marls and calcarenites). The present structural setting of the Apennines is the result of the superimposition, since the Early Miocene, of two concurrent tectonic processes, compression in the foreland and extension in the hinterland. The two processes are co-axial and both are characterized by an eastward migration over time. As the result of this tectonic evolution, the Apennines show two different crustal domains: a western Tyrrhenian domain, where extensional deformation destroyed the pre-existing compressional belt, and an eastern Adriatic domain, where the compressional structures are still preserved (Barchi 2010; Cosentino et al. 2010).

Crustal thinning related to extension in the Tyrrhenian area was followed by intense magmatic activity along the western margin of the Italian peninsula. Igneous rocks, ranging in age from Late Miocene to Recent, are widespread through the Tyrrhenian extensional zone (Barberi et al. 1971; Civetta et al. 1978; Peccerillo 1999; Peccerillo & Frezzotti 2015).

The western Tyrrhenian domain is characterized by a thin crust (0–25 km), high heat flow, locally more than 200 mW m−2 (Della Vedova et al. 1984), positive magnetic anomalies (Arisi Rota & Fichera 1985), shallow earthquakes and positive gravity anomaly. The most important Italian geothermal fields (e.g. Larderello, Amiata, Latera, Torre Alfina and Campi Flegrei) are located in this region (Fig. 1).

The eastern Adriatic domain is characterized by lower heat flow values, negative gravity anomalies and an average crustal thickness of about 35 km. This area is also characterized by high seismicity, which has caused several strong earthquakes in recent times (Amato et al. 1998; Chiarabba et al. 2009; Chiaraluce et al. 2017).

The hydrogeological settings of the western Tyrrhenian and eastern Adriatic domains reflect the differences described above: the western region is characterized by small shallow aquifers hosted by Quaternary volcanic rocks and small isolated carbonate structures deriving from the disruption of the pre-existing compressional structures; in contrast, the eastern region is characterized by large regional aquifers hosted by Mesozoic permeable limestones (e.g. Boni et al. 1986).

Several regional-scale studies showed that the western Tyrrhenian domain is affected by intense CO2 degassing, which results in numerous cold, CO2-rich gas emissions at the surface (e.g. Chiodini et al. 2000, 2004b, 2011; Rogie et al. 2000; Minissale 2004). Different techniques were used to measure the gas flux depending on the specific features of each gas emission. For a detailed description of these methods, the reader is referred to the specific literature (Chiodini et al. 1998, 2010; Rogie et al. 2000; Cardellini et al. 2003; Costa et al. 2008). The measured gas flow rates in some cases are very high, as, for example, at Mefite d'Ansanto, the largest cold CO2 emission ever measured on Earth, which releases amounts of CO2 (c. 1.7 × 1010 mol a−1) approaching those released from the most active volcanoes (Chiodini et al. 2010).

The idea of measuring CO2 flux at regional scale derived from the observation that the main Apennine springs emit large amounts of dissolved CO2. In principle, the spring discharge can be used to estimate the specific CO2 fluxes over large areas; that is, the areas of the hydrogeological basins. However, before computing the CO2 fluxes from TDIC it is essential to define the origin of the dissolved gas in water because in addition to a possible deeply derived component, atmospheric sources, biogenic activity and water–rock interactions contribute to TDIC. This aspect is solved by applying the mass-balance approach described below.


The dataset used in this study (supplementary material) refers to 160 Apennine springs (Fig. 1) that were sampled in different surveys (Civita 1977; Celico et al. 1980; Chiodini et al. 2000, 2004b, 2013; Frondini et al. 2012). Water temperature, pH and total alkalinity were determined directly in the field. The chemical (Ca, Mg, Na, K, Cl, SO4, NO3) and stable isotope (H2O and carbon) composition was determined according to the analytical methods reported in the original studies. The TDIC was computed with the code PHREEQC (Parkhurst & Appelo 1999), using the field determinations of T, pH and alkalinity, and the major ion concentrations as input data. For each sample two additional variables (Cext and δ13Cext; supplementary material) were computed according to the method explained in the next section.

The strategy of the survey was to sample the high-discharge springs to obtain a dataset highly representative of the Apennines groundwater. The flow rates of the springs and the areal coverage for the corresponding hydrogeological basins (supplementary material) were taken from the literature (Celico et al. 1979; Celico & Corniello 1980; Celico 1983; Boni et al. 1986).

Carbon mass balance of the aquifers

To compute the flux of endogenous CO2 transported by groundwater it is necessary to discriminate the different sources of dissolved carbon dioxide. A theoretical treatment of the evolution of the carbon isotopes in natural waters has been given by Wigley et al. (1978). The effects of an arbitrary number of sources (such as dissolution of carbonate minerals, input of deeply derived CO2, oxidation of organic material) and sinks (such as mineral precipitation and CO2 degassing), and equilibrium fractionation between solid, gas and aqueous phases are considered. The results are expressed as equations relating changes in isotopic composition to changes in carbonate chemistry.

The evolution of the TDIC and 13C/12C ratio of natural water systems is described by the following equations:
where Ii and Oi are the molalities of the ith input and the ith output of carbon in the solution, R is the 13C/12C isotopic ratio of the solution as a whole, R* is the 13C/12C isotopic ratio of the ith input species, αis is the fractionation factor between the ith output and the solution, and N and M are the numbers of incoming and outgoing carbon species. For a finite-difference solution, equation (1) becomes
and equation (2), considering the isotopic values in terms of delta units, becomes

In equations (2) and (4) it is assumed that isotopic equilibrium exists between all dissolved carbon species in solution and that all outputs are in isotopic equilibrium with the parent solution. Furthermore, it is assumed that no fractionation occurs during the input of any source to the solution, and fractionation occurs only during the output of carbon from the parent solution.

In the case of the regional Apennines aquifers, equations (3) and (4) were rearranged assuming that no carbon sinks are present in the aquifer, an assumption that was demonstrated to be realistic for most of the studied springs (Chiodini et al. 2000, 2011). In the case of carbonate aquifers, the no-sink assumption means that carbonate precipitation and CO2 degassing do not significantly occur before the emergence of groundwater at springs. In this case the amount of carbon dissolved by groundwater is given by the sum of the carbon deriving from carbonate dissolution and the carbon deriving from sources external to the aquifer. It is possible to eliminate the contribution of carbonate dissolution (Ccarb) to TDIC by computing two new variables, Cext and δ13Cext:
where Cext is given by the sum of the carbon deriving from a deep source (Cdeep) and the carbon content of infiltrating waters (Cinf; i.e. the atmospheric CO2 plus that from biogenic sources active in the soils during the infiltration):
from which derives
where δ13Cinf and δ13Cdeep correspond to the isotopic composition of Cinf and Cdeep respectively.

In the above set of equations, only the TDIC and δ13CTDIC are analytically determined whereas the other variables are derived from general considerations and from the geochemical interpretation of the data.

For example, in the carbonate aquifers of the Apennines, Ccarb is given by
Equation (9) states that the moles of Ccarb are equivalent to the moles of Ca and Mg entering the solution from carbonate mineral dissolution, and the negative term –mSO4 corrects for the moles of Ca in solution that were derived from gypsum or anhydrite dissolution. In the case of aquifers composed of silicate crystalline rocks or volcanic rocks, the contribution of carbonate dissolution to TDIC may be assumed negligible and Cext is equal to TDIC.

δ13Ccarb is the carbon isotopic composition of the carbonate rocks and preferably should be determined based on known values for the carbonate rocks that host the aquifers (e.g. in the Apennine aquifers it was assumed to be +2.21‰, based on hundreds of isotopic analyses of carbonate rocks; Chiodini et al. 2004b).

The values of Ccarb and δ13Ccarb allow one to compute Cext and δ13Cext (from equations (5) and (6)). The other unknown variables (Cinf, Cdeep, δ13Cinf and δ13Cdeep) are computed from equations (7) and (8), by making suitable assumptions constrained by the geochemical interpretation of the data.

Once all the terms of the carbon budget are determined (Ci), the total output of CO2(QCO2,i) and the corresponding CO2 flux (FCO2,i) can be computed for any carbon source using the following equations:
where q is the spring flow rate and A is the surface area of the hydrogeological basin of the spring.

The case of the Apennine aquifers is described below.

Carbon mass balance of Apennines aquifers

The computation of Cdeep is based on the analysis of the Cext v. δ13Cext diagram (Fig. 2) where data from the Apennines springs are distributed along mixing lines consistent with expected theoretical behaviour and calculated using equations (7) and (8). A cluster of samples (blue circles in Fig. 2; described as blue samples below) is characterized by low Cext values (<3–4 mmol l−1) and very negative δ13Cext (from −16 to −25‰). These are the springs where Cdeep is absent (or negligible) and where Cext practically coincides with Cinf. Most of the other samples (red circles in Fig. 2) show an enrichment in Cext and higher δ13Cext values, indicative of the addition of isotopically heavier carbon, typical of endogenous sources (i.e. the deeply derived Cdeep). The Cinf values of the blue samples range in a narrow interval and the mean value (2.31 ± 0.61 mmol l−1) can be assumed as representative of the recharge waters of the Apennines. Using this value in equation (7) we can compute Cdeep for each spring affected by the input of deep CO2 (red circles in Fig. 2). In this computation, for the few samples where Cext is lower than the mean Cinf value, we assumed Cdeep = 0.

Alternatively, Cdeep can be computed assuming (or evaluating) the isotopic compositions of the pure components. This method, described by Chiodini et al. (2011), is generally applicable in areas where one can assume the deep carbon source has a unique isotopic signature.

To characterize the isotopic signature of the deep carbon sources in the Apennines, we have to consider that the no-sink assumption, which is necessary for application of this method, is realistic only for samples characterized by TDIC <40 mmol l−1 (which roughly corresponds to a Cext of c. 30 mmol l−1; dashed line in Fig. 2). This threshold was defined by Chiodini et al. (2011) modelling the effects of CO2 input on the dissolved carbon species of water at equilibrium with calcite. Above this threshold samples could be affected by CO2 degassing and calcite precipitation, which can cause significant isotopic fractionation and a decrease of Cext (i.e. it would not be representative of mixing between Cinf and Cdeep). Figure 2 shows that samples with Cext lower than the no-sink threshold plot in the theoretical mixing field computed by adding to the infiltrating water (Cinf=2.31 ± 0.61 mmol l−1; δ13Cinf = −21.6 ± 2.9‰) variable amounts of Cdeep with δ13Cdeep ranging from −5 to +1‰, the typical range of deeply derived CO2 emitted in Italy by volcanoes, geothermal systems and cold CO2 emissions (Chiodini et al. 2004b).

The estimated concentrations of the different carbon sources multiplied by the flow rate of the springs (equation (10)) give the total output of CO2 associated with the Apennines groundwater as follows: QCO2,carb=2.2×1010mola1, QCO2,inf=1.6×1010mola1 and QCO2,deep=2.9×1010mola1. It is noteworthy that the results show that the deep source of CO2 is the main source supplying inorganic carbon to the Apennine groundwaters.

Regional map of CO2 Earth degassing

To derive the map of the regional CO2 flux, following the approach of Chiodini et al. (2004b), the dataset of the Apennine springs was integrated with literature data for springs from smaller carbonate aquifers in the western part of Italy (with flow rates >10 l s−1). Furthermore, we also considered in the computation the average CO2 flux affecting the volcanic aquifer of Mt. Albani, as previously estimated by Chiodini & Frondini (2001). Starting from the Cext values derived from equation (5), the CO2 flux of each spring (in mol a−1 km2) was computed by multiplying its Cext by the corresponding flow rate and dividing the result by the surface area of its hydrogeological basin (equations (10) and (11)). The CO2 fluxes computed for each spring were subsequently determined using a geostatistical approach based on the sequential Gaussian simulation procedure (sGs; Deutsch & Journel 1998), which frequently is applied to derive maps of soil CO2 diffuse degassing from field measurements (e.g. Cardellini et al. 2003). In this study, each spring was considered as a single CO2 flux measurement point for the application of the sGs algorithm and the derivation of the regional map. We computed 200 maps of CO2 flux by sGs according to the geostatistical parameters of the dataset (a spherical variogram model with nugget = 0.35, sill = 1 and range = 50 km), weighting each datum point for the spring flow rate. The average values computed from the 200 simulations at each location are shown in the map of Figure 3. The map highlights the presence of two large regional anomalies of CO2 flux on the western side of the Italian peninsula: the Tuscan Roman Degassing Structure and the Campanian Degassing Structure (TRDS and CDS in Fig. 3; Chiodini et al. 2004b). On the map the areas degassing deeply derived CO2 are shown by colours from light green to red, corresponding to Cext flux values higher than 3 × 106 mol a−1 km−2. This threshold was computed considering a reasonable maximum value for Cinf of biological origin (c. 0.004 mol l−1, Chiodini et al. 2004b) and using equations (10) and (11) to convert concentration into flux.

The northern degassing structure (TRDS) partially overlaps the Tuscany, Roman magmatic provinces whereas the southern structure (CDS) is related to areas of Campanian volcanism. The two volcanic provinces are characterized by Quaternary potassic and ultrapotassic magmas rich in fluids with high CO2/H2O ratios (Foley 1992). Magma geochemistry is consistent with the melting of a mantle source metasomatized by the addition of subducted crustal material (e.g. Peccerillo 1999; Frezzotti et al. 2009). In agreement with Chiodini et al. (2004b, 2011) we suggest that the TRDS and CDS mostly reflect degassing of this metasomatized mantle.

It should be noted that the real existence of the TRDS and CDS is confirmed by the presence of numerous (c. 140) widespread CO2-rich gas emissions consisting of cold gas vents, areas of soil diffuse degassing and bubbling waters (Fig. 3; www.magadb.net). In particular, the free gas emissions occur in the western sectors of the TRDS and CDS, which are characterized by outcrops of relatively low-permeability formations, whereas in the eastern sectors, where permeable carbonate formations prevail, the deeply derived gas is almost completely dissolved by groundwater circulating in the Apennine aquifers.

The total amount of deeply derived carbon released from the TRDS and CDS was estimated by integrating the simulated CO2 flux values (i.e. Cext) over the area and subtracting the contribution from the atmospheric and organic sources (i.e. Cinf). This computation assumes a mean recharge water infiltration of 20 l s−1 km−2, and gives a deep carbon flux of 1.4 × 1011 mol a−1 for the TRDS and 0.7 × 1011 mol a−1 for the CDS, corresponding to deeply derived specific CO2 fluxes of c. 3.7 × 106 mol a−1 km−2 and c. 4.7 × 106 mol a−1 km−2, respectively. The total amount of deeply derived CO2 discharged by this sector of Italy (c. 2.1 × 1011 mol a−1) corresponds to about 2–15% of the estimated present-day global CO2 fluxes from active volcanoes (Burton et al. 2013), highlighting the global relevance of the regional CO2 Earth degassing process in the region. Finally, the computed specific CO2 fluxes are 4–5 times higher than the value of c. 106 mol a−1 km−2 suggested as the reference value for the convective hydrothermal CO2 emission calculated for high heat flow areas (Kerrick et al. 1995).

Advective heat transport and CO2 fluxes

On the heat flux map of Italy the Apennine portions of the TRDS and CDS correspond to zones of apparently very low heat flux (Cataldi et al. 1995). However, this map refers only to the conductive heat flux because the map is based on temperature gradients measured in wells. In regions characterized by permeable rocks that host large aquifers, such as the Apennines, groundwater circulation invalidates use of a purely conductive model to estimate heat flux. Instead, because of the abundant groundwater circulation, advective heat flow can be the dominant form of heat transfer, and the temperature of spring water can be used to estimate the values of geothermal heat flux more realistically with respect to the estimations made using the thermal gradients measured in deep wells. A method based on the energy balance of groundwater was applied to the US Cascade Range, allowing the identification of a deep thermal source active in the area (e.g. Ingebritsen et al. 1989; Manga 1998; Ingebritsen & Mariner 2010). More recently, Chiodini et al. (2013) used the same method to estimate the advective geothermal heat flux (together with the CO2 flux) for 11 aquifers of the central Apennines belonging to the TRDS (Fig. 4).

Briefly, the geothermal heat flux (FH) is computed starting from ΔT, the temperature difference between the recharge water and the water discharged from the springs. In particular, the method considers the effect on water temperature for both geothermal heating and dissipation of gravitational potential energy (Manga & Kirchner 2004). Assuming that the effect of the conductive heat transfer from the aquifer to the surface is negligible, a reasonable assumption for aquifers with a water table at depths higher than 100 m, the temperature difference between the discharge (Ts) and the infiltration water (Tr), ΔT (in K), is given by
where FH is the geothermal heat flux (in W m−2), ρw (kg m−3) and Cw (J kg−1 K−1) are the density and the heat capacity of water respectively, A is the areal extent of the hydrogeological basin covered by each spring (m2), q is the spring flow rate (m3 s−1), Δz (m) is the difference in elevation between the water recharge area (Zr) and the spring (Zs), and g is the gravitational acceleration (m s−2). In equation (12), the first term relates to the geothermal heating and the second relates to the dissipation of gravitational potential energy.
The total amount of geothermal heat released by an aquifer is given by

In Table 1 we show the total geothermal heat (QH) and the CO2

(QCO2,deep) outputs computed for 11 aquifers of the TRDS by Chiodini et al. (2013) and for one aquifer of the CDS (Matese aquifer; Di Luccio et al. 2018). For further details on the method and on the computation of the variables in equations (12) and (13) we refer the reader to the original studies.

The QH and QCO2,deep values estimated for the 12 aquifers show a positive correlation (Fig. 5), suggesting that the gas and the heat are advectively transported into the aquifers by hot CO2-rich fluids. It is noteworthy that the original fluids have a QH/QCO2 ratio similar to that of the hot liquids circulating in the nearby geothermal systems of Latium (Latera and Torre Alfina; Figs 4 and 5). The total CO2 output from the two geothermal systems (Table 1) was estimated through hundreds of CO2 flux measurements performed with the accumulation chambers technique (Chiodini et al. 2007; Lucidi 2010). Considering the temperature and the CO2 contents of the geothermal liquids at Torre Alfina (120°C, 0.37 mol kg−1; Gambardella et al. 2004) and Latera (212°C, 0.72 mol kg−1; Chiodini et al. 2007) we estimated QH of the two systems, computing the total amount of liquid necessary to supply the measured QCO2 and multiplying it by the specific enthalpy of the liquid at the measured temperature.

For comparison, in Figure 5 are also reported the data from the geothermal Taupo Volcanic Zone (TVZ) that were used by Kerrick et al. (1995) to derive the specific CO2 flux value of c. 106 mol a−1 km−2 assumed by those researchers as the baseline for the convective hot zones of the Earth. Figure 5 explains well why in central Italy the specific CO2 fluxes are 4–5 times higher than those derived from TVZ data: the central Italy original fluids have a QH/QCO2 ratio that is much lower than that in the TVZ.

A main goal of modern geoscience is to better define and quantify the geological components of the global carbon cycle. In particular, the amount of CO2 derived from inorganic sources (mantle, crust) and released to the atmosphere is a central parameter that at the moment is known only with large uncertainties. A special study has been made recently to quantify the CO2 emitted by the plumes of volcanoes (DCO DECADE initiative, https://deepcarboncycle.org/about-decade). However, no effort has been devoted to understanding and quantifying diffuse degassing, whose contribution to the total CO2 flux to the atmosphere is practically unknown, despite some appreciable but local research.

Here we have shown that knowledge of the carbon balance of an aquifer is a powerful tool to quantify the diffuse degassing of deep inorganic carbon sources. The results of the carbon balance of aquifers can give the average CO2 flux over large areas (typically of tens to hundreds of square kilometres in the Apennines).

The application of the method in Italy shows that the specific CO2 fluxes (i.e. c. 3.7 × 106 mol a−1 in the TRDS and c. 4.7 × 106 mol a−1 in the CDS) are higher than the reference baseline specific CO2 flux for hot, hydrothermal zones of the Earth (106 mol a−1; Kerrick et al. 1995). This could reflect a real strong CO2 anomaly in Italy and/or a difference in the method of computation. It should be noted that the method based on the carbon and heat balance of the aquifers (Table 1, Fig. 5) refers to a natural steady-state condition of the process of ascent of hot, CO2-rich fluids and their consequent mixing with shallow groundwater. On the other hand, the method based on data from deep geothermal wells (i.e. Kerrick et al. 1995) could be affected by the perturbation of the steady-state condition caused by the extraction of the geothermal fluids. It should be noted that the two geothermal systems of central Italy (Latera and Torre Alfina), which give results comparable with those returned by the aquifers, were exploited only for short periods and are currently unexploited.

Finally, the total amount of deeply derived CO2 (c. 2.1 × 1011 mol a−1) released by the TRDS and CDS, the two degassing structures located in peninsular Italy (Fig. 3), is the main component (55%) of the geological CO2 budget of Italy (Fig. 6) including the SO2-bearing plumes of active volcanoes (Etna, Stromboli, Vulcano; Burton et al. 2013) and the volcanoes with hydrothermal activity where SO2 is practically absent (Campi Flegrei and Vesuvio; Frondini et al. 2004; Cardellini et al. 2017). The inclusion in the budget of other hydrothermal systems, such as Ischia and Pantelleria, where CO2 flux estimates are minor with respect to the total budget (and may also be affected by large uncertainties; Favara et al. 2001; Chiodini et al. 2004a; Pecoraino et al. 2005; Granieri et al. 2014), would not change this picture of the diffuse degassing as the main natural producer of inorganic geological CO2 in Italy.

The results obtained in Italy demonstrate that total CO2 flux estimates cannot be reliably quantified without the investigation of groundwaters, which in permeable orogens of tectonically young and active areas can dissolve most, if not all, the CO2 rising from depth. Although it has long been recognized (Barnes et al. 1978; Irwin & Barnes 1980) that seismically active regions worldwide are characterized by the occurrence of CO2 degassing, quantitative data on CO2 fluxes are practically missing for most of these. We believe that investigation of diffuse degassing in these regions is crucial to better constrain the global carbon flux.

This work was funded by the Istituto Nazionale di Geofisica e Vulcanologia.

Scientific editing by Igor Maria Villa

This is an Open Access article distributed under the terms of the Creative Commons Attribution 3.0 License (http://creativecommons.org/licenses/by/3.0/)