|Home | About | Journals | Submit | Contact Us | Français|
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.4km2), 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.
Volcanoes emit volatiles through active plumes, fumarolic vents and zones of diffuse soil degassing1, 2. Emitted volatiles may represent the surface manifestation of magma degassing2–6 providing useful information for the better understanding of processes occurring at depth, for assessing the state of activity of a volcano and, potentially, for forecasting the likelihood of a volcano erupting. Because of the relatively low solubility of CO2 in silicate melt7–9, CO2 is particularly useful as it exsolves from magma at greater depths than other volatile species, and therefore can reflect deep processes10–13. Diffuse CO2 degassing may represent the dominant mode of volcano degassing at calderas and volcanoes with hydrothermal activity (see for example refs 14–19), Several calderas have shown signs of unrest (ref. 20 and refs therein), however in some cases is problematic to understand if these are driven by magmatic activity (e.g., magma intrusion) or are related to hydrothermal dynamics (e.g., pressurization of the hydrothermal system)3, 12, 21–23.
A reliable technique for measurement of soil diffuse CO2 degassing (accumulation chamber, AC, see Methods) was developed at the end of 20th century, rapidly becoming extensively used in volcanological sciences26, 27. Solfatara is one of the first sites in the world where this technique, together with those used in soil CO2 diffuse degassing data analysis, were tested and improved throughout the 1990s28, 29. In general, Solfatara has become a natural laboratory for testing new types of measurements for the gas flux from hydrothermal sites based on the in situ and remote sensing determination of CO2 12, 30–33.
Hydrothermal activity at Solfatara results in numerous fumaroles and in widespread hot soils and diffuse gas emissions. The thermal energy released by diffuse degassing at Solfatara is by far the main mode of energy release from the entire Campi Flegrei caldera15.
The diffuse degassing at Solfatara is fed by a 1.5–2km-deep subterranean vapor plume, the presence of which was first hypothesised based on geochemical conceptual models of the fumaroles15, 34–40 and subsequently highlighted by the re-interpretation of seismic tomography of CFc25, 41, 42. The same concept, i.e. the presence of a subterranean vapor plume, is returned by TOUGH2 modelling of the hydrothermal system feeding the Solfatara fumarolic field3, 4 (Fig. 2).
The vapor plume connects the surface to a hydrothermal zone at about 2km depth15, 25, 37, 39, 43, 44, where the meteoric fluids mix with magmatic gases coming from a deeper zone at about 3–6km depth45–47.
The emitted CO2 is thought to derive mainly from magma degassing34, even if we cannot exclude a minor contribution from decarbonation of hydrothermal calcite48. A relatively positive (−1.3‰±0.4‰ ref. 48) carbon isotope signature of the fumarolic CO2, as well of the CO2 involved in the past deposition of hydrothermal calcites48 indicates a primary origin of the CO2 from a mantle metasomatised by crustal fluids34, 49, 50.
In this work, the results of 30 diffuse CO2 flux surveys performed at Solfatara from 1998 to 2016 are presented and discussed. The CO2 soil fluxes were measured over an area of ~1.2×1.2km, including the Solfatara crater and the hydrothermal site of Pisciarelli (Fig. 1), using the AC (see Method). Each survey consisted of a number of CO2 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 anywhere18, 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 CO2 emissions varied during the progress of the CFc volcanic unrest that is characterised by accelerating geophysical and geochemical signals4, 25, 53, 54. Since the 1950’s, CFc underwent several episodes of ground uplift and deformation55, generally accompanied by seismic swarms and abruptly followed by significant changes in the composition of fumaroles4. The largest bradyseismic episode occurred from 1982 to 1984 with a total uplift of 1.79m. 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 system3, 25, and, regarding the 2012–2013 accelerated ground uplift, by the intrusion of magma at shallow depths53, 56.
In the following discussion, measured CO2 fluxes will be compared with the results of a recently published model3 that simulated the effects of the injection of magmatic fluids into a virtual hydrothermal system, ideally representing the system feeding CO2 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 CO2 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 CO2 flux data from Solfatara, a hydrothermally active volcano, with both measured geothermal systems, and CO2 fluxes from active volcanic plumes on which global volcanic CO2 emissions are largely based.
The measured CO2 flux of each survey is distributed in a wide range of values from >0.4gm−2 d−1 up to 72,000gm−2 d−1 (Table 1). The logarithmic probability plot of Fig. 3 reveals that the CO2 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 CO2 flux log-normal populations (see Methods), which could indicate the occurrence of two CO2 flux sources29, 57, different mechanisms of gas transport, different permeability of the soil, etc.
The statistical distribution of CO2 flux of each survey was modelled (see Methods) with the combination of a population characterised by a high mean CO2 flux value, HF population, and one characterised by a lower mean CO2 flux value, LF population (Fig. 3).
The mean CO2 flux value of LF and HF populations range from 14gm−2 d−1 to 135gm−2 d−1 and from 1,270gm−2 d−1 to 6,580gm−2 d−1 respectively (Table 1). The high mean CO2 flux values of HF populations clearly indicate that they are fed by the CO2 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 CO2 flux from a wide variety of ecosystems ranges from 0.2gm−2 d−1 to 21gm−2 d−1 (ref. 57 and reference therein). Roughly in agreement with these typical biogenic-derived fluxes, the coupled analyses of the CO2 flux and the isotopic composition of the CO2 efflux performed at Solfatara in March 2007 (see Supplementary Information Fig. S1) indicated a mean biogenic CO2 flux of 26 (±3) g m−2 d−1 (ref. 58). Since 2004, the mean CO2 flux of LF populations increased to 136gm−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 CO2 within the soil, resulting in soil CO2 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 CO2 flux of LF population (Table 1). Therefore, even though mean CO2 flux values of LF populations are consistently 1 to 2 orders of magnitude lower than mean CO2 flux values of HF populations, since 2004, the LF population clearly represents a mixture of biogenic and deeply-derived CO2. This suggests that deep-sourced CO2 degassing widely affects both Solfatara crater and its surroundings.
To characterise the spatial distribution of CO2 fluxes, all 13,158 CO2 flux measurements from 1998 to 2016 have been used to produce a map of CO2 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 CO2 flux as a cut-off. The threshold value selected was 50gm−2 d−1, the 95th percentile of the biogenic soil CO2 fluxes defined on the basis of the isotopic compositions of the CO2 efflux of the March 2007 survey58.
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 CO2 (see Methods). Solfatara DDS is defined in yellow to red in Fig. 4.
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 CO2 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 deposits60 that can act as a permeability barrier for the rising deep CO2. In fact, Astroni deposits are locally constituted by massive and fall out and ash surge deposits with thickness from ~10 to ~30m, affected only by mesoscale normal faults with lengths of tens of centimeters and displacements of a few centimeters60.
However, where normal faults deform and cut the Astroni deposits with up to metric dip separation61 (e.g., at Via Antiniana area and Pisciarelli), high CO2 fluxes occur at the surface (Fig. 4), suggesting that anomalous CO2 pressures can be present below this low permeability layer. The roughly NE-SW-elongated CO2 flux anomalies in the northern part of Pisciarelli further support the probable continuity of the CO2 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 CO2 fluxes, and of the total amount of CO2 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 CO2 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 CO2 flux is greater than 50gm−2 d−1 is higher than 50% (see Methods). The DDS extent ranges from 0.45×106m2 in 1998 to more than 1×106m2 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, CO2 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 CO2 flux passed from typical background values (10–40gm−2 d−1) in the pre‐2003 period to values higher up to 1 order of magnitude. The 2003 increase of CO2 fluxes and the DDS enlargement were already interpreted as due to the “first arrival” of hydrothermal CO2 at the surface, in the peripheral areas of Solfatara62.
The total amount of CO2 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 CO2 sources active in the area, it provides a good estimate of the hydrothermal CO2 release because CO2 fluxes from the hydrothermal source are much higher than those from biogenic sources. For example, assuming a constant biogenic background CO2 flux of 26gm−2 d−1 over the entire surveyed area, the biogenic CO2 output would result in ~38 t d−1. Considering this biogenic CO2 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.
In this section, variations in the CO2 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” CO2 emission (LF population) and the extent of DDS (Fig. 6a and b), and in another, the total CO2 output (Fig. 7a).
The mean CO2 flux of the LF population, which will be referred here as SBF (Solfatara background flux), increases from 20–30gm−2 d−1 in 1998 to 136gm−2 d−1 in the last campaign of June 2016, following a scattered, but almost continually increasing trend (Fig. 6a).
It is worth noting that typical mean CO2 fluxes generated by a biogenic source in the soil (20–30gm−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 CO2, indicating that areas previously unaffected by an anomalous CO2 degassing started to release deeply-derived gas. In fact, the SBF increase is correlated with an important change of the spatial pattern of CO2 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.5km2 to ~1km2 (Table 1 and Fig. 6b) and SBF increased from 20–30gm−2 d−1 to 40–50gm−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.2km2 (Table 1) and SBF increased to 70–100gm−2 d−1. We observe that the DDS extent reached 80–90% of the 1.4 km2 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 plume25. Repeated episodes of mud emission3, 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 CO2 fluxes increasing soil water content, changing soil permeability, dissolving soil CO2, 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 DTO measured during the first two campaigns (1998 and 2000) is relatively high (1,300–1,500 t d−1) while the two following campaigns, performed in 2003, are characterised by the lowest values (750–800 t d−1). Since then, a trend of increasing DTO begins, accompanying the ongoing unrest of CFc. The highest CO2 fluxes were measured in the last six campaigns (after September 2014) when DTO values, excluding November 2015, are above 1,500 t d−1 and reach peak values of 2,800 t d−1 and 1,900 t d−1 in January and March 2015, respectively.
The increase of DTO is mainly due to the increase of the degassing rate in areas of high CO2 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 CO2 flux increased from ~500gm−2 d−1 and ~400gm−2 d−1 during 2003–2013 to ~800gm−2 d−1 and ~1,100gm−2 d−1 after 2013, respectively (Supplementary Information Fig. S3). The most important increase occurred in Pisciarelli area where the CO2 flux passed progressively from ~100gm−2 d−1 to ~300 with a peak value of ~1,000gm−2 d−1 in 2015.
The relative increase of CO2 fluxes is particularly relevant at Pisciarelli, where the CO2 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 activity25, 53.
The DTO represents only a fraction of the total CO2 output, as it does not include the contribution from fumarolic vents. The vent contribution to the Solfatara CO2 budget most probably increased in the last years concurrently with the abovementioned increase in the fumarolic activity. Measurements of the CO2 flux performed in the 1980’s (refs 63 and 64) suggest, in fact, that the CO2 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 CO2 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 CO2 flux from the vents of Pisciarelli area, that was evaluated at only ~18 t d−1 in March 2009, reaches values of 170 t d−1 to 490 t d−1 in 2012–2015 (ref. 31). Finally, a high total (diffuse+vents) CO2 flux of 2,300 to 4,600 t/d from the Solfatara crater was estimated in March 2016 based on a new portable DIAL laser system, while the same method applied at Pisciarelli gave 266±211 t d−1 (ref. 33).
In Fig. 7b the CO2 release is compared with the seismicity4 and ground deformation54. Since 2005, the occurrence of earthquakes (located at depth <3km, ref. 66) increased together with the beginning of an uplift phase. Since 2005, a general correlation emerges between the increase of CO2 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 CO2 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 CO2 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 2000s3, 4. Here we show, for the first time, the modelled CO2 output returned by the last published simulation of the process3. The modelled CO2 output from the entire simulation domain (Fig. 2) shows roughly the same increasing trend of the measurements (Fig. 7a), even if the observed CO2 fluxes are systematically lower than the modelled CO2 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 CO2. 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 CO2 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 modelling3.
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 system67, a process that would be still ongoing. This depressurization is clearly in contrast with the evolution of total CO2 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 CO2 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 CO2 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 CO2 flux released by Solfatara is compared with the mean volcanic plume CO2 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 CO2 flux to the atmosphere similar to a “medium-large” volcanic plume. If, instead, we consider the highest total CO2 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 CO2 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 CO2 from hydrothermal sites is potentially globally relevant considering also that can reach very high values up to the 10–60kt d−1 estimated for Yellowstone17, 70.
The long-term volcanic CO2 degassing at Solfatara, and degassing from hydrothermal systems in general, can contribute to refining estimates of the volcanic CO2 contribution to the atmosphere, and can aid the ability to assess the possible role of CO2 degassing from hydrothermal systems.
CO2 flux measurements were performed using the accumulation chamber method, which is based on the measurement of the CO2 concentrations over time, inside an inverted chamber placed on the ground27–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,000ppm of CO2. 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.8L volume, to the IR and vice versa by a pump (~0.0167Ls−1). The AC is internally equipped with ring-shaped perforated manifold that re-injects the gas favouring the mixing in the chamber. The CO2 concentration in the circulating gas is acquired every 0.025s 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 CO2 flux is computed in real time from the rate of CO2 concentration increase in the chamber (dCCO2/dt) according to the relation CO2 flux=cf×dCCO2/dt (ref. 28). The proportionality factor (cf) between dCCO2/dt and the CO2 flux was determined before each survey by laboratory tests, during which imposed CO2 fluxes typically in the range 10 to 10,000gm−2 d−1, were measured over a “synthetic soil” made of dry sand (10cm-thick) placed inside a plastic box. In each calibration test, the cf value was computed from the linear best-fit line of CO2 flux vs. dCCO2/dt.
The measured soil CO2 flux and the measurement locations of each survey are reported in Supplementary Information Dataset S1, together with the soil temperature at ~10cm depth measured at the same time and at the same site of CO2 flux measurement. The soil temperature data is provide for completeness but is not discussed in this work.
The CO2 flux data was analysed by statistical methods to define the statistical parameters of the flux, which offer insight into the origin of degassed CO2.
In volcanic-hydrothermal areas, the CO2 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 CO2 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 points71. 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 CO2 fluxes28, 29, 57. Since the so computed M i value refers to the logarithm of the CO2 flux values, the mean value of CO2 flux was estimated using a Montecarlo simulation procedure (Table 1).
The estimated mean CO2 flux values have been used in literature to compute the CO2 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 CO2 output obtained from the GSA can be affected by arbitrary choices in the partitioning procedure29. 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.
In order to obtain a reliable estimation of the total CO2 output and to produce maps of the CO2 flux, the CO2 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 CO2 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 CO2 flux in this study), here performed using the sgsim algorithm of the GSLIB software library73. Since the sGs assumes multi-Gaussian distribution of the attribute to be simulated, the CO2 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 CO2 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×10m (except for the map reported in Fig. 2 for which a 2×2m 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 CO2 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 CO2 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 CO2. According to ref. 29 in this work, the DDS is considered as the area where the probability, among the 100 realizations, that the simulated CO2 flux is higher than the biogenic CO2 flux threshold over 50%.
The simulated flux values by sGs were used also to compute the total CO2 release by diffuse degassing. The total CO2 release is computed for each realization by summing the products of the simulated CO2 flux value at each grid cell by the cell surface. The mean of the values of total CO2 release computed for the 100 realizations are assumed as the characteristic diffuse total CO2 output (DTO, see above) for each period. The standard deviation of the 100 estimates is assumed as the DTO uncertainty.
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), that in turn causes the condensation of the steam25, would cause the heating of the system. The occurrence of condensation-heating was successively returned by a physical modelling approach3. The ref. 25 proposed a method to compute the fraction of the condensed steam starting from the fumarolic compositions. Because the computation in the original work is not detailed, here we illustrate the mass balance on which the computation is based on. The mass balance of a batch of gas which passes from reservoir to fumaroles, can be expressed in terms of number of moles of CO2 and H2O in the reservoir (n CO2,eq, n H2O,eq), in the fumarole (n CO2,fu, n H2O,fu) and the number of moles of water removed (n H2O,rem) and/or added (n H2O,add) during the fluid transfer. Assuming that secondary processes do not affect CO2, the mass balance equations are:
We can express the number of moles of water removed and/or added as the product of a proportionality factor (f rem, f add) by the original number of moles (n H2O,eq):
Considering that n CO2,fu=n CO2,eq, equation (6) can be written as:
From equation (7), the variable f is given by:
Finally, equation (8) is solved considering and :
The variable XCO2,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 H2O-H2-CO2-CO gas system. In detail, we used the following expressions38, 74:
where the partial pressures of water and CO2 are assumed equal to their fugacities.
This study has benefited from funding provided by INGV (project COHESO), DECADE (an initiative from the Deep Carbon Observatory) and by the Italian Presidenza del Consiglio dei Ministri Dipartimento della Protezione Civile (DPC), INGV-DPC Research Agreement 2012–2015, Progetto V2 “Precursori di eruzioni”. This paper does not necessarily represent DPC official opinion and policies. We thank Dr. Rebecca Astbury for comments and language editing that improved the quality of the manuscript. We wish to thank J. Lowenstern, A. Di Muro and an anonymous reviewer for the helpful suggestions which improved the quality and the clarity of the manuscript.
C.C., G.C. designed and leaded the fieldwork through the time, and with F.F. drafted the manuscript with input from all of the authors. C.C., A.R., B.E., C.S., L.M. and R.A. participated in the fieldwork, and reviewed the manuscript.
The authors declare that they have no competing interests.
Electronic supplementary material
Supplementary information accompanies this paper at doi:10.1038/s41598-017-06941-2
Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.