Bagni San Filippo area is characterized by the discharge of thermal waters and deeply produced CO2-rich gases both from vents and soil diffuse degassing. The thermal waters are the results of the mixing between meteoric waters and hot fluids deriving from the condensation, at depth, of vapours uprising from a deep hydrothermal reservoir. This process gives rise to a relatively shallow thermal system at temperature close to 50°C, characterized by SO4-rich and Cl-poor waters and elevated PCO2 (~7 bar). Most of the incondensable gas of deep originated vapour is released as a free gas phase forming cold gas vents and localized spots of anomalous CO2 diffuse degassing. The location and the shape of these degassing zones are strongly controlled by the main tectonic structures of the area. Through detailed soil diffuse degassing surveys and hydrogeochemical modelling, we estimate at 226-326 t d-1 and at 965 t d-1 the deep CO2 emission and the amount of condensates discharged by the thermal springs, respectively. The thermal energy associated to the process results at ~ 29 MW, most of which (~ 25 MW) is associated with condensation occurring at depths greater than groundwater circulation.

Central-southern Italy is affected by an active and intense process of CO2 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 CO2 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 CO2 degassing in central and southern Italy, showing the occurrence of two large scale CO2 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×104 t d-1 (i.e., 6.16 Mt y-1) and 8.44×103 t d-1 (i.e., 3.08 Mt y-1) of CO2, respectively. The total amount of CO2 released to the atmosphere through diffuse regional degassing (9.24 Mt y-1) is the main component of the geological CO2 budget in Italy, being higher than the amount of CO2 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 CO2-rich groundwater, both TRDS and CDS are also characterized by the presence of many CO2 vents and areas of anomalous soil CO2 diffuse degassing (e.g., Chiodiniet alii, 1999, 2007, 2011; Rogieet alii, 2000, 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 CO2 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 CO2 degassing, CO2-rich thermal waters with temperatures from 23 to 48°C and large travertine deposits (Frondiniet alii, 2009). The areas of stronger CO2 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; LF302.html; Chiodiniet alii, 2008).

In this work, we combine i) the results of a detailed survey on the CO2 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 CO2 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 Mwe (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 CO2, chemical and isotopic composition of gas discharged, chemical and isotopic composition of groundwater and chemical composition of the gas dissolved by groundwater.

Soil diffuse degassing

From December 2003 to May 2004 four CO2 flux surveys were performed at BSF area using the accumulation chamber method (Chiodiniet alii, 1998), using a LI-800 infrared spectrometer as CO2 detector operating in the range 0-20000 ppm. During December 2003, an area of about 1.7 km2 (Fig. 1b) was investigated with 595 measurements of CO2 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 CO2 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 CO2 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 CO2 fluxes are reported in the Supplementary Material.

The CO2 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 CO2 flux and the sources of the diffuse degassing. In fact, in volcanic and geothermal areas CO2 soil diffuse degassing is frequently fed by multiple gas sources, such as the biological and the deeply produced CO2. This multiple origin can result in a polymodal distribution of the CO2 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 (fi) between different populations, their mean (Mi) and standard deviation (σi). Since the calculated Mi refers to the logarithm of CO2 flux values, the mean of the CO2 flux was then estimated using a Monte Carlo simulation procedure. The uncertainty of this estimate is given in terms of 5th and 95th percentiles of the mean CO2 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 CO2 degassing and to quantify the total CO2 output when possible. The sGs method consists of the production of numerous equiprobable realizations of the spatial distribution of the CO2 flux, here performed using the sgsim algorithm of the GSLIB software library (Deutsch & Journel, 1998). The CO2 flux is simulated on a regular grid, according to the variogram model defined fitting the experimental variogram of the normal scores of the CO2 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 CO2 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 CO2 flux at any location (i.e., at grid nodes) is above a cut-off value. Selecting the threshold value of the biogenic CO2 flux as a cut-off, the probability map is used to highlight the area interested by the emission of deeply derived CO2 (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 CO2 flux is higher than the biogenic CO2 flux threshold is over 50%.

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

Groundwater and gas emissions sampling and analyses

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 HCO3 were measured in the field. HCO3 concentration was determined by acid titration with 0.01N HCl using methyl orange as indicator. Water samples for chemical and isotopic (δD and δ18O) 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 SrCO3 by adding SrCl2 and NaOH in solid state to the water sample. Carbonate precipitate was recovered by filtering the sample and washed with distilled water in a CO2-free atmosphere in the laboratory. Carbon isotopic analyses of the TDIC (δ13CTDIC) were performed on the CO2 gas liberated from the SrCO3 precipitate by reaction with 100% H3PO4.

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 SO4 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‰, δ18O ± 0.08‰ and δ13C ± 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 H2O2, by acid-base titration and by ion-chromatography, respectively. The gas species not absorbed by the NaOH solution (i.e., He, Ar, O2, N2, H2, and CH4) 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 CO2 was determined on 50 mL dry gas sample by the same instrumentation used for δ13CTDIC analyses

CO2 diffuse degassing

General survey (GS)

The measured CO2 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 CO2 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 CO2 flux is the result of the overlapping of 2 log-normal CO2 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 CO2 deep source. On contrary the mean flux of LFP-GS well falls within the range of the CO2 fluxes produced by the biological activity in the soil in a wide variety of ecosystems (mean CO2 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 CO2 flux for the period of the GS.

The variogram of the normal scores of the CO2 flux of the GS shows a high nugget value (0.75, Fig. 4a) suggesting that the CO2 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 CO2 output of 849 t d-1 was computed for the GS survey (Table 3) by subtracting to the total CO2 output estimated by sGs approach (865 t d-1, see Methods) the biological CO2 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 CO2 fluxes, in comparison to the spacing between the measurements. As already noted in other areas, a large measurement spacing of highly spatially variable soil CO2 fluxes can in fact cause an overestimation of the total CO2 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 CO2 flux probability map (Fig. 2) which has been drawn assuming a CO2 flux cut-off value of 25 g m-2 d-1 (i.e., the 95th percentile of the LFP-GS population, Fig. 3a).

Detailed surveys (BO, PA, LA)

The measured CO2 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 CO2 fluxes is interpreted as the results of the overlapping of three CO2 flux populations (Fig. 3b-d). The populations with the lowest mean CO2 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 CO2 contribution that in any area resulted higher than the biological background CO2 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 CO2 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 CO2 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 CO2 flux too high for the biogenic source (Table 2). The presence of two deeply derived CO2 flux populations could reflect the relative prevalence of diffusive or advective CO2 transport mechanisms.

The computation of the CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 CO2 degassing is evident in the southern part of the area and along the NNE-SSW direction are located numerous CO2 vents and a large morphological depression (Fig. 5a) where a strong CO2 vent at its bottom causes high CO2 air concentrations (up to 50 % vol CO2 was measured during the survey). At Lavinate the CO2 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 CO2 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 CO2 degassing suggest the main role of the tectonic structures on the upward transfer of the deeply produced CO2.

Based on the sGs results, the total CO2 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 CO2 output was then estimated subtracting the biological CO2 contribution (i.e., the mean CO2 flux of the respective LFP populations) to the total CO2 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 CO2 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 CO2 output of 47 t d-1 is obtained for the entire Bollore area.

Total CO2 output from BSF

The total deep CO2 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 CO2 from CLV and LE results at 30 t d-1 to 80 t d-1. Therefore, our best estimation of the total deep CO2 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.

Processes controlling the chemical and isotopic compositions of thermal springs and gas emissions

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-SO4 component with a Ca–HCO3 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 SO4 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 CO2 of the geothermal vapours and the sulphate content possibly from both dissolution of anhydrite and oxidation of deeply derived H2S.

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 δ18O 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 CO2 (i.e., of condensed steam and CO2). In one hand the steam produced by boiling is depleted in 18O and 2H with respect to the parent liquid, in the other hand the CO2 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 CO2 (PCO2) 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 CO2. For the spring 33, assumed as representative of the system (see above), we estimated an original PCO2 of ~7 bar by adding CO2 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 PCO2. 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 CO2 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 CO2 dominated composition (Table 6) and relatively high contents of CH4 and nonatmospheric N2. The CH4 concentrations range from 1.4% to 1.9% while the N2 concentrations are of 2% - 2.5% with N2/Ar ratios from 2200 to 6300 (excluding the samples BSF4 and BSF10 whose N2/Ar of ~ 100 indicates an evident air contamination during sampling). The relative contents of N2, He and Ar of BSF gases fall in the range of other gas emissions of the TRDS (Fig. 9) that overall are enriched in N2 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 N2/He (and N2/Ar) ratios of the TRDS gases are ascribed to N2 enrichment due to the thermal decomposition of organic matter recycled through subduction (e.g., Chiodiniet alii, 2007).

Also the carbon isotopic composition of the CO213CCO2 from -2.8‰ to -2.3‰, Table 6) is in the typical range of TRDS gases (δ13CCO2 from -4‰ to -2‰, Chiodiniet alii, 2007) which have been interpreted as the result of a mixing between a CO2 produced at depth from decarbonation of carbonates and a more negative mantle CO2 (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 (δ13CTDIC 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 CO2, and the carbon from calcite dissolution (δ13C = +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., N2/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 PCO2 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 PCO2 values (i.e. 0.2 bar, curve B in Fig. 10). In general the measured δ13CTDIC and TDIC of BSF thermal springs follow this trend of depressurizing solutions that have reached relatively high PCO2 at depth. In particular, the model results fit the data of the main and hottest spring 33 (Fig. 10).

Conceptual model of the system and flow rates of deep CO2 and steam condensate

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 CO2 emissions. This is also the zone that was surveyed for diffuse soil CO2 fluxes (Figs.1 and 2), whose results indicate a total deep CO2 flux of 160-210 t d-1 (see section 4.1). Considering the CO2 emission from the gas vents (estimated in 50-100 t d-1, Tassiet alii, 2009) and the amount of CO2 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 CO2 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 CO2 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 CO2 molar fraction XCO2 of 0.08-0.12 (XCO2 = molCO2/(molcond + molCO2)). This is plausible XCO2 range for a geothermal vapour (e.g., XCO2 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, XCO2 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 HCO3-SO4 composition and temperatures up to 48°C, and numerous cold emissions of large amounts of CO2, both from soil diffuse degassing and vents, occur in Bagni San Filippo area, an area of few km2. 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 CO2 emissions point to a unique deep carbon source (δ13C = -2.3‰). In a first stage, the condensates and the deep-originated CO2 enter the shallow aquifer causing a temperature increase to 50°C and elevated PCO2 (~7 bar). In a second phase, the ascending thermal solutions depressurize and release a great part of the gas contributing to the observed CO2 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 CO2 gas emission, both vents and diffuse degassing structures, is strongly controlled by the tectonic lineaments. Diffuse CO2 flux surveys were performed to quantify the deep CO2 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 CO2 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 CO2 flux of 130 t d-1 from diffuse degassing. Our best estimation of the total deep CO2 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 CO2 and condensate we infer that the original pre-condensed vapour has a CO2 molar fraction of 0.08-0.12, a possible range for a geothermal vapour that, e.g., is similar to the CO2 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 CO2 spatial distribution and a reliable estimation of the CO2 output an appropriate sampling density must be adopted in relation to the spatial structure of the CO2 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 CO2 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”.