Monitoring diffuse volcanic degassing during volcanic unrests: the case of Campi Flegrei (Italy)

In volcanoes with active hydrothermal systems, diffuse CO2 degassing may constitute the primary mode of volcanic degassing. The monitoring of CO2 emissions can provide important clues in understanding the evolution of volcanic activity especially at calderas where the interpretation of unrest signals is often complex. Here, we report eighteen years of CO2 fluxes from the soil at Solfatara of Pozzuoli, located in the restless Campi Flegrei caldera. The entire dataset, one of the largest of diffuse CO2 degassing ever produced, is made available for the scientific community. We show that, from 2003 to 2016, the area releasing deep-sourced CO2 tripled its extent. This expansion was accompanied by an increase of the background CO2 flux, over most of the surveyed area (1.4 km2), with increased contributions from non-biogenic source. Concurrently, the amount of diffusively released CO2 increased up to values typical of persistently degassing active volcanoes (up to 3000 t d−1). These variations are consistent with the increase in the flux of magmatic fluids injected into the hydrothermal system, which cause pressure increase and, in turn, condensation within the vapor plume feeding the Solfatara emission.

the location of all the CO 2 flux measurements (yellow dots) and, as example, the locations of CO 2 measurements of January 2016 survey (blue dots); the location of Bocca Nuova (BN), Bocca Grande (BG) and Pisciarelli (Pi) fumaroles; the main tectonic structures 60 ; the area considered for the mapping of CO 2 fluxes (white box); the area considered for the computation of the total CO 2 output from Pisciarelli area (PIS, box indicated by the dashed yellow line). Coordinates are reported as meters projection UTM European Datum 50. All the maps were realized with the software Surfer, Version 11.0.642 (http://www.goldensoftware.com/products/surfer). subterranean vapor plume, is returned by TOUGH2 modelling of the hydrothermal system feeding the Solfatara fumarolic field 3,4 (Fig. 2).
The emitted CO 2 is thought to derive mainly from magma degassing 34 , even if we cannot exclude a minor contribution from decarbonation of hydrothermal calcite 48 . A relatively positive (−1.3‰ ± 0.4‰ ref. 48) carbon isotope signature of the fumarolic CO 2 , as well of the CO 2 involved in the past deposition of hydrothermal calcites 48 indicates a primary origin of the CO 2 from a mantle metasomatised by crustal fluids 34,49,50 .
In this work, the results of 30 diffuse CO 2 flux surveys performed at Solfatara from 1998 to 2016 are presented and discussed. The CO 2 soil fluxes were measured over an area of ~1.2 × 1.2 km, including the Solfatara crater and the hydrothermal site of Pisciarelli (Fig. 1), using the AC (see Method). Each survey consisted of a number of CO 2 flux measurements varying from 372 to 583 (Table 1) resulting in a total of 13,158 measurements.
This data set, entirely reported in Supplementary Information Dataset S1, is one of the largest datasets made anywhere 18, 51, 52 on a single degassing volcanic-hydrothermal system. It is particularly relevant in the framework of volcanological sciences because it was acquired during a long period of unrest at CFc. Aside from making this large data set available, our main aim is to investigate how CO 2 emissions varied during the progress of the CFc volcanic unrest that is characterised by accelerating geophysical and geochemical signals 4,25,53,54 . Since the 1950's, CFc underwent several episodes of ground uplift and deformation 55 , generally accompanied by seismic swarms and abruptly followed by significant changes in the composition of fumaroles 4 . The largest bradyseismic episode occurred from 1982 to 1984 with a total uplift of 1.79 m. After about twenty years of prevailing subsidence, a new unrest phase started at the beginning of the new millennium. This crisis, still ongoing, deviates from the previous episodes due to the long duration and the clear acceleration of signals, which were recently interpreted as being mainly caused by an ongoing heating of the system 3,25 , and, regarding the 2012-2013 accelerated ground uplift, by the intrusion of magma at shallow depths 53,56 .
In the following discussion, measured CO 2 fluxes will be compared with the results of a recently published model 3 that simulated the effects of the injection of magmatic fluids into a virtual hydrothermal system, ideally representing the system feeding CO 2 emissions at Solfatara. With respect to previous physical models of the system, the injected magmatic fluids in ref. 3 are progressively richer in water, thus explaining the heating of the system. Here, for the first time, the modelled CO 2 fluxes are compared with those observed during the ongoing crisis at CFc.
Finally, a further objective of the work is the comparison of the long series of CO 2 flux data from Solfatara, a hydrothermally active volcano, with both measured geothermal systems, and CO 2 fluxes from active volcanic plumes on which global volcanic CO 2 emissions are largely based.

Results
Statistical distribution of CO 2 flux. The measured CO 2 flux of each survey is distributed in a wide range of values from >0.4 g m −2 d −1 up to 72,000 g m −2 d −1 ( Table 1). The logarithmic probability plot of Fig. 3 reveals that the CO 2 fluxes of each survey plot as a curve with an inflection point. Such curves correspond to a bimodal statistical distribution, resulting from the combination of two CO 2 flux log-normal populations (see Methods), which could indicate the occurrence of two CO 2 flux sources 29,57 , different mechanisms of gas transport, different permeability of the soil, etc.
The statistical distribution of CO 2 flux of each survey was modelled (see Methods) with the combination of a population characterised by a high mean CO 2 flux value, HF population, and one characterised by a lower mean CO 2 flux value, LF population (Fig. 3).
The mean CO 2 flux value of LF and HF populations range from 14 g m −2 d −1 to 135 g m −2 d −1 and from 1,270 g m −2 d −1 to 6,580 g m −2 d −1 respectively ( Table 1). The high mean CO 2 flux values of HF populations clearly indicate that they are fed by the CO 2 up-rising from the underlying hydrothermal system (see for example refs 15, 28, 29, 57 and 58). The mean flux value of LF populations varies within a range that precludes the possibility that they have a purely biologic source for the entire period. In fact, the mean biogenic CO 2 flux from a wide variety of ecosystems ranges from 0.2 g m −2 d −1 to 21 g m −2 d −1 (ref. 57 and reference therein). Roughly in agreement with these typical biogenic-derived fluxes, the coupled analyses of the CO 2 flux and the isotopic composition of the CO 2 efflux performed at Solfatara in March 2007 (see Supplementary Information Fig. S1) indicated a mean biogenic CO 2 flux of 26 (±3) g m −2 d −1 (ref. 58). Since 2004, the mean CO 2 flux of LF populations increased to 136 g m −2 d −1 (Table 1), i.e. values much higher than those associated with a purely biological source. The seasonal variability of the biologic production of CO 2 within the soil, resulting in soil CO 2 flux lower in autumn-winter and higher in spring-summer seasons (see for example ref. 59), cannot account for the temporal variation of the mean CO 2 flux of LF population (Table 1). Therefore, even though mean CO 2 flux values of LF populations are consistently 1 to 2 orders of magnitude lower than mean CO 2 flux values of HF populations, since 2004, the LF population clearly represents a mixture of biogenic and deeply-derived CO 2 . This suggests that deep-sourced CO 2 degassing widely affects both Solfatara crater and its surroundings.
Mapping of diffuse degassing and total CO 2 output. To characterise the spatial distribution of CO 2 fluxes, all 13,158 CO 2 flux measurements from 1998 to 2016 have been used to produce a map of CO 2 soil degassing using the sGs approach (see Methods). The result of the sGs simulations is reported as a probability map (Fig. 4) created using the threshold value for the biological background CO 2 flux as a cut-off. The threshold value selected was 50 g m −2 d −1 , the 95 th percentile of the biogenic soil CO 2 fluxes defined on the basis of the isotopic compositions of the CO 2 efflux of the March 2007 survey 58 .
The map shows the presence of a well-defined diffuse degassing structure (Solfatara DDS, Fig. 4), that is the area characterised by degassing of deeply derived CO 2 (see Methods). Solfatara DDS is defined in yellow to red in Fig. 4 Table 1. Summary statistics of the CO 2 flux dataset, statistical parameters of the partitioned CO 2 flux populations, areal extents of the Solfatara DDS and total CO 2 output estimates. The measured CO 2 flux mean are presented as mean ± SD. The mean CO 2 flux of statistically partitioned population are presented as mean ± SD of the mean. The total CO 2 output is presented as mean ± SD.
As this map was created considering all the measurements, it reasonably highlights the areas that have been affected more frequently by the degassing of hydrothermal CO 2 during the last eighteen years.
The Solfatara DDS includes both the area inside the Solfatara crater and neighbouring areas outside the crater, in particular the Pisciarelli area to the E, the Monte Olibano to the S and the area of via Antiniana to the SE (Fig. 4). The shape of the Solfatara DDS is well correlated with the location of volcanic and extensional tectonic structures (faults and fractures, Fig. 4) that allow the gas to transfer towards the surface. The Solfatara DDS is bounded to the NW, and interrupted along a NW-SE band between Solfatara crater and Pisciarelli, by low flux areas corresponding to the outcrop of volcanic products belonging to Astroni deposits 60 that can act as a permeability barrier for the rising deep CO 2 . In fact, Astroni deposits are locally constituted by massive and fall out and ash surge deposits with thickness from ~10 to ~30 m, affected only by mesoscale normal faults with lengths of tens of centimeters and displacements of a few centimeters 60 .
However, where normal faults deform and cut the Astroni deposits with up to metric dip separation 61 (e.g., at Via Antiniana area and Pisciarelli), high CO 2 fluxes occur at the surface (Fig. 4), suggesting that anomalous CO 2 pressures can be present below this low permeability layer. The roughly NE-SW-elongated CO 2 flux anomalies in the northern part of Pisciarelli further support the probable continuity of the CO 2 anomaly below the Astroni deposits. These anomalies match the directions of the drainage network that eroded the Astroni deposits in the eastern flank of the Solfatara cone, resulting in higher fluxes along the valleys. Finally, the geometry of Solfatara DDS is also affected, especially at the Pisciarelli and via Antiniana areas, by intense urbanization (e.g., presence of roads, buildings, excavations and paved squares). In particular, an evident low-flux zone in the via Antiniana area (area A in the map of Fig. 4) coincides with the presence of a paved terrace, while high flux zones in the southern part of Pisciarelli (area B in the map of Fig. 4) correspond to recent excavated zones. This behaviour suggests the high impact of anthropogenic intervention in the natural degassing zones.
To investigate the changes over the time of the spatial distribution of CO 2 fluxes, and of the total amount of CO 2 released by diffuse degassing, each dataset was processed by the sGs method. The results are reported in Fig. 5 as probability maps.
The reliability of the produced maps is supported by the good spatial continuity and well-defined spatial structure of all the datasets indicated by experimental variograms of the CO 2 flux n-score (see Methods, Supplementary Information Fig. S2 and Table S1).
The areal extent of the DDS (1) was computed for each of the maps relating to different surveys (Fig. 5), as the area where the probability that the CO 2 flux is greater than 50 g m −2 d −1 is higher than 50% (see Methods). The DDS extent ranges from 0.45 × 10 6 m 2 in 1998 to more than 1 × 10 6 m 2 in many surveys after 2012. The expansion of DDS mainly interests the areas external to the Solfatara crater, that before 2003 were characterised by low, biogenic, CO 2 fluxes (Fig. 5). The largest expansion occurred in the Pisciarelli area and particularly in correspondence with the NE-SW fault network of Pisciarelli area and along a band connecting Pisciarelli with the degassing area of via Antiniana in the south. Examples of variations occurred over time in different areas are reported in Supplementary Information Fig. S3. In these areas (areas 1, 2 and 3 in Supplementary Information Fig. S3) the median of the CO 2 flux passed from typical background values (10-40 g m −2 d −1 ) in the pre-2003 period to values higher up to 1 order of magnitude. The 2003 increase of CO 2 fluxes and the DDS enlargement were already interpreted as due to the "first arrival" of hydrothermal CO 2 at the surface, in the peripheral areas of Solfatara 62 . The total amount of CO 2 released by diffuse degassing (diffuse total output, DTO) was estimated from the results of the sGs (see Methods) for the different surveys and ranges from 745 (±61) t d −1 to 2,815 (±318) t d −1 (Table 1). Even if the DTO refers to the sum of all the CO 2 sources active in the area, it provides a good estimate of the hydrothermal CO 2 release because CO 2 fluxes from the hydrothermal source are much higher than those from biogenic sources. For example, assuming a constant biogenic background CO 2 flux of 26 g m −2 d −1 over the entire surveyed area, the biogenic CO 2 output would result in ~38 t d −1 . Considering this biogenic CO 2 output for all the surveys, it is only from 1% to 5% of DTO and always within the DTO uncertainty (from 8% to 13%, Table 1). In the following, we will consider the DTO as a representative of the deep hydrothermal source, also considering that the low contribution of biogenic sources is likely an overestimation, as portions of the surveyed area are characterised by the absence of vegetation.

Discussion
In this section, variations in the CO 2 degassing during the investigated period are discussed and compared with other geochemical and geophysical signals. We refer in particular to the variations that affected, in one aspect, the "background" CO 2 emission (LF population) and the extent of DDS ( Fig. 6a and b), and in another, the total CO 2 output (Fig. 7a). The mean CO 2 flux of the LF population, which will be referred here as SBF (Solfatara background flux), increases from 20-30 g m −2 d −1 in 1998 to 136 g m −2 d −1 in the last campaign of June 2016, following a scattered, but almost continually increasing trend (Fig. 6a). The maps report the probability that the simulated CO 2 flux is greater than 50 g m −2 d −1 , selected as the threshold for a pure biogenic CO 2 flux. Yellow to red colours i.e., probability of CO 2 Flux >50 g m −2 d −1 is higher than 50% define the Solfatara DDS where degassing deeply derived CO 2 occurs. Tectonic structures are from ref. 60. The maps were realized with the software Surfer, Version 11.0.642 (http://www.goldensoftware.com/products/surfer).
It is worth noting that typical mean CO 2 fluxes generated by a biogenic source in the soil (20-30 g m −2 d −1 , see above) were only measured in the first three campaigns, whilst from 2003 onwards, SBF were 2-5 times higher. This increase implies that after 2003, the SBF began to represent a mixture of biogenic and deeply-derived CO 2 , indicating that areas previously unaffected by an anomalous CO 2 degassing started to release deeply-derived   12,31,32,65 , is reported together with the estimated CO 2 release by diffuse degassing. The diffuse degassing computed for the Pisciarelli area (area PIS in Fig. 1) is also reported. The CO 2 output simulated by the physical-numerical simulation is reported with a 0.7 scaling factor. (b) Earthquake magnitudes, number of earthquakes and ground elevation at RITE GPS station are reported for comparison 4,54 (2016 data update from http://www.ov.ingv.it/ov/it/banchedati.html). (c) Comparison of total CO 2 output with the P CO2 estimations based on two recently published alternative approaches to gas equilibria 25,67 (see the text for explanations).
Scientific RepoRts | 7: 6757 | DOI:10.1038/s41598-017-06941-2 gas. In fact, the SBF increase is correlated with an important change of the spatial pattern of CO 2 fluxes as shown by maps of the Solfatara DDS (Fig. 5). Chronograms in Fig. 6 show an initial growth of Solfatara DDS and an increase of SBF from February 2003 to March 2004, when the DDS doubled its extent from ~0.5 km 2 to ~1 km 2 ( Table 1 and Fig. 6b) and SBF increased from 20-30 g m −2 d −1 to 40-50 g m −2 d −1 . After this expansion, the DDS extent and the SBF remained quite stable until July-October 2012 when, after a relative reduction, DDS reached an extent of ~1.2 km 2 ( Table 1) and SBF increased to 70-100 g m −2 d −1 . We observe that the DDS extent reached 80-90% of the 1.4 km 2 of the investigated area, implying that this parameter is now close to saturation and further important variation in the future are not possible while SBF could continue to increase.
The SBF increase and the extension of Solfatara DDS may reflect variations in the vapor plume (Fig. 2) feeding the Solfatara and Pisciarelli emissions. According to ref. 25, steam condensation and temperature increase affects the hydrothermal system of Solfatara. This process could lead to a generalised increase of fluid pressure and to the formation of batches of condensates within the vapor plume 25 . Repeated episodes of mud emission 3,25 at Pisciarelli, occurred in April 2006, July 2012, October 2013, July 2014, February 2015, May 2016 and February 2017, confirm the increased amount of condensates produced by the system. It is worth noting that the enlargement of the Solfatara DDS and the increase of SBF proceed concurrently with the increase in the fraction f of condensation estimated from the composition of the main fumaroles of Solfatara ( Fig. 6a and b, see Methods). Therefore, it is highly probable that the SBF increase, as well as the enlargement of Solfatara DDS, is linked to the deep dynamics of the vapor plume and in particular to the ongoing condensation and heating processes.
Despite the good agreement between variations of DDS, SBF and f values, deviations from the long-term trend depicted in Fig. 6 can be caused by a short-term process and by the uncertainty of f estimation. For example, the DDS extent could be partially controlled by the rain that affect CO 2 fluxes increasing soil water content, changing soil permeability, dissolving soil CO 2 , etc. These processes can affect the DDS extent because are more effective in low-flux areas, that are those mainly contributing DDS extent variations.
Other relevant changes regarded the DTO (Fig. 7), that is well representative of deep source degassing. The increase of DTO is mainly due to the increase of the degassing rate in areas of high CO 2 flux that reached a peak in 2014-2015. For example, considering the measured fluxes inside the Solfatara crater and in the fumarolic area southern of BG-BN flumaroles (areas 4 and 5 in Supplementary Information Fig. S3) the median values of CO 2 flux increased from ~500 g m −2 d −1 and ~400 g m −2 d −1 during 2003-2013 to ~800 g m −2 d −1 and ~1,100 g m −2 d −1 after 2013, respectively ( Supplementary Information Fig. S3). The most important increase occurred in Pisciarelli area where the CO 2 flux passed progressively from ~100 g m −2 d −1 to ~300 with a peak value of ~1,000 g m −2 d −1 in 2015.
The relative increase of CO 2 fluxes is particularly relevant at Pisciarelli, where the CO 2 diffuse emission rose from 90 t d −1 to 260 t d −1 from 2003 to 2016, with a large peak of 428 t d −1 in 2015 (area PIS in Fig. 1). Furthermore, at Pisciarelli, the increase of diffuse emission corresponds to an evident macroscopic intensification of the hydrothermal activity 25,53 .
The DTO represents only a fraction of the total CO 2 output, as it does not include the contribution from fumarolic vents. The vent contribution to the Solfatara CO 2 budget most probably increased in the last years concurrently with the abovementioned increase in the fumarolic activity. Measurements of the CO 2 flux performed in the 1980's (refs 63 and 64) suggest, in fact, that the CO 2 emission from fumaroles was relatively low in the past, also during the large crisis of 1982-1984. On the contrary, more recent measurements (2012-2015 period) performed with different techniques point to significant CO 2 emissions from the fumarolic vents that, from BN-BG and Pisciarelli fumaroles, were estimated to be 343 to 858 t d −1 (refs 12, 31, 32 and 65). It is also worth noting that the CO 2 flux from the vents of Pisciarelli area, that was evaluated at only ~18 t d In Fig. 7b the CO 2 release is compared with the seismicity 4 and ground deformation 54 . Since 2005, the occurrence of earthquakes (located at depth <3 km, ref. 66) increased together with the beginning of an uplift phase. Since 2005, a general correlation emerges between the increase of CO 2 release and the increase in seismicity and ground uplift, although the gas flux time series is more discontinuous than geophysical ground uplift and seismicity records (Fig. 7b). In particular, in the period from 2012 to 2016, two main peaks in the CO 2 output coincide with episodes of accelerating ground uplift and of intensification in the seismicity. It is also worth noting, that these two periods of stronger degassing are separated by a period of relatively low CO 2 output, which corresponds to a period of low seismicity and a pause in the deformation.
These correlations support the link between the degassing and the deformation-seismicity, both controlled by repeated episodes of magmatic fluids input into the hydrothermal system, episodes that increased in frequency after 2005 (refs 3, 4 and 25). These episodes of magmatic fluid injections, modelled with a geothermal simulator, involved increasing amounts of fluids since the mid 2000s 3, 4 . Here we show, for the first time, the modelled CO 2 output returned by the last published simulation of the process 3 . The modelled CO 2 output from the entire simulation domain (Fig. 2) shows roughly the same increasing trend of the measurements (Fig. 7a), even if the observed CO 2 fluxes are systematically lower than the modelled CO 2 output (~70%). Considering the approximations of the model, this discrepancy is very low. It is possibly caused by either a general overestimation in the modelling, or the fact that our measurements do not include the entire amount of the released CO 2 . For example, the measurements neither include the contribution of the many low flux fumaroles (which are not measured and not measurable with the AC, see Method), nor the CO 2 emissions from zones outside the investigated area. In any case, the match between modelled emission and measurements provides convincing evidence that the visible escalation in degassing activity is caused by repeated injections of magmatic gases into the Campi Flegrei hydrothermal system as simulated in the modelling 3 .
Recently, two alternative approaches to gas equilibria leaded to opposite results regarding the temporal evolution of T-P conditions of the hydrothermal system feeding the Solfatara emission. Ref. 25 points to a generalized increase of pressures and temperatures while the other interpretation indicates the progressive depressurization of the system 67 , a process that would be still ongoing. This depressurization is clearly in contrast with the evolution of total CO 2 output that almost tripled from early 2000's in agreement with the increase of the pressure in the hydrothermal system as returned by the model of ref. 25 (Fig. 7c).
The CO 2 total output in recent years from the Solfatara hydrothermal system can be reasonably assumed to be at least 2000-3000 t d −1 . For comparison, the estimations for the last period are much higher than the CO 2 emitted by geothermal power plants around the world. They are, for example, 4-6 times higher than those emitted by Icelandic geothermal power plants (435 t d −1 , ref. 68) and ~3 times more than those associated with geothermal plants in New Zealand (~800 t d −1 ; www.nz.geothermal.org.nz/emissions/).
In Fig. 8, the CO 2 flux released by Solfatara is compared with the mean volcanic plume CO 2 fluxes from persistently degassing volcanoes reviewed by ref. 69. Considering only the contributions from soil diffuse degassing, it emerges that Solfatara DDS on average (i.e., mean DTO = 1,309 t d −1 ) sustains a daily CO 2 flux to the atmosphere similar to a "medium-large" volcanic plume. If, instead, we consider the highest total CO 2 release (vents + diffuse) of 3,420 t d −1 , measured in January 2015, this would constitute the eighth highest value ever measured at volcanic plumes.
This finding queries the reliability of the actual estimates of the natural flux of CO 2 from volcanic activity, considering that many calderas around the world which are affected by hydrothermal sites, similar to Campi Flegrei, are not normally included in most of the global estimates. We think that the flux of CO 2 from hydrothermal sites is potentially globally relevant considering also that can reach very high values up to the 10-60 kt d −1 estimated for Yellowstone 17,70 .
The long-term volcanic CO 2 degassing at Solfatara, and degassing from hydrothermal systems in general, can contribute to refining estimates of the volcanic CO 2 contribution to the atmosphere, and can aid the ability to assess the possible role of CO 2 degassing from hydrothermal systems. Methods CO 2 flux measurement. CO 2 flux measurements were performed using the accumulation chamber method, which is based on the measurement of the CO 2 concentrations over time, inside an inverted chamber placed on the ground [27][28][29] . The used instruments were developed, assembled and tested at the laboratories of Università di Perugia and INGV of Naples. Each instrument consists of: (1) a metal cylindrical vessel (the accumulation chamber, AC), (2) an Infra-Red (IR) spectrophotometer, (3) an analog-digital (AD) converter, and (4) a palmtop computer ( Supplementary Information Fig. S4). Since 2003, the measuring apparatuses have been equipped with LI-COR IR sensors, Li-800 or LICOR Li-820, operating in the range 0-20,000 ppm of CO 2 . The instrument used for the surveys of 1998 and 2000 was equipped with a Dräeger Polytron IR sensor with adjustable measurement range from 2,000 ppm up to 100% vol. The gas is circulated from the AC, of ~2.8 L volume, to the IR and vice versa by a pump (~0.0167 L s −1 ). The AC is internally equipped with ring-shaped perforated manifold that re-injects the gas favouring the mixing in the chamber. The CO 2 concentration in the circulating gas is acquired every 0.025 s and transmitted to a palmtop computer where is plotted versus the measuring time. A specific applicative (Gasdroide, https://bitbucket.org/moovida/gasdroide), capable of acquiring the digital signals from the AD and to elaborate concentration vs time diagram, was developed in 2012 for the Android operating system and released with a General Public License.
The CO 2 flux is computed in real time from the rate of CO 2 concentration increase in the chamber (dC CO2 /dt) according to the relation CO 2 flux = cf × dC CO2 /dt (ref. 28). The proportionality factor (cf) between dC CO2 /dt and the CO 2 flux was determined before each survey by laboratory tests, during which imposed CO 2 fluxes typically in the range 10 to 10,000 g m −2 d −1 , were measured over a "synthetic soil" made of dry sand (10 cm-thick) placed inside a plastic box. In each calibration test, the cf value was computed from the linear best-fit line of CO 2 flux vs. dC CO2 /dt.
The measured soil CO 2 flux and the measurement locations of each survey are reported in Supplementary Information Dataset S1, together with the soil temperature at ~10 cm depth measured at the same time and at the same site of CO 2 flux measurement. The soil temperature data is provide for completeness but is not discussed in this work.
Statistical data elaboration. The CO 2 flux data was analysed by statistical methods to define the statistical parameters of the flux, which offer insight into the origin of degassed CO 2 .
In volcanic-hydrothermal areas, the CO 2 soil diffuse degassing is frequently fed by multiple gas sources, such as biological and volcanic (ref. 29 and references therein). The multiple origin of the gas can result in a bimodal statistical distribution of CO 2 flux values, which plots as a curve with an inflection point on a logarithmic probability plot (see Fig. 3). In fact, while a single log-normal population plots as a "straight" line on a logarithmic probability plot, a polimodal distribution resulting from the overlapping of n log-normal populations plots as a curve with n-1 inflection points 71 . The partition of such complex statistical distributions into individual log-normal populations and the estimation of the proportion (f i ), the mean (M i ), and the standard deviation of each population were performed following the graphical-statistical procedure proposed by ref. 71. which is largely applied to soil CO 2 fluxes 28,29,57 . Since the so computed M i value refers to the logarithm of the CO 2 flux values, the mean value of CO 2 flux was estimated using a Montecarlo simulation procedure ( Table 1).
The estimated mean CO 2 flux values have been used in literature to compute the CO 2 output pertinent to each population associating a fraction of the degassing area (i.e., S i = f i S where S is the total extent of the surveyed area, GSA approach) 28 . However, the reliability of the CO 2 output obtained from the GSA can be affected by arbitrary choices in the partitioning procedure 29 . In particular, the interpretations of the tails of the distributions are critical as, especially at the high flux values, they are generally defined based on a low number of measured values. Furthermore, the GSA approach does not consider the spatial distribution of the measurements implicitly assuming a homogeneous sampling density.
Geostatistical data elaboration. In order to obtain a reliable estimation of the total CO 2 output and to produce maps of the CO 2 flux, the CO 2 fluxes were elaborated with a geostatistical approach proposed by ref. 29 based on sequential Gaussian simulations (sGs). According to refs 29 and 72, sGs yields a realistic representation of the spatial distribution of CO 2 fluxes reproducing both the statistics and the spatial features of the experimental data.
The sGs method consists of the production of numerous equiprobable realizations of the spatial distribution of an attribute (i.e., maps of CO 2 flux in this study), here performed using the sgsim algorithm of the GSLIB software library 73 . Since the sGs assumes multi-Gaussian distribution of the attribute to be simulated, the CO 2 flux values were transformed into a normal distribution (n-scores of data) using the nscore algorithm of the GSLIB software library. The n-scores are then used in the simulation procedure and transformed back into flux values at the end of the simulation process, using the inverse of the normal score transform. The CO 2 flux values are simulated at locations defined by a regular grid, here consisting in a grid of 14,520 squared cells (121 cells in the EW and 120 cells in the NS direction) with a cell size of 10 × 10 m (except for the map reported in Fig. 2 for which a 2 × 2 m simulation cell was used). The n-scores are simulated by random sampling of a Gaussian cumulative distribution function (cdf), defined at each location on the basis of a mean value and variance computed at each grid node by means of simple kriging method. Simple kriging estimate and variance are computed considering the measured data and those previously simulated during the procedure, according to the variogram model of n-scores and to the statistical distribution of the data. The variogram model is defined fitting the experimental variogram of n-scores and provides a description of how the data are spatially correlated. The variogram models are given in terms of nugget, range and sill parameters, where the nugget represents the small-scale variation (including measurement errors), 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 simulation was run in order to produce 100 realizations for each dataset.
The produced realizations were post-processed to produce the E-type estimate map and the probability map. The E-type estimate map, i.e. the map of the "expected" value at any location, is obtained through a pointwise linear average of all the realizations. The probability map consists in a map of the probability that, among all the realizations, the simulated CO 2 flux at any location (i.e., at grid nodes) is above a cut-off value. The probability map drawn for each survey is reported in Fig. 5 while the maps of E-type estimate are reported in the Supplementary  Information Fig. S5.
According to ref. 29, selecting the threshold value of biogenic CO 2 flux as a cut-off, the probability maps were used to define the extent of the diffuse degassing structures (DDS) 28 , that is the area interested by the release of deeply derived CO 2 . According to ref. 29 in this work, the DDS is considered as the area where the probability, among the 100 realizations, that the simulated CO 2 flux is higher than the biogenic CO 2 flux threshold over 50%.
The simulated flux values by sGs were used also to compute the total CO 2 release by diffuse degassing. The total CO 2 release is computed for each realization by summing the products of the simulated CO 2 flux value at each grid cell by the cell surface. The mean of the values of total CO 2 release computed for the 100 realizations are assumed as the characteristic diffuse total CO 2 output (DTO, see above) for each period. The standard deviation of the 100 estimates is assumed as the DTO uncertainty.
Fraction of condensed steam. A previous work by ref. 25 illustrated several evidences of heating of the hydrothermal system feeding the Solfatara emissions based on the compositional variations of the main fumaroles of the area (BG and BN). The pressurization of the vapor plume (Fig. 2) The variable X CO2,fu and X CO2,fu are the analytical values, while P H2O,eq and P CO2,eq refer to the equilibrium values computed considering gas equilibria in the H 2 O-H 2 -CO 2 -CO gas system. In detail, we used the following expressions 38 where the partial pressures of water and CO 2 are assumed equal to their fugacities. Equation (9) was used to compute the condensed fraction f for each BG and BN fumarolic samples reported in Fig. 6. The fumarolic compositions are available in the literature 3 .