The hydrothermal system of Bagni San Filippo (Italy): fluids circulation and CO₂ degassing

Central-southern Italy is affected by an active and intense process of CO₂ Earth degassing from both active volcanoes and volcanically not active areas, and is one of the few regions of the world where detailed mapping and quantification of the CO₂ fluxes have been performed in the last two decades (Chiodiniet alii, 1999; Chiodiniet alii, 2000; Chiodiniet alii, 2004; Frondiniet alii, 2008; Frondiniet alii, 2019). In particular, Chiodiniet alii (2004), based on the carbon balance of regional aquifers, produced the first regional map of CO₂ degassing in central and southern Italy, showing the occurrence of two large scale CO₂ anomalies on the Tyrrhenian (western) side of the Italian peninsula. The two structures, named Tuscan Roman Degassing Structure (TRDS) and Campanian Degassing Structure (CDS), discharge to the surface 1.69×10⁴ t d^-1 (i.e., 6.16 Mt y^-1) and 8.44×10³ t d^-1 (i.e., 3.08 Mt y^-1) of CO₂, respectively. The total amount of CO₂ released to the atmosphere through diffuse regional degassing (9.24 Mt y^-1) is the main component of the geological CO₂ budget in Italy, being higher than the amount of CO₂ discharged by active Italian volcanoes (7.32 Mt y^-1, Frondiniet alii, 2019 and references therein). According to Frezzottiet alii, (2009), the degassing process is related to decarbonation or melting, at pressure greater than 4 GPa (130 km), of carbonate-rich lithologies subducted in the mantle.

In addition to CO₂-rich groundwater, both TRDS and CDS are also characterized by the presence of many CO₂ vents and areas of anomalous soil CO₂ diffuse degassing (e.g., Chiodiniet alii, 1999, 2007, 2011; Rogieet alii, 2000, www.magadb.net) that are generally fed by buried carbonate reservoirs, covered by low permeability formations, where the gas produced at depth accumulates before the expulsion to the surface (e.g., Chiodiniet alii, 2007; Collettiniet alii, 2008).

These gas emissions can represent a severe risk for humans and animals: at ambient temperature carbon dioxide is about 1.5 times heavier than air. During stable atmospheric conditions the gas can accumulate in topographic depressions reaching high, often lethal, concentrations (e.g., Chiodiniet alii, 2008; Costaet alii, 2008; Chiodiniet alii, 2010).

The Bagni San Filippo hydrothermal system (BSF) is a paradigmatic example of the CO₂ degassing process affecting western central Italy. BSF is located within TRDS, in the Mt. Amiata volcanic and geothermal area (AVGA, Fig. 1a). The degassing process at BSF results at the surface with focused gas vents, zones of diffuse soil CO₂ degassing, CO₂-rich thermal waters with temperatures from 23 to 48°C and large travertine deposits (Frondiniet alii, 2009). The areas of stronger CO₂ emission can be very dangerous both for humans and animals, as testified by a lethal accident occurred on November, 2003 in one of the topographic depressions of the area (Fig. 1b; https://ricerca.gelocal.it/iltirreno/archivio/iltirreno/2003/11/20/LF3PO LF302.html; Chiodiniet alii, 2008).

In this work, we combine i) the results of a detailed survey on the CO₂ soil diffuse degassing, ii) a hydrogeochemical investigation of the BSF thermal waters, and iii) a study of the chemical and isotopic features of the gas emissions in order to derive a comprehensive conceptual and quantitative model of the system. Besides this main aim, the definition of the most intensely degassing areas may contribute to mitigate the risk associated to the CO₂ release.

The study is based on original data collected from December 2003 to May 2004, after the fatal accident of November 2003.

BSF is located in southern Tuscany (Italy), in the NE sector of the AVGA, on the western side of the Siena-Radicofani graben (Fig.1a). The geological setting of the region is the result of the post-thickening evolution of the inner sector of the Northern Apennine, characterized by extensional processes accompanied by coeval magmatism active since early-middle-Miocene (e.g., Carmignaniet alii, 1995; Batiniet alii, 2003; Brogiet alii, 2005; Marroniet alii, 2015a). During Cretaceous to early Miocene the study area was affected by compressional tectonics producing the stacking of the tectonic units belonging to the Northern Apennine paleogeographic domain (Brogi & Liotta, 2006). Subsequently, from Miocene to present, the compressional structures were deformed by extensional tectonics producing the NNW-SSE basins (e.g., Siena-Radicofani graben), filled by marine to continental sediments, and the associated magmatism (Ferrariet alii, 1996; Brogiet alii, 2005, 2015; Brogi & Liotta, 2006).

AVGA is characterized by the following tectonic stratigraphic units, from bottom to top: i) the Paleozoic crystalline metamorphic basement made up of graphitic phyllites, meta-sandstones, siltites, meta-carbonates and graphite-rich metasediments (Elter & Pandeli, 1991). The formations of the basement are not exposed at the surface but have been drilled in the Piancastagnaio area (Chiodiniet alii, 1988; Elter & Pandeli, 1991); ii) the Tuscan nappe, that includes sedimentary rocks ranging from Late Triassic evaporites to Jurassic carbonate platform, Cretaceous-Oligocene pelagic sediments and late Oligocene-early Miocene turbidites (Brogiet alii, 2015; Marroniet alii, 2015a); iii) the Sub-Ligurian and Ligurian Tectonic Units, composed by ophiolites, shales, sandstones, calcarenites and marls (Pandeliet alii, 2005; Marroniet alii, 2015b); iv) the Neogene sedimentary deposits, mainly represented by sands and clays; v) the Mt. Amiata volcano, consisting on trachytic to olivine latitic lava flows and domes emplaced in a very short time between 305 and 231 ka (Ferrariet alii, 1996; Conticelliet alii, 2015).

In the AVGA are present two water-dominated geothermal reservoirs sited in the southern part of the Amiata volcano. The shallower geothermal system is hosted by the carbonate-evaporite formations at depths of 500-1000 m and shows fluids with temperature ranging from 200 to 240 °C, whereas the deeper one is localized in the metamorphic basal complex at depths of 2000-4500 m and temperatures between 300 and 350 °C (Calamaiet alii, 1970; Bertiniet alii, 1995). The geothermal reservoirs are tapped by the wells of two active geothermal fields (Piancastagnaio and Bagnore, Fig.1) with a total installed capacity of 121 Mw_e (Manzellaet alii, 2019). The geothermal gradient of the area is generally between 80 and 150 °C km^-1 with local peaks up to 300 °C km^-1 (Baldiet alii, 1994) and the heat flow is 200 mW m^-2, on average, reaching values of about 600 mW m^-2 in correspondence of the geothermal fields (Chiodiniet alii, 1988; Brogiet alii, 2005).

The hydrogeological setting of AVGA is dominated by a regional aquifer hosted by the Mesozoic carbonate-evaporite formations of the Tuscan nappe and by a shallower aquifer hosted by fractured volcanic rocks (Boniet alii, 1986; Chiodiniet alii, 1988). The two aquifers are separated by a low permeability complex represented by Ligurian units and Neogene sediments, acting as an aquiclude. In AVGA the carbonate-evaporite aquifer is almost completely confined under the low permeability complex and is characterised by discontinuous recharge areas corresponding to the outcrops of the carbonate-evaporite formations. On the southern and western side of Mt. Amiata the carbonate-evaporite structures represent the main geothermal reservoir and their lateral continuity permits the circulation of deep fluids over a large area (Celatiet alii, 1990; Minissaleet alii, 2002; Frondiniet alii, 2009).

One of the few zones of AVGA where the Tuscan nappe sedimentary formations outcrop is the BSF area. In particular the carbonate and evaporite formations of the Tuscan Nappe outcrop on the western side of the study area (Mt. Poggio Zoccolino, Fig.1) while its western side is characterized by the Ligurian tectonic Units and its central part is characterized by Pliocene marine and continental deposits topped by a large Quaternary travertine plate, deposited by the thermal waters of the area. The aquifer hosted by the Tuscan Nappe, that is locally recharged by the Mt. Poggio Zoccolino (1035 m a.s.l.) structural high, discharge in the BSF area more than 35 L s^-1 of thermal water with temperatures up to 48 °C.

A multi-parametric geochemical study of the BSF area was performed from December 2003 to May 2004. During this period the area was surveyed for the diffuse soil degassing of CO₂, chemical and isotopic composition of gas discharged, chemical and isotopic composition of groundwater and chemical composition of the gas dissolved by groundwater.

From December 2003 to May 2004 four CO₂ flux surveys were performed at BSF area using the accumulation chamber method (Chiodiniet alii, 1998), using a LI-800 infrared spectrometer as CO₂ detector operating in the range 0-20000 ppm. During December 2003, an area of about 1.7 km² (Fig. 1b) was investigated with 595 measurements of CO₂ flux, homogeneously distributed and on average spaced 50 m (General survey, GS; Table 1). Based on the results of the GS data analysis, 3 detailed CO₂ flux surveys were performed at Palazzo (PA), Lavinate (LA) and Bollore (BO) areas (Fig. 1b), by using a measurement spacing from 5 to 20 m. In detail, 423 CO₂ flux measurements were performed at PA in March-April 2004, 325 at Lavinate in May 2004 and 303 at Bollore in March 2004 (Table 1). The complete datasets of the measured CO₂ fluxes are reported in the Supplementary Material.

The CO₂ flux data were elaborated using both statistical and geostatistical tools. The graphical statistical method (GSA method, Chiodiniet alii, 1998), based on the Sinclair’s partitioning method (Sinclair, 1974), was used to characterise the statistical distribution of the CO₂ flux and the sources of the diffuse degassing. In fact, in volcanic and geothermal areas CO₂ soil diffuse degassing is frequently fed by multiple gas sources, such as the biological and the deeply produced CO₂. This multiple origin can result in a polymodal distribution of the CO₂ flux, which plots on a log-probability plot as a curve with n-1 inflection points when n lognormal populations overlap. Conversely, a single log-normal population would result in a “straight” line (Sinclair, 1974). The GSA allows to partition the polymodal distribution into individual log-normal populations and to estimate the relative proportion (f_i) between different populations, their mean (M_i) and standard deviation (σ_i). Since the calculated M_i refers to the logarithm of CO₂ flux values, the mean of the CO₂ flux was then estimated using a Monte Carlo simulation procedure. The uncertainty of this estimate is given in terms of 5^th and 95^th percentiles of the mean CO₂ flux values returned by the Monte Carlo simulation.

A geostatistical approach based on the sequential Gaussian simulations (sGs method; Cardelliniet alii, 2003) was used to map the diffuse CO₂ degassing and to quantify the total CO₂ output when possible. The sGs method consists of the production of numerous equiprobable realizations of the spatial distribution of the CO₂ flux, here performed using the sgsim algorithm of the GSLIB software library (Deutsch & Journel, 1998). The CO₂ flux is simulated on a regular grid, according to the variogram model defined fitting the experimental variogram of the normal scores of the CO₂ flux, which provides a description of how the data are spatially correlated. The variogram model is given in terms of nugget, range and sill parameters, where the nugget represents the small-scale variability, the range represents the distance within which data are correlated and the sill is the plateau the variogram reaches for a distance equal to the range. The produced realizations were then used to draw the maps of the CO₂ diffuse degassing, which are here provided as probability maps. The probability map consists of a map of the probability that, among all the realizations, the simulated CO₂ flux at any location (i.e., at grid nodes) is above a cut-off value. Selecting the threshold value of the biogenic CO₂ flux as a cut-off, the probability map is used to highlight the area interested by the emission of deeply derived CO₂ (i.e., the diffuse degassing structure, DDS; Chiodiniet alii, 2001, Cardelliniet alii, 2003). According to Cardelliniet alii (2003), the DDS is here considered as the area where the probability that the simulated CO₂ flux is higher than the biogenic CO₂ flux threshold is over 50%.

The total CO₂ output is defined for each area as the mean of the total CO₂ output computed for each realization, which is computed by summing the products of the simulated CO₂ flux value at each grid cell by the cell surface. The uncertainty of the estimate is assessed from the 5^th and 95^th percentiles of the total CO₂ output values computed for all the realizations.

Ten springs, representing the main thermal groundwater outflows of BSF, were sampled in February 2004. Location of sampling points is shown in Fig. 1b. For each water sample, temperature, pH, Eh, electrical conductivity and HCO₃ were measured in the field. HCO₃ concentration was determined by acid titration with 0.01N HCl using methyl orange as indicator. Water samples for chemical and isotopic (δD and δ¹⁸O) analyses were collected in 100 mL and 250 mL HDPE (high-density polyethylene) bottles. One of the 100 mL aliquots was filtered upon sampling through 0.45 μm membrane filters and then acidified with 1% of 1:1 diluted HCl. Water samples for isotopic composition of total dissolved inorganic carbon (TDIC) were collected in 1000 mL glass bottles. The dissolved carbon species were precipitated directly in the field as SrCO₃ by adding SrCl₂ and NaOH in solid state to the water sample. Carbonate precipitate was recovered by filtering the sample and washed with distilled water in a CO₂-free atmosphere in the laboratory. Carbon isotopic analyses of the TDIC (δ¹³C_TDIC) were performed on the CO₂ gas liberated from the SrCO₃ precipitate by reaction with 100% H₃PO₄.

Major ions were determined at the laboratory of Perugia University. Calcium and Mg concentrations were determined by atomic absorption (AA) flame spectroscopy on the acidified sample while Na and K were determined by atomic emission (AE) flame spectroscopy, using an Instrumentation Laboratory aa/ae spectrophotometer 951. Cl and SO₄ were determined by ion chromatography using a Dionex DX-120.

Carbon, oxygen and hydrogen isotopes analyses were performed at the Geochemistry Laboratory of INGV-Naples using a Finnigan Delta plusXP continuous flow mass spectrometer coupled with GasbenchII device. Isotopic compositions are given in δ‰ notation per mill versus V-PDB for carbon and versus V-SMOW for hydrogen and oxygen (analytical errors are: δD ± 1‰, δ¹⁸O ± 0.08‰ and δ¹³C ± 0.06‰).

The main accessible gas emissions were sampled in November 2003 and February 2004 (Fig. 1b). Two samples were collected at each sampling site using a 50 mL and a 250 mL glass vials, both equipped with a Teflon stopcock fixed with a rubber O-ring. The 250 ml glass vials were pre-evacuated and filled with about 50 ml of a 4N NaOH solution (Giggenbach, 1975; Giggenbach & Goguel, 1989). Chemical analysis of gas samples were performed at the Geochemistry Laboratory of INGV-Naples.

Carbon dioxide and sulphur species, absorbed in the NaOH solution in the vial, were analysed after oxidation with H₂O₂, by acid-base titration and by ion-chromatography, respectively. The gas species not absorbed by the NaOH solution (i.e., He, Ar, O₂, N₂, H₂, and CH₄) present in the vial headspace, were analysed using a gaschromatograph (Agilent Technologies 6890N), make up of two channels equipped with two PLOT columns (MS 5Å, 30m x 0.53mm x 50μm; He and Ar as carrier gas) and Thermal Conductivity Detectors (TCD). Carbon monoxide was determined on the 50 mL dry gas sample by means of gas-chromatographic separation (Agilent Technologies 6890N) with a PLOT column MS 5 Å 30m x 0.53mm x 50μm (He as carrier gas) using a high-sensitivity reduced gas detector (Trace Analytical RGD2; detection limit 0.05 ppm).

The isotopic composition of CO₂ was determined on 50 mL dry gas sample by the same instrumentation used for δ¹³C_TDIC analyses

The measured CO₂ flux in the GS survey (Fig. 2) varies of orders of magnitude from a minimum of ~ 0.5 g m^-2 d^-1 to ~ 70,000 g m^-2 d^-1 (average value of 533 g m^-2 d^-1). The high CO₂ fluxes (i.e., higher than 200-300 g m^-2 d^-1) are strongly variable, represent a relatively low proportion of the entire dataset (~ 5%) (Fig. 3a) and are generally located nearby the main gas vents (Fig. 2).

The probability plot (Fig. 3a) shows that the statistical distribution of the CO₂ flux is the result of the overlapping of 2 log-normal CO₂ flux populations. For the low flux population (low flux population general survey, LFP-GS, Table 2) we compute a mean value of 9.4 g m^-2 d^-1 (8.8 - 10.0 g m^-2 d^-1) while for the high flux population (high flux population general survey HFP-GS) the Monte Carlo approach returns a mean value of ~14000 g m^-2 d^-1 with a very large variability of the estimation (2000 -41000 g m^-2 d^-1). Such high fluxes are clearly fed by a CO₂ deep source. On contrary the mean flux of LFP-GS well falls within the range of the CO₂ fluxes produced by the biological activity in the soil in a wide variety of ecosystems (mean CO₂ flux from 0.2 g m^-2 d^-1 to 21 g m^-2 d^-1, Cardelliniet alii, 2017; Viveiroset alii, 2010 and reference therein). For this reason the LFP-GS population can be considered as representative of the biological CO₂ flux for the period of the GS.

The variogram of the normal scores of the CO₂ flux of the GS shows a high nugget value (0.75, Fig. 4a) suggesting that the CO₂ flux has a large spatial variability at a scale smaller than the spacing between the measurements. This is possibly due to the low permeability of the soils that forces the gas to flow mainly through fractures.

A total deep CO₂ output of 849 t d^-1 was computed for the GS survey (Table 3) by subtracting to the total CO₂ output estimated by sGs approach (865 t d^-1, see Methods) the biological CO₂ contribution of 9.4 g m^-2 d^-1 (Table 2). However, the reliability of this estimate is low because the high small-scale spatial variability of the CO₂ fluxes, in comparison to the spacing between the measurements. As already noted in other areas, a large measurement spacing of highly spatially variable soil CO₂ fluxes can in fact cause an overestimation of the total CO₂ output (e.g., Viveiroset alii, 2010) and/or a large uncertainty of the estimate (Cardelliniet alii, 2003).

In order to have more reliable estimations, we decided to perform detailed surveys of the most anomalous zones. We selected the anomalous areas of Bollore (BO), Palazzo (PA) and Lavinate (LA) based on the CO₂ flux probability map (Fig. 2) which has been drawn assuming a CO₂ flux cut-off value of 25 g m^-2 d^-1 (i.e., the 95^th percentile of the LFP-GS population, Fig. 3a).

The measured CO₂ fluxes of the detailed surveys distribute from few g m^-2 d^-1 to more than 100000 g m^-2 d^-1, similar to that of the GS (Table 1). For all the three surveys the probability distribution of CO₂ fluxes is interpreted as the results of the overlapping of three CO₂ flux populations (Fig. 3b-d). The populations with the lowest mean CO₂ fluxes are characterised by very similar and reliably estimated mean values: 23.1 g m^-2 d^-1, 22.9 g m^-2 d^-1 and 19.4 g m^-2 d^-1 for LFP-PA, LFP-LA and LFO-BO, respectively (Table 2). These values are our best estimate of the biogenic CO₂ contribution that in any area resulted higher than the biological background CO₂ flux for the GS survey (9.4 g m^-2 d^-1), likely because the different season of the surveys (the GS was performed during the winter while PA, LA and BO surveys were performed during the spring). The biological soil CO₂ flux in fact varies with the environmental conditions (e.g., air/soil temperature, precipitation and soil water content) and is generally higher in the spring-summer when can be up to of 4-5 times higher respect to the winter (e.g., Davinsonet alii, 1998; Raich & Schlesinger, 1992). The deeply derived CO₂ flux is clearly identified by the other two populations, HFP and IFP (intermediate flux population), which in each survey are both characterized by a mean CO₂ flux too high for the biogenic source (Table 2). The presence of two deeply derived CO₂ flux populations could reflect the relative prevalence of diffusive or advective CO₂ transport mechanisms.

The computation of the CO₂ output from the detailed surveys was computed, when possible, using the sGs approach that is based on the spatial correlations among the measurements. The variograms of the LA and PA surveys dataset (Fig. 4b, c) show nugget effect (0.6 and 0.54) lower than that of the GS variogram model (0.75) indicating that the spacing between the measurements of the detailed surveys (smaller than that used in the GS) is more adequate to catch the spatial structure of the CO₂ flux. Conversely, for the BO survey, even if a relatively high measurements density was used, the experimental variogram shows a high nugget component (higher than 0.8) suggesting a strong small-scale variability of the CO₂ flux. This variability is particularly strong in the NE sector of the surveyed area (sector A, Fig. 2) as suggested by the “pure nugget effect” variogram of CO₂ fluxes measured in this area Fig. 4d) which indicates the lack of any spatial correction. Sector A is in fact dominated by the presence of a hill made of travertine whit null or scarce soil, and the CO₂ mainly flow through fractures resulting in a variability of the flux of orders of magnitude at a very small scale. Contrariwise, the variogram of sector B (i.e., the remaining BO area, sector B in Fig. 2) shows a relatively low nugget (Fig. 4d).

The CO₂ flux was mapped by sGs for LA, PA and only for the sector B of BO performing 200 simulations for each area (using 2 m × 2 m cells) and are reported in Fig. 5 as probability maps. At Palazzo the CO₂ flux map (Fig. 5a) highlights the presence of a series of DDSs roughly distributed along a NNE-SSW band, a direction coincident with the orientation of the fault bordering Mt. Poggio Zoccolino (Fig 1b). A fault damage zone interested by a strong CO₂ degassing is evident in the southern part of the area and along the NNE-SSW direction are located numerous CO₂ vents and a large morphological depression (Fig. 5a) where a strong CO₂ vent at its bottom causes high CO₂ air concentrations (up to 50 % vol CO₂ was measured during the survey). At Lavinate the CO₂ flux map shows a well-defined DDS that mainly develops along the creeks directions (e.g., Torrente Lavinate and Fosso Rondinaia, Fig. 5b) which follow the directions of the main faults (Fig. 1b). In particular, the NE-SW direction corresponds to Bagni San Filippo Fault which represents the prosecution of the Mt. Amiata fault along which developed the main volcanic centres (Brogi & Fabbrini, 2009, Brogiet alii, 2010). For the sector B of Bollore the CO₂ flux map (Fig. 5c) highlights a roughly E-W elongated DDS in continuity with the orientation of the travertine hill, where the highest fluxes were measured. All these features of the CO₂ degassing suggest the main role of the tectonic structures on the upward transfer of the deeply produced CO₂.

Based on the sGs results, the total CO₂ output was computed for PA, LA and BO sector B resulting in 21.4 t d^-1, 65 t d^-1 and 6.6 t d^-1 respectively (Table 3). The deep CO₂ output was then estimated subtracting the biological CO₂ contribution (i.e., the mean CO₂ flux of the respective LFP populations) to the total CO₂ output and results of 19.8 t d^-1, 63.8 t d^-1 and 5.4 t d^-1 for PA, LA and BO sector B respectively (Table 3). For sector A of BO the total deep CO₂ output was estimated by the GSA method (Chiodiniet alii, 1998) because of the lack of spatial correlation among the measured fluxes. Summing the estimates for the two sectors (Table 3), a total deep CO₂ output of 47 t d^-1 is obtained for the entire Bollore area.

The total deep CO₂ output from the BSF area is computed by summing the above BO, PA and LA results (130 t d^-1) to the flux estimated for other two degassing areas not measured in detail, i.e. Campo la Villa (CLV) and Lavinate-East (LE) (Fig. 2). For these two areas, rough estimations were obtained by the sGs applied to the GS data multiplied by an empirical correction factor (cf) estimated as the ratio between the results of the detailed campaigns at BO, PA and LA and the results over the same areas using the GS data. The GS estimates are systematically higher than those of the detailed surveys and the cf varies from 0.3 to 0.85. Using these values, the cumulative contribution of deeply derived CO₂ from CLV and LE results at 30 t d^-1 to 80 t d^-1. Therefore, our best estimation of the total deep CO₂ output diffusively emitted at BSF is of 160-210 t d^-1. This value is the highest among the measured gas manifestations of AVGA including Selvena (17 t d^-1; Rogieet alii, 2000), Banditella (10 t d^-1; Frondiniet alii, 2009) and Ermeta (12 t d^-1; Nisiet alii, 2014).

As a general consideration, we remark that the comparison of the detailed surveys with the GS data indicates a large overestimation of the GS-based total emission possibly arising from the large spacing between the GS measurements respect to the dimension of the DDSs.

The sampled thermal springs (BSF springs, Table 4) are characterized by flow rates ranging from 0.01 to 24 L s^-1 and discharge temperatures from 23 to 48 °C. Their position in the Langelier Ludwig diagram (Fig. 6) indicates that the BSF thermal waters can be interpreted as mixtures of a Ca-SO₄ component with a Ca–HCO₃ hydro-type typical of the sedimentary aquifers of the region (both limestones and sandstones, Chiodiniet alii, 1988; Frondiniet alii, 2009). Fig. 6 suggests also that there are no significant mixing processes neither with groundwaters of the volcanic aquifer nor with a liquid phase from the high enthalpy geothermal systems of AVGA.

The temperature of the springs is positively correlated with SO₄ and total dissolved carbon (TDIC) while it is negatively correlated with Cl (Fig. 7), suggesting that the groundwaters are heated up by hot fluids rich in carbon and sulphur and poor in Cl. Such characteristics are typical of hydrothermal-geothermal vapours: most of the carbon content of the waters would derive from the deep CO₂ of the geothermal vapours and the sulphate content possibly from both dissolution of anhydrite and oxidation of deeply derived H₂S.

The enthalpy-chloride diagram (Fig. 8a) supports this interpretation because the BSF thermal springs fit the theoretical mixing of groundwaters with condensates of hydrothermal vapours, or directly with steam (red arrow in Fig. 8a), rather than with geothermal liquids (for comparison in the figure is considered the mean values of several Amiata geothermal wells, Minissaleet alii, 1997). In the following we investigate in detail the mixing condensates-groundwaters because the alternative interpretation of the direct input of steam into the aquifer feeding the BSF thermal spring, even if possible, is less supported by the enthalpy-chloride data of the thermal springs (Fig. 8a).

The mixing between groundwaters and condensates (red line in Fig. 8a) has been drawn as the best linear fit of the BSF thermal waters and the vapour condensate at 100 °C (enthalpy = 419 J g^-1) assumed with no chloride. At low enthalpy values, the regression line points to a groundwater sample (green point in Fig. 8a) that, among the many samples from sedimentary aquifers reported in Frondiniet alii, (2009), is that closest to Bagni San Filippo area. The mixing model assumes that the enthalpy (temperature) of the BSF thermal waters, and in particular of the hottest spring 33 whose flow rate is 70% of the total BSF discharge (Fig. 8a), is representative of the system at depth. The estimations of the deep temperatures were attempted using the log (Ca/Mg) geothermometer based on the equilibrium of the solutions with calcite and dolomite (Mariniet alii, 1986). The geothermometer provides estimations of 47 °C - 52 °C practically coincident with the measured values of the springs of highest temperatures. This coincidence regards, in particular, the spring 33 suggesting that only a minor loss of heat occurs during the ascent of the water from depth, corroborating the reliability of the mixing model.

Interestingly in the classical δD vs δ¹⁸O diagram (Fig. 8b), the BSF springs plot close to the possible meteoric lines but with an evident shift towards lighter oxygen isotopic compositions. This behaviour can be explained by the input into groundwaters of meteoric origin of a vapour rich in CO₂ (i.e., of condensed steam and CO₂). In one hand the steam produced by boiling is depleted in ¹⁸O and ²H with respect to the parent liquid, in the other hand the CO₂ can exchange the oxygen with the water causing a fractionation again toward lighter isotopic compositions.

Aqueous speciation calculations, computed by means of the PHREEQC code (Parkhurst & Appelo, 1999), show that the BSF waters have partial pressures of CO₂ (P_CO2) ranging from 0.1 to 1.5 bar (Table 5) and that are all slightly undersaturated with respect to gypsum and saturated or supersaturated with respect to calcite and dolomite. This supersaturation with respect to carbonate minerals suggests that the waters at depth have higher contents of CO₂. For the spring 33, assumed as representative of the system (see above), we estimated an original P_CO2 of ~7 bar by adding CO₂ to the solution (PHREEQC model) until reaching equilibrium conditions with calcite.

Summing up, we infer the presence at shallow depth of a system recharged by a mixture of groundwaters and condensates, a system with a temperature close to 50 °C and high P_CO2. Most of the incondensable gases of the original vapour phase would be released as a free gas phase in the zones where the vapour enters the aquifer while a main part of the dissolved fraction would be released during the ascent and the depressurization of the waters. The sounding of this model is suggested also by the presence in the area of numerous gas emissions (Fig. 2, Table 6) and areas of CO₂ soil diffuse degassing (Figs. 2 and 5) that would be fed from both a free gas phase present at depth and from the ascending degassing solutions.

The gas emissions have a CO₂ dominated composition (Table 6) and relatively high contents of CH₄ and nonatmospheric N₂. The CH₄ concentrations range from 1.4% to 1.9% while the N₂ concentrations are of 2% - 2.5% with N₂/Ar ratios from 2200 to 6300 (excluding the samples BSF4 and BSF10 whose N₂/Ar of ~ 100 indicates an evident air contamination during sampling). The relative contents of N₂, He and Ar of BSF gases fall in the range of other gas emissions of the TRDS (Fig. 9) that overall are enriched in N₂ respect to He when compared with the gases from of the active quiescent volcanoes of the Italian Peninsula (e.g., Vesuvius, Campi Flegrei and Ischia). The relatively high N₂/He (and N₂/Ar) ratios of the TRDS gases are ascribed to N₂ enrichment due to the thermal decomposition of organic matter recycled through subduction (e.g., Chiodiniet alii, 2007).

Also the carbon isotopic composition of the CO₂ (δ¹³C_CO2 from -2.8‰ to -2.3‰, Table 6) is in the typical range of TRDS gases (δ¹³C_CO2 from -4‰ to -2‰, Chiodiniet alii, 2007) which have been interpreted as the result of a mixing between a CO₂ produced at depth from decarbonation of carbonates and a more negative mantle CO₂ (Chiodiniet alii, 2000; Chiodiniet alii, 2004). We note that the carbon isotopic compositions of the BSF gas emissions are significantly lighter than the carbon dissolved in the thermal waters (δ¹³C_TDIC from -0.8‰ to +2.2‰, Tab 4). In order to investigate the processes controlling such differences we applied a gas-water-rock interaction model (GWRI) based on i) the PHREEQC code to manage the chemical composition of the waters, ii) the Henry law to compute the composition of the dissolved gases, and iii) data from Wigleyet alii (1978) to compute the carbon isotopic fractionations during the simulated processes. The model simulates the composition of a normal groundwater affected both by the input of a deep gas and by degassing at a temperature of 50°C, that is the temperature measured in the main springs (sample 33) and returned by the log (Ca/Mg) geothermometer.

According to the results of many investigation in the carbonate aquifers of the Apennines (e.g., Frondiniet alii, 2019 and references therein) the total dissolved inorganic carbon (TDIC) of the normal groundwater is set at 0.0035 mol L^-1 with a carbon isotopic composition of -10‰, a value controlled by the mixing among the atmospheric carbon, the biogenic carbon of the soil CO₂, and the carbon from calcite dissolution (δ¹³C = +2.5‰). The chemical and isotopic composition of the deep gas has been assumed as that of BSF9, the gas emission less affected by air contamination (i.e., N₂/Ar = 6260, 2 orders of magnitude higher than air and the highest among the sampled gas emissions). In the model the deep gas is added until the solution reaches a P_CO2 of 7 bars (see above) and allowing the solution to dissolve calcite at equilibrium conditions (curve A in Fig. 10). The solution is then depressurised (without allowing calcite precipitation as suggested by the supersaturation of the springs) down to the lowest measured P_CO2 values (i.e. 0.2 bar, curve B in Fig. 10). In general the measured δ¹³C_TDIC and TDIC of BSF thermal springs follow this trend of depressurizing solutions that have reached relatively high P_CO2 at depth. In particular, the model results fit the data of the main and hottest spring 33 (Fig. 10).

Our conceptual model of the BSF system is sketched in Fig. 11. The thermal springs of BSF are recharged in the Mt. Poggio Zoccolino area where permeable carbonate rocks outcrop. The recharge area is separated from the discharge area, that is less than 2 kilometres distant (Fig. 1), by a zone (gas emissions zone in Fig. 11) where the hydrothermal condensates likely enter into the shallow aquifer.

Note that in this model the vapour condensation should occur at depths greater than that of groundwater circulation because our evidence is the input of a hot liquid plus a separated gas phase rather than of a high enthalpy vapour phase (Fig. 8a). This zone, separating the recharge and discharge areas, is the zone where most of the incondensable gases are emitted as testified by the presence of numerous free CO₂ emissions. This is also the zone that was surveyed for diffuse soil CO₂ fluxes (Figs.1 and 2), whose results indicate a total deep CO₂ flux of 160-210 t d^-1 (see section 4.1). Considering the CO₂ emission from the gas vents (estimated in 50-100 t d^-1, Tassiet alii, 2009) and the amount of CO₂ dissolved and transported by the waters (6.4 t d^-1 computed as the sum of the TDIC of the springs multiplied by their flow rates) our best estimation of the total deep CO₂ emission at BSF is at 226-326 t d^-1.

The condensate flow rate is evaluated from the fraction of condensates that feed the BSF thermal springs. The enthalpy and chloride mass balance applied to the mean springs value (weighted by the respective flow rate), assuming as pure end-members the condensate and the normal groundwater (Fig. 8a), indicates that ~11.2 kg s^-1 (965 t d^-1), of the total 35 kg s^-1 thermal water discharge, are of condensates. Summing up, the deep CO₂ emission is of 226-326 t d^-1 and the condensate flow rate is of 965 t d^-1, values that point to the input of a vapour phase with a CO₂ molar fraction X_CO2 of 0.08-0.12 (X_CO2 = mol_CO2/(mol_cond + mol_CO2)). This is plausible X_CO2 range for a geothermal vapour (e.g., X_CO2 measured in several geothermal wells of Amiata region range from 0.07 to 0.12; Chiodini & Marini, 1998). The total thermal energy release involved in the condensation (latent heat of condensation) and in the cooling of the condensates would result at ~ 29 MW. Of this total thermal energy release, ~ 4 MW are associated with the BSF thermal discharge while the remaining 25 MW would refer to deeper processes, i.e. the steam condensation occurring at depths greater than groundwater circulation.

These results are dependent on the selected mixing model, if, e.g., we assume the alternative model of the input of steam (condensation directly into the aquifer, dashed red arrow in Fig. 8a) the estimated flux of condensates formed into the aquifer, X_CO2 of the pre-condensed vapour and the thermal energy release would be of 146 t d^-1, 0.38-0.48, and 4.5 MW, respectively.

Thermal waters of Ca HCO₃-SO₄ composition and temperatures up to 48°C, and numerous cold emissions of large amounts of CO₂, both from soil diffuse degassing and vents, occur in Bagni San Filippo area, an area of few km². We individuate in ascending geothermal vapours the source of both the gas emission and the high temperatures of the waters. The thermal waters are in fact mixtures of condensates of hydrothermal steam (or, less likely, of steam) with locally recharged waters, whereas the cold gas emissions are fed by the incondensable gases originally contained in the ascending vapours. This conceptual model is supported by the temperature of the springs that is negatively correlated with the chloride and positively with sulphate and dissolved carbon suggesting that groundwaters are heated up by hot fluids rich in carbon and sulphur and poor in Cl, such as typical hydrothermal vapours. In agreement with this model, the isotopic compositions of the carbon dissolved in the thermal waters and that of the cold CO₂ emissions point to a unique deep carbon source (δ¹³C = -2.3‰). In a first stage, the condensates and the deep-originated CO₂ enter the shallow aquifer causing a temperature increase to 50°C and elevated P_CO2 (~7 bar). In a second phase, the ascending thermal solutions depressurize and release a great part of the gas contributing to the observed CO₂ emission. We computed at 11.2 kg s^-1 the amount of condensate transported by the springs, whose total flow rate is of 35 kg s^-1, and at ~ 29 MW the total thermal energy release. Most of this heat flux (~ 25 MW) is associated with the condensation occurring at depths greater than groundwater circulation. Assuming alternatively that the condensation process occurs directly into the aquifer (i.e., input of steam instead of condensate, less probable but possible) the flux of condensates would reduce at 1.7 kg s^-1, with a thermal energy release at 4.5 MW.

The location of the deep CO₂ gas emission, both vents and diffuse degassing structures, is strongly controlled by the tectonic lineaments. Diffuse CO₂ flux surveys were performed to quantify the deep CO₂ emission from this area. The analysis of the data of a first general measurement campaign indicates a spatial variably of the fluxes larger than the spacing between the measurements, a fact that implies a large overestimation of the deep CO₂ emission (~ 850 t d^-1). For this reason we performed specific, more detailed, campaigns focussed on three of the most anomalous zones obtaining a total deep CO₂ flux of 130 t d^-1 from diffuse degassing. Our best estimation of the total deep CO₂ emission becomes of 226-326 t d^-1 when the fluxes estimated for other two anomalous degassing areas, from gas vents and that transported as dissolved carbon by the thermal waters are considered. Combining the estimated fluxes of deep CO₂ and condensate we infer that the original pre-condensed vapour has a CO₂ molar fraction of 0.08-0.12, a possible range for a geothermal vapour that, e.g., is similar to the CO₂ molar fractions measured in the nearby Amiata geothermal wells.

Finally, this work highlights, as already evidenced in previous works (e.g., Cardelliniet alii, 2003; Viveiroset alii, 2010; Biniet alii, 2019), that to obtain a reliable picture of the CO₂ spatial distribution and a reliable estimation of the CO₂ output an appropriate sampling density must be adopted in relation to the spatial structure of the CO₂ flux. In particular, it is necessary to use a measurement spacing suitable for the extent of the degassing anomalies, which can be estimated through the variographic analysis of the data (i.e., through the range of the variogram, Cardelliniet alii, 2003).

The Supplementary Material includes the datasets of the measured CO₂ flux from the soil.

We wish to thank F. Tassi and an anonymous reviewer for the helpful comments and suggestions which improved the quality of the manuscript. This study was financially supported by the MIUR project n. PRIN2017-2017LMNLAW “Connect4Carbon”.