Measuring and interpreting CO2 fluxes at regional scale: the case of the Apennines, Italy

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 CO 2 is controlled by the balance between CO 2 consumed by chemical weathering and CO 2 degassed by metamorphism and magmatism (Walker et al. 1981;Berner et al. 1983;Berner 1993Berner , 2004Kerrick & 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 CO 2 degassing rate over geological times assuming that present-day CO 2 degassing equals the flux of CO 2 consumed by chemical weathering. However, the global estimates of present-day CO 2 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 CO 2 consumed by chemical weathering and those released by metamorphic-magmatic degassing should be computed separately (Fresia & Frezzotti 2015).
The number of CO 2 flux measurement studies at volcanoes is continuously increasing but the CO 2 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 CO 2 flux of about 12 × 10 12 mol a −1 and a total flux, including mid-ocean ridge degassing, of 14 × 10 12 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 CO 2 degassing based on direct measurements (e.g. Evans et al. 2008;Girault et al. 2016).
The so-called non-volcanic diffuse CO 2 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. 1999Chiodini et al. , 2000Chiodini et al. , 2004bKerrick 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 CO 2 sink through silicate weathering (Gaillardet et al. 1999), also aid in the production of CO 2 -rich fluids, mainly related to regional metamorphism. More recently, a very large amount of mantle CO 2 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 CO 2 regional degassing. In the geothermal areas of the Taupo Volcanic Zone, New Zealand, and the Salton Trough region, California, USA, the convective hydrothermal CO 2 emission was computed using data on convective heat flow, temperature and CO 2 concentration of reservoir fluids (Kerrick et al. 1995;Seward & Kerrick 1996). The results showed that the two areas were characterized by a similar specific CO 2 flux (c. 10 6 mol a −1 km −2 ), and according to Kerrick et al. (1995) this value could be used as a baseline to compute convective hydrothermal CO 2 emission from other areas of high heat flow. However, using a similar approach, Frondini et al. (2008) computed a specific CO 2 flux about five times higher (4.8 × 10 6 mol a −1 km −2 ) for Tuscany and Northern Latium (Italy).
In the areas where a flux of deep CO 2 directly affects regional aquifers, most of the deeply generated gas can be dissolved by groundwater because the relatively high solubility of CO 2 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 CO 2 .
The Earth degassing process was extensively studied in Italy and numerous techniques currently used for the estimation of nonvolcanic CO 2 degassing at regional scale were originally developed for the study of the Italian territory.
At Colli Albani, central Italy (Chiodini & Frondini 2001), the CO 2 emission was computed by summing the contributions of CO 2 transported by groundwater and the CO 2 released by two large zones of diffuse soil degassing. The flux of the deeply derived CO 2 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 CO 2 production rate is 2.8 × 10 6 mol a −1 km −2 , of the same order of magnitude as the CO 2 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 centralsouthern Italy. Chiodini et al. (2000Chiodini et al. ( , 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 CO 2 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 CO 2 associated with groundwater circulation, (2) the relationships between regional flux of CO 2 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 CO 2 degassing and advective heat flow.

Geological, hydrogeological and geophysical setting
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   (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 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).

CO 2 degassing in peninsular Italy
Several regional-scale studies showed that the western Tyrrhenian domain is affected by intense CO 2 degassing, which results in numerous cold, CO 2 -rich gas emissions at the surface (e.g. Chiodini et al. 2000Chiodini et al. , 2004bChiodini et al. , 2011Rogie 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(Chiodini et al. , 2010Rogie 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 CO 2 emission ever measured on Earth, which releases amounts of CO 2 (c. 1.7 × 10 10 mol a −1 ) approaching those released from the most active volcanoes (Chiodini et al. 2010).
The idea of measuring CO 2 flux at regional scale derived from the observation that the main Apennine springs emit large amounts of dissolved CO 2 . In principle, the spring discharge can be used to estimate the specific CO 2 fluxes over large areas; that is, the areas of the hydrogeological basins. However, before computing the CO 2 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.

Data
The dataset used in this study (supplementary material) refers to 160 Apennine springs ( Fig. 1) that were sampled in different surveys (Civita 1977;Chiodini et al. 2000Chiodini et al. , 2004bChiodini et al. , 2013Frondini et al. 2012). Water temperature, pH and total alkalinity were determined directly in the field. The chemical (Ca, Mg, Na, K, Cl, SO 4 , NO 3 ) and stable isotope (H 2 O 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 (C ext and δ 13 C ext ; 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 CO 2 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 CO 2 , oxidation of organic material) and sinks (such as mineral precipitation and CO 2 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 13 C/ 12 C ratio of natural water systems is described by the following equations: where I i and O i are the molalities of the ith input and the ith output of carbon in the solution, R is the 13 C/ 12 C isotopic ratio of the solution as a whole, R* is the 13 C/ 12 C isotopic ratio of the ith input species, α i-s 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. , 2011. In the case of 410 F. Frondini et al. carbonate aquifers, the no-sink assumption means that carbonate precipitation and CO 2 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 (C carb ) to TDIC by computing two new variables, C ext and δ 13 C ext : where C ext is given by the sum of the carbon deriving from a deep source (C deep ) and the carbon content of infiltrating waters (C inf ; i.e. the atmospheric CO 2 plus that from biogenic sources active in the soils during the infiltration): from which derives where δ 13 C inf and δ 13 C deep correspond to the isotopic composition of C inf and C deep respectively.
In the above set of equations, only the TDIC and δ 13 C TDIC 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, C carb is given by Equation (9) states that the moles of C carb are equivalent to the moles of Ca and Mg entering the solution from carbonate mineral dissolution, and the negative term -mSO 4 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 C ext is equal to TDIC.
δ 13 C carb 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 C carb and δ 13 C carb allow one to compute C ext and δ 13 C ext (from equations (5) and (6)). The other unknown variables (C inf , C deep , δ 13 C inf and δ 13 C deep ) 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 (C i ), the total output of CO 2 (Q CO 2 ,i ) and the corresponding CO 2 flux (F CO 2 ,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 C deep is based on the analysis of the C ext v. δ 13 C ext 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 C ext values (<3-4 mmol l −1 ) and very negative δ 13 C ext (from −16 to −25‰). These are the springs where C deep is absent (or negligible) and where C ext practically coincides with C inf . Most of the other samples (red circles in Fig. 2) show an enrichment in C ext and higher δ 13 C ext values, indicative of the addition of isotopically heavier carbon, typical of endogenous sources (i.e. the deeply derived C deep ). The C inf 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 C deep for each spring affected by the input of deep CO 2 (red circles in Fig. 2). In this computation, for the few samples where C ext is lower than the mean C inf value, we assumed C deep = 0. Alternatively, C deep 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 C ext of c. 30 mmol l −1 ; dashed line in Fig. 2). This threshold was defined by Chiodini et al. (2011) modelling the effects of CO 2 input on the dissolved carbon species of water at equilibrium with calcite. Above this threshold samples could be affected by CO 2 degassing and calcite precipitation, which can cause significant isotopic fractionation and a decrease of C ext (i.e. it would not be representative of mixing between C inf and C deep ). Figure 2 shows that samples with C ext lower than the no-sink threshold plot in the theoretical mixing field computed by adding to the infiltrating water (C inf = 2.31 ± 0.61 mmol l −1 ; δ 13 C inf = −21.6 ± 2.9‰) variable amounts of C deep with δ 13 C deep ranging Fig. 2. C ext v. δ 13 C ext diagram. Blue circles refer to groundwater with a C ext and δ 13 C ext compatible with dissolution of biological CO 2 only, whereas red circles refer to groundwater with a high C ext and δ 13 C ext compatible with some input of deeply derived CO 2 . The theoretical field resulting from the addition of C deep with δ 13 C deep ranging from −5 to +1‰ to infiltrating water with C inf = 2.31 ± 0.61 mmol l −1 and δ 13 C inf = −21.6 ± 2.9‰ is indicated in grey. The dashed vertical line represents the C ext upper limit for the applicability of the no-sink assumption.

411
Measurement of CO 2 fluxes at regional scale from −5 to +1‰, the typical range of deeply derived CO 2 emitted in Italy by volcanoes, geothermal systems and cold CO 2 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 CO 2 associated with the Apennines groundwater as follows: Q CO 2 ,carb ¼ 2:2 Â 10 10 mol a À1 , Q CO 2 ,inf ¼ 1:6 Â 10 10 mol a À1 and Q CO 2 ,deep ¼ 2:9 Â 10 10 mol a À1 . It is noteworthy that the results show that the deep source of CO 2 is the main source supplying inorganic carbon to the Apennine groundwaters.

Regional map of CO 2 Earth degassing
To derive the map of the regional CO 2 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 CO 2 flux affecting the volcanic aquifer of Mt. Albani, as previously estimated by Chiodini & Frondini (2001). Starting from the C ext values derived from equation (5), the CO 2 flux of each spring (in mol a −1 km 2 ) was computed by multiplying its C ext by the corresponding flow rate and dividing the result by the surface area of its hydrogeological basin (equations (10) and (11)). The CO 2 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 CO 2 diffuse degassing from field measurements (e.g. Cardellini et al. 2003). In this study, each spring was considered as a single CO 2 flux measurement point for the application of the sGs algorithm and the derivation of the regional map. We computed 200 maps of CO 2 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 CO 2 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 CO 2 are shown by colours from light green to red, corresponding to C ext flux values higher than 3 × 10 6 mol a −1 km −2 . This threshold was computed considering a reasonable maximum value for C inf 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 CO 2 /H 2 O 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. (2004bChiodini et al. ( , 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 CO 2rich 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 CO 2 flux values (i.e. C ext ) over the area and subtracting the contribution from the atmospheric and organic sources (i.e. C inf ). This computation assumes a mean recharge water infiltration of 20 l s −1 km −2 , and gives a deep carbon flux of 1.4 × 10 11 mol a −1 for the TRDS and 0.7 × 10 11 mol a −1 for the CDS, corresponding to deeply derived specific CO 2 fluxes of c. 3.7 × 10 6 mol a −1 km −2 and c. 4.7 × 10 6 mol a −1 km −2 , respectively. The total amount of deeply derived CO 2 discharged by this sector of Italy (c. 2.1 × 10 11 mol a −1 ) corresponds to about 2-15% of the estimated presentday global CO 2 fluxes from active volcanoes (Burton et al. 2013),
highlighting the global relevance of the regional CO 2 Earth degassing process in the region. Finally, the computed specific CO 2 fluxes are 4-5 times higher than the value of c. 10 6 mol a −1 km −2 suggested as the reference value for the convective hydrothermal CO 2 emission calculated for high heat flow areas (Kerrick et al. 1995).

Advective heat transport and CO 2 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 CO 2 flux) for 11 aquifers of the central Apennines belonging to the TRDS (Fig. 4). Briefly, the geothermal heat flux (F H ) 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 (T s ) and the infiltration water (T r ), ΔT (in K), is given by where F H is the geothermal heat flux (in W m −2 ), ρ w (kg m −3 ) and C w (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 (m 2 ), q is the spring flow rate (m 3 s −1 ), Δz (m) is the difference in elevation between the water recharge area (Z r ) and the spring (Z s ), 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 (Q H ) and the CO 2 (Q CO 2 ,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 Q H and Q CO 2 ,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 CO 2 -rich fluids. It is noteworthy that the original fluids have a Q H =Q CO 2 ratio similar to that of the hot liquids circulating in the nearby geothermal systems  413 Measurement of CO 2 fluxes at regional scale of Latium (Latera and Torre Alfina; Figs 4 and 5). The total CO 2 output from the two geothermal systems (Table 1) was estimated through hundreds of CO 2 flux measurements performed with the accumulation chambers technique (Chiodini et al. 2007;Lucidi 2010). Considering the temperature and the CO 2 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 Q H of the two systems, computing the total amount of liquid necessary to supply the measured Q CO 2 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 CO 2 flux value of c. 10 6 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 CO 2 fluxes are 4-5 times higher than those derived from TVZ data: the central Italy original fluids have a Q H =Q CO 2 ratio that is much lower than that in the TVZ.

Conclusion
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 CO 2 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 CO 2 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 CO 2 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 CO 2 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 CO 2 fluxes (i.e. c. 3.7 × 10 6 mol a −1 in the TRDS and c. 4.7 × 10 6 mol a −1 in the CDS) are higher than the reference baseline specific CO 2 flux for hot, hydrothermal zones of the Earth (10 6 mol a −1 ; Kerrick et al. 1995). This could reflect a real strong CO 2 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, CO 2 -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 CO 2 (c. 2.1 × 10 11 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 CO 2 budget of Italy (Fig. 6) including the SO 2 -bearing plumes of active volcanoes (Etna, Stromboli, Vulcano; Burton et al. 2013) and the volcanoes with hydrothermal activity where SO 2 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 CO 2 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 CO 2 in Italy.
The results obtained in Italy demonstrate that total CO 2 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 CO 2 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 CO 2 degassing, quantitative data on CO 2 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.
Funding This work was funded by the Istituto Nazionale di Geofisica e Vulcanologia. Fig. 6. Output of deeply derived CO 2 in Italy from volcanoes with hydrothermal systems (yellow), volcanoes emitting SO 2 (red), and from regional diffuse degassing structures in central Italy (blue).