3D Numerical Simulations of Non-Volcanic CO₂ Degassing in Active Fault Zones Based on Geophysical Surveys

Abstract

An integrated approach that combines geophysical surveys and numerical simulations is proposed to study the processes that govern the fluid flow along active fault zones. It is based on the reconstruction of the architecture of the investigated fault system, as well as the identification of possible paths for fluid migration, according to the distribution of geophysical parameters retrieved by multi-methodological geophysical prospecting. The aim is to establish, thanks to constraints deriving from different types of data (e.g., geological, geochemical and/or hydrogeological data), an accurate 3D petrophysical model of the survey area to be used for simulating, by numerical modelling, the physical processes likely taking place in the imaged system and its temporal evolution. The effectiveness of the proposed approach is tested in an active fault zone in the Matese Mts (southern Italy), where recent, accurate geochemical measurements have registered very high anomalous values of non-volcanic natural emissions of CO₂. In particular, a multi-methodological geophysical survey, consisting of electrical resistivity tomography, self-potential and passive seismic measurements, integrated with geological data, was chosen to define the 3D petrophysical model of the investigated system and to identify possible source geometries. Three different scenarios were assumed corresponding to three different CO₂ source models. The one that hypothesizes a source located along the fault plane at the depth of the carbonate basement was found to be the best candidate to represent the test site. Indeed, the performed numerical simulations provide CO₂ flow estimates comparable with the values observed in the investigated area. These findings are promising for gas hazards, as they suggest that numerical simulations of different CO₂ degassing scenarios could forecast possible critical variations in the amount of CO₂ emitted near the fault.

Article Highlights

• Non-volcanic CO₂ degassing processes that govern the fluid flow along active fault zones are studied

• An approach based on high-resolution 3D ERT surveys and numerical simulations is proposed to model the fault architecture and the fluid uprising

• The methodological procedure applied to a test site allowed us to reproduce the gas saturation and CO₂ flow values measured in situ

1 Introduction

Faults and related cataclasites and damage zones can control the fluid migration in the upper crust (e.g., Caine et al. 1996; Agosta and Kirschner, 2003; Faulkner et al. 2010). Depending on the nature of the fault rocks, the texture and thickness of the cataclastic zone, degree of cementation and amount of displacement, they can act as conduits or barriers for the fluid migration (Caine et al. 1996; Bense and Person 2006; Annunziatellis et al. 2008; Smeraglia et al. 2016). Recent studies demonstrate a close correlation between gas releases and extensional tectonic regimes and confirm that both play a key role in creating pathways for the rising of gases from the micro- to the macro-scale (e.g., Kämpf et al. 2013; Brune et al. 2017; Tamburello et al. 2018; Italiano et al. 2019).

In the last 20 years, a growing interest is noticed in studying the process of carbon dioxide degassing in active fault zones with the aim to define the relationship between gas flux on the ground surface and faults (Miller et al. 2004; Cappa et al. 2009; Tamburello et al. 2018; Chiodini et al. 2020). In particular, an interesting challenge is to quantify the deeply derived CO₂ released into the atmosphere, which could be significant and hence represent an additional feature to be considered in terms of human health, ecosystems, groundwater and soils. Therefore, the characterization and modelling of fault zones, where large amounts of CO₂ are released, are fundamental for understanding the processes that govern the flow of fluids along the fault zones, which convey the gases towards the surface. This is especially true for active faults where fluid migrations (gases and liquids) can trigger earthquakes (Irwin and Barnes 1980; Noir et al. 1997; Miller et al. 2004; Cappa et al. 2009; Waldhauser et al. 2012; Fischer et al. 2017; Chiarabba et al. 2018; Girault et al. 2018).

In this framework, we propose an integrated approach based on geophysical surveys and numerical modelling to reconstruct the tectonically active system architecture and to simulate the rising fluid dynamics, respectively. In the last decade, numerous geophysical methods have been successfully applied to the study of CO₂ degassing phenomena, in both volcanic and non-volcanic areas (e.g., Pettinelli et al. 2010; Revil et al. 2011; Byrdina et al. 2014), while numerical modelling has proven to be a useful tool for analysing the temporal evolution of complex geological systems (e.g., Secchi and Schrefler 2012; Peiffer et al. 2015). In particular, our approach proposes a combined use of electrical resistivity tomography, self-potential and horizontal-to-vertical spectral ratio surveys to build a 3D accurate petrophysical model of the investigated system with the help of geological information available for the study area. Then, multi-dimensional discrete dynamical models are suggested to simulate the temporal evolution of CO₂ flows and the degree of gas saturation along the reconstructed fault system. The effectiveness of the proposed procedure is tested in an area located in the Matese Mts (southern Apennines, Italy), where, associated with the segmented Colle Sponeta fault, massive non-volcanic CO₂ degassing occurs (Ascione et al. 2018). Specifically, an innovative high-resolution 3D electrical resistivity tomography (ERT) prospecting, integrated with a self-potential measurement profile and two passive seismic soundings, has been employed to obtain a detailed imaging of the Colle Sponeta fault zone and to identify possible paths for CO₂ flow migration. The retrieved 3D geophysical model has been combined with geological and geochemical data available for the survey area to characterize the rock properties and the stratigraphic sequence of the investigated volume. Then, by using the TOUGH2 (Transport Of Unsaturated Groundwater and Heat) numerical simulator and its graphical interface Petrasim, we built a 3D petrophysical model of the test area and simulated natural CO₂ degassing for different choices of gas sources, with the purpose to reproduce the observed CO₂ flux. In the following, before presenting the results of these analyses, two sections are dedicated to the state of the art and the illustration of the proposed approach.

2 State of Art

The natural degassing of CO₂ from the geosphere occurs along major tectonic structures, hydrothermal and volcanic systems. In active fault zones, this transport takes place through permeable cataclastic rocks and damage zones formed by fractures that commonly create complex networks (Caine et al. 1996; Sibson 1996, 2000; Guerriero et al. 2013). The processes that govern the fluid flow along fault zones and convey gases towards the Earth’s surface depend on the permeability of the system and are related to rock–water–gas interactions. The fault zone structure can change with the faulting activity, and the fractured rock and damage zone permeability can therefore be highly variable over time. For example, fault cores made of porous breccias and cataclasites can behave as conduits for gas and/or fluid escape or as barriers when they are formed of very fine rocks or the open pore space is filled with mineral precipitation (Caine et al. 1996). However, the type of host rock, the degree of structural heterogeneity and the spatial distribution of fractures within fault zones are primary control factors on the migration pathways of non-volcanic carbon dioxide and, more in general, of fluid flows (e.g., Guerriero et al. 2013). Therefore, the detailed reconstruction of the architecture of an active fault zone is fundamental for identifying possible in-depth circuits for non-volcanic carbon dioxide migration, whose knowledge is useful to simulate the time evolution of the system for gas hazard mitigation.

In the last 30 years, many laboratory and field studies have been performed to identify and characterize in particular seismically active faults (e.g., Caine et al. 1996; Seront et al. 1998; Cello et al. 2000; Wibberley et al. 2008; Yönlü et al. 2017). Paleo-seismological analyses, for example, have been used for interpreting the seismogenic behaviour of active faults and defining magnitude, movement rates, return times and elapsed time (Galadini and Galli 1999; Michetti et al. 2005; Camelbeeck et al. 2007; Galli and Naso 2009; Rockwell et al. 2010; Galli and Peronace 2014; Zygouri et al. 2015). Moreover, geological studies along fault zones allowed the identification of optimal paleo-seismological sites, i.e. areas that potentially subtend not only the surficial fault exposition but also a good continuity and stratigraphic quality (Papanikolaou et al. 2015 and references therein). However, in specific conditions, these elements may fail for erosive phenomena or, more simply, due to a lack of competition between movement rates and exogenous dynamics (Papanikolaou et al. 2015). In these cases, the use of indirect methodologies, such as geophysical techniques, plays a relevant role as they allow identifying possible expositions and defining geometry and size of the potential deposits associated with the accommodation space resulting from the rock movements, such as the basins related to the activity of Quaternary normal faults in the Apennines (e.g, Gori et al. 2017). Over the last decades, indeed, the geophysical methods have been successfully applied to the study of fault zones allowing accurate volumetric reconstructions of buried fracture networks and detection of degassing processes in volcanic and non-volcanic areas (e.g., Chow et al. 2001; Caputo et al. 2003; Reiss et al. 2003; Revil et al. 2008; Byrdina et al. 2009, 2014; Pettinelli et al. 2010; McCalpin et al. 2011). In particular, the electrical and electromagnetic investigations appear to be among the most appropriate methods for revealing spatial distributions of carbon dioxide in active fault zones due to the possibility to describe the investigated subsurface in terms of electrical parameters, such as electrical resistivity and natural electrical charge polarization, which strongly depend on the physical properties (e.g., porosity, fluid content, grain size) of the geological formations. For example, electrical resistivity imaging identifies zones of influence of gas vents as both conductive and resistive areas depending on the geological setting and the physical, chemical and biological soil environment of the investigated vents, while self-potential data usually show a strong correlation between CO₂ degassing and negative self-potential anomalies in areas with focused release (e.g., Rogie et al. 2000; Arts et al. 2009; Byrdina et al. 2009; Pettinelli et al. 2010; Sauer et al. 2014).

Although geophysical surveys, integrated with geological, geochemical and hydrogeological data, allow the definition of consistent petrophysical models of active fault zones, as well as the identification of anomalous sectors correlated to possible fluid migration pathways, the physical processes likely driving the underground fluid flow circulation can be confirmed/inferred only by numerical simulation of flow and transport processes in permeable systems. Over the past two decades, numerical modelling of underground fluid flows has proven to be a very powerful tool for studying a wide range of geological processes, such as the generation and migration of hydrocarbons and gases along active fault zones (e.g., Caine and Forster 1999; Caillet and Batiot 2003; Yang 2006; Gessner et al. 2009; Monreal et al. 2009; Chi et al. 2010; Wang et al. 2015; Changgui et al. 2019). This is thanks to its ability to define and estimate the primary control parameters of the flow of fluids and their interactions, to simulate the velocity trend of underground fluids and to describe the space–time variations of pressure, hydraulic load, fluid and gas saturation in porous, fractured and/or faulted rocks. The multiphase and multicomponent fluids transport in such systems can be described by a set of partial differential equations that express the principles of mass and energy conservation. Except in special cases, these equations cannot be solved analytically but must be solved numerically by discretization in space and time to obtain an equivalent system of linear algebraic equations, which can therefore be solved with direct or iterative approaches. Anyway, it is worth to note that the modelling reliability is strongly dependent on the initial conceptual model developed for the investigated system, or rather on its schematization in terms of thermophysical and chemical properties. This means that the reliability of the simulations is linked to quantity and quality of surface and sub-subsurface data available for the study area coming, for example, from geological, geochemical and hydrogeological investigations.

3 A New Approach for Studying CO₂ Degassing Based on Geophysical Surveys and Numerical Simulations

The proposed approach aims to reproduce non-volcanic CO₂ flows observed along active fault zones through the integration of geophysical surveys and numerical simulations. This study is part of a broader research project whose final goal is the development of a useful tool for toxic gas hazard assessment that can be generalized for other kinds of fluids or gases according to their specific properties.

The methodological path is described schematically in Fig. 1. The first step of the procedure consists of building a 3D realistic petrophysical model of the investigated system and of identifying possible geometries of the system gas source. To this aim, a multi-methodological geophysical investigation consisting of 3D high-resolution electrical resistivity tomography (ERT), self-potential (SP) and horizontal-to-vertical spectral ratio (HVSR) techniques is proposed. Specifically, the ERT and SP measurements are addressed to define the architecture of the fault damage zones and to localize possible preferential CO₂ migration paths in terms of distributions of the electrical parameters. The choice of such methodologies lies on the recognized strong dependence of the electrical properties of water-bearing rock systems on many factors relevant to the storage of CO₂ (porosity, fracturing, water saturation, quantity and type of ions present in the water, pH, cation exchange capacity of clay minerals and temperature). But the HVSR measurements are suggested as further support to the interpretation of the ERT survey, specifically to characterize the thickness of the superficial deposits and to estimate the bedrock depths (e.g., Sauret et al. 2015). It is worth noting that the use of other geophysical methodologies could be equally effective for the purposes of the proposed approach. For example, electromagnetic (e.g., magnetotellurics, time-domain electromagnetics) and gravimetric methods would be successful in the case of large volumes and/or exploration depths. Indeed, such methods allow the characterization of the buried volume to be investigated in terms of resistivity value distribution and density contrasts, respectively, which are suggestive of fluid and/or gas occurrence within the surveyed faulted system.

Fig. 1

Methodological path to study the dynamics of non-volcanic CO₂ flux and gas saturation observed in active fault zones

The second step of the proposed procedure is focused on the reconstruction of the 3D petrophysical model of the examined complex geological system and on the identification of key parameters for studying variations in the system response due to CO₂ degassing. This is achieved by integrating the geophysical model with field data from geological, geochemical and hydrogeological studies available in the literature or specifically planned for the surveyed area.

The third phase is based on the numerical modelling, which aims at simulating the dynamics of CO₂ along the faulted system and at estimating the flow values within the investigated underground volume. To this end, the EOS2 module of TOUGH2 numerical simulator (see Appendix A for details) is proposed because it solves the thermodynamic equation for a mixture of water and CO₂ and accounts for the non-ideal behaviour of gaseous CO₂ and its dissolution in the aqueous phase. EOS2 (see Appendix B) is an updated version of the module originally developed by O’Sullivan et al. (1985) for the description of fluids in gas-rich geothermal reservoirs, which may contain masses of CO₂ with fractions ranging from a few per cent to occasionally 80% or more (Atkinson et al. 1980). In general, the workflow used for the numerical modelling is based on the following steps: discretization of the model geometry; definition of input parameters and boundary conditions; simulation of CO₂ flux distribution and gas saturation for the gas sources assumed for the investigated faulted system; and interpretation and validation of the simulated results by comparison with the experimental data observed in the survey area.

4 Application to the Test Site of Ciorlano

The proposed approach has been applied to the Ciorlano site in the Matese Mts (southern Apennines, Italy), which recent accurate geological and geochemical analyses classify as one of the areas with the highest non-volcanic natural emissions of CO₂ ever measured on Earth (Ascione et al. 2018).

4.1 Geological Background

The Ciorlano area is characterized by the Meso-Cenozoic sedimentary succession of the Matese Mts (Fig. 2). It consists of Triassic–Cretaceous shallow-water carbonates, covered by Miocene shallow-water limestones, slope marls, and finally, siliciclastic sediments (Vitale and Ciarcia 2018). The area is defined by numerous high-angle faults with NW–SE, E–W and N–S directions, related to the extensional tectonic activity during the late Quaternary (Galli et al. 2017; Ascione et al. 2018). These faults control the continental depocentres and karst basins distribution and deform the exposed Quaternary deposits. The survey area lies on an alluvial terrace at about 280 m above sea level (a.s.l.) composed of a succession of silts and clays interbedded with pyroclastic layers dissected by N–S-trending normal faults (Ascione et al. 2018). The area is characterized by an advective and localized gas leakage with dominant CO₂ associated with gas vents recognized on the field by scarce and/or total absence of vegetation (Ascione et al. 2018). The alignments of the gas vents are consistent with the mapped fault segments active during the late Quaternary. In particular, the anomalous values of concentration and flow of CO₂ at the Ciorlano site, which also correspond to high He concentrations in soil gas, clearly follow the Colle Sponeta fault escarpment (Fig. 3). Therefore, the gas vents, as well as the anomalous diffuse gas emissions, are strongly controlled by the fault network, which is able to maintain a high permeability along the fault direction.

Fig. 2

Geological map of the N–W sector of Matese Mts (southern Apennines) showing the non-volcanic gas emission sites and the location of the study area (red rectangle) (from Vitale and Ciarcia 2018, modified)

Fig. 3

Spatial distribution of CO₂ soil concentration and flux at the Ciorlano survey area (Matese Mts, southern Italy). The anomalous concentrations and flow values, which also correspond to high He concentrations in soil gas, are aligned along the N–S fault escarpment of Colle Sponeta (from Ascione et al. 2018, modified). The four parallel white lines indicate the profiles, 20 m apart, used for the 3D ERT investigation; the yellow triangles (T1 and T2) show the location of the two HVSR soundings

4.2 Geophysical Model

The black rectangle in Fig. 3 shows the geophysical survey area, where the four parallel white lines and the two yellow triangles indicate, respectively, the profiles along which the electrical resistivity tomography and the passive seismic (HVSR) surveys were performed. The branch ERT 3 in Fig. 3 also represents the profile investigated by the SP measurements.

4.2.1 3D Electrical Resistivity Tomography Survey

As is well known, the ERT technique aims to obtain detailed images of the underground buried structures by inversion of a great number of apparent resistivity measurements acquired on the Earth’s surface according to 2D, 2.5D or 3D electrode disposals (e.g., Drahor et al. 2007; Chambers et al. 2012; Di Maio and Piegari 2012; Di Maio et al. 2016a). Conventional designs for 3D ERT surveys consist of multi-electrode cables distributed along parallel profiles (2.5D mode) or arranged according to 2D regular meshes (3D mode). The latter usually have small sizes being limited by the number of electrodes that can be managed by the modern georesistivimeters (generally 48, 72 or 96); therefore, they allow high-resolution investigations only for the shallowest portions of the subsurface. To obtain both high-resolution resistivity distribution and greater investigation depth, in this work, an innovative unconventional 3D acquisition scheme was used for the electrical characterization of the survey area. Specifically, the data acquisition was carried out along the four parallel spreads in Fig. 3 and implemented in three steps as shown in Fig. 4a. Quadrupoles sequence generation for such geometry included both “common cable” quadrupoles (i.e. transmitting and receiving electrodes belonging to the same cable) and “cross cable” quadrupoles (i.e. all the possible combinations given by transmitting and receiving dipoles not belonging to the same cable) (Santarato et al. 2011). Such a procedure is thus able to provide a real volumetric distribution of apparent resistivity values and not an interpolated one, like that obtained by the pseudo-3D (i.e. 2.5D) ERT technique. The apparent resistivity data were collected by the IRIS Syscal Pro Switch multi-channel georesistivimeter using the pole-dipole as measurement array. According to the target of the survey, 48 electrodes spaced of 10 m were used, which allowed us to investigate an area of 470 m × 60 m. As shown in Fig. 4a, each 3D acquisition overlaps the adjacent spreading to ensure optimal coverage of the whole survey area. To attain the optimal 3D acquisition configuration for the study area, the ERT-Lab Sequencer software (distributed by Geostudi Astier S.r.l.) was used, which allows the setting of complex acquisition geometries. In particular, a sequence consisting of 1116 current injections and 10,722 quadrupoles was provided for each of the three acquisitions, which required a data collection time of about 50 min. As example, Fig. 4b shows the distribution of the measurement points in the volume corresponding to the 3D acquisition no. 3 in Fig. 4a, where the position of the pole–dipole array corresponding to the first measurement of the “common cable” sequence along the ERT 4 profile (Fig. 3) is indicated.

Fig. 4

a Data acquisition scheme used for the 3D ERT survey in the Ciorlano area. b Distribution of the measurement points in the volume corresponding to the 3D acquisition no. 3 in (a); A2, M25 and N48 indicate the position of the pole-dipole array corresponding to the first measurement of the “common cable” sequence along the ERT 4 profile in Fig. 3

The collected apparent resistivity values (about 32,000) were implemented in a single 3D dataset, which considers the elevation of the measurement electrodes. The very high number of well-distributed measurements allowed the definition of the high-resolution apparent resistivity pseudo-volume (470 × 60 × 105 m³) shown in Fig. 5 after highly noisy data filtering.

Fig. 5

Apparent resistivity pseudo-volume from the whole 3D dataset consisting of the three acquisition steps shown in Fig. 4a (see text for details)

The 3D inversion of the whole apparent resistivity dataset was performed by the ERTLab64 inversion software (produced by Geostudi Astier Srl, Livorno, Italy, and Multi-Phase Technologies LLC, Nevada, USA), which uses a finite elements approach to model the subsoil by adopting tetrahedral meshes to correctly include surface topography. The inversion process, based on a least squares smoothness constrained method (Constable et al. 1987; LaBrecque et al. 1996), needs an initial guess which is generally given by the homogeneous half-space derived from the statistical analysis of the acquired apparent resistivity values. Throughout the inversion iterations, the effect of non‐Gaussian noise is appropriately managed using a robust data weighting approach (Morelli and LaBrecque 1996), which controls the noise through the standard deviation reweighting algorithm on the outliers and the remodulation of the roughness factor to ensure a final model robust enough against possible artefacts (Fischanger et al. 2019). The latter occurring in the presence of strong resistivity contrasts between the features of interest and the hosting medium (Constable et al. 1987; Francese et al. 2009) as in our survey area, where two main lithologies characterized by a sharp contrast in their electrical properties (i.e. very conductive clays overlapped on high resistive carbonate basement) are present. The final 3D resistivity model of the investigated subsurface was obtained by regular tetrahedral mesh with a resolution of 2.5 m starting from an initial homogeneous model of 100 Ωm. A fast convergence of the model response to the field data was reached with a RMS (root mean square) error lower than 5%. Finally, by using the ViewLab3D software (Geostudi Astier S.r.l.), a three-dimensional visualization of the investigated volume in terms of electrical resistivity (ρ) distribution was retrieved as shown in Fig. 6. As it can be seen, the whole volume is characterized by resistivity values ranging from 1.50 Ωm to 10,000 Ωm (Fig. 6a). Based on the lithology of the test area, the range 1.50 Ωm ≤ ρ ≤ 20 Ωm can be associated with a clay layer (blue volume in Fig. 6b) that, in the easternmost part of the survey area, reaches a thickness of about 60 m below the ground level (b.g.l.). Within this layer, more resistive patterns are observed (30 Ωm ≤ ρ ≤ 70 Ωm, green volumes in Fig. 6b) potentially attributable to a mixture of water and CO₂, which would cause an increase in the resistivity values (Revil et al. 1999; Byrdina et al. 2009). Such a mixture likely also involves the top of the carbonate basement (Fig. 6c, d) supposed at a depth of about 40–50 m b.g.l. in the western part of the survey area and of about 60–70 m b.g.l. in its central and eastern sectors (dashed black line in Fig. 6d). The basement, in fact, is recognized as the deep electro-layer characterized by a wide range of resistivity values (200 Ωm ≤ ρ ≤ 10,000 Ωm) likely corresponding to a different fracturing degree of the limestone bedrock, which appears to decrease with increasing depth, and/or to the amount of fluids filling the fractures. Such assumptions would justify the relatively low resistivity values observed in the uppermost part of the basement, and the very high resistivity anomalies (ρ > 2000 Ωm) visible in the western and central portions of the investigated volume (dark red zones in Fig. 6), which are probably ascribable to the presence of CO₂ gas accumulations inside the fractured basement. In particular, for the central resistive anomaly (deep black circle in Fig. 6d), this hypothesis could be supported by the small degassing pipe, with groundwater bubbling phenomenon due to CO₂ leakage, clearly visible on the field at the time of the survey in correspondence with ERT lines 3 and 4 (shallow black circle in Fig. 6d). This phenomenon is likely due to the contact between aquifer and CO₂, the latter coming from fractures of the carbonate basement.

Fig. 6

a 3D resistivity model of the survey area (1.50 Ωm ≤ ρ ≤ 10,000 Ωm). The images b–d highlight volumes characterized by specific ranges of resistivity values: b 1.50 Ωm ≤ ρ ≤ 20 Ωmn (blue shade) c 30 Ωm ≤ ρ ≤ 70 Ωmn (green shade) and d 1000 Ωm ≤ ρ ≤ 10,000 Ωm (red shade). T1 and T2 in d show the location of the two HVSR surveys. The dashed black line in (a) sketches the hypothesized hidden fault

The high resolution of the 3D ERT survey in the Ciorlano area also highlighted the possible in-depth geometry of the fault structure. Indeed, by cutting the resistivity volume in Fig. 6a according to slices at different altitudes a.s.l. along the Z-axis (Fig. 7), the resistive-conductive contact, which is observed near the progressive 235 m along the X-axis, could be indicative of the fault plane, which appears clearly identified in the depth range from 34 to 84 m from the ground surface. In this regard, it is worth noticing that the fault is well recognized as a central plane of a resistive zone at about 70 m b.g.l. (altitude 246 m, Fig. 7a). Instead, at a depth of 20 m from the surface level (altitude 296 m, Fig. 7d) the fault is detected as a centre planar volume of an inhomogeneous conductive zone, where the less conductive sectors could suggest the presence of fractures filled by CO₂.

Fig. 7

ERT slices at depth a Z = 246 m a.s.l. b Z = 266 m a.s.l., c Z = 276 m a.s.l. and d Z = 296 m a.s.l.

To support the interpretative hypotheses suggested by the ERT model, ambient seismic noise and self-potential measurements along profile were performed.

4.2.2 Horizontal-to-Vertical Spectral Ratio Soundings

The horizontal-to-vertical spectral ratio (HVSR) technique (e.g., Nakamura 1989; Bonnefoy-Claudet et al. 2006) is based on passive measurements of ambient seismic noise consisting of small-amplitude continuous oscillations due to both natural sources (large-scale weather perturbations, wind, ocean waves, etc., frequency < 1 Hz) and anthropic sources (vehicular traffic or industrial activities, frequency > 1 Hz). The HVSR method is applied to vertical and horizontal records of ambient vibrations measured by a three-component broadband or short-period seismometer. The ratio of the averaged horizontal-to-vertical frequency spectrum is used to determine the fundamental site resonance frequency, fr, which, in the simplest case (soil layer above a half-space), depends on layer thickness (H) and shear wave velocity (Vs), according to the following equation (Nakamura 1989)

$$f_{r} = \frac{{V_{s} }}{4H}$$

(1)

Therefore, the method highlights the occurrence of resonance phenomena and provides an estimation of the frequencies at which the ground motion can be amplified, as a result of site effects induced by the presence of stratigraphic discontinuities, and thus by variations in the geophysical parameters that characterize the subsurface. To obtain the Vs value, constraints on the depth of a stratigraphical contact marked by a seismic impedance contrast are needed, or, vice versa, to assess the depth of the interface between two layers the knowledge of Vs for the upper layer is required (e.g., Lachet and Bard 1994; Delgado et al. 2000; D’Amico et al. 2004; Haefner et al. 2010; Paudyal et al. 2013; Castellaro 2016). Compared to the simple case indicated by Eq. (1), more complex equations are taken into account to transform a HVSR curve into a Vs profile, according to specific inversion software.

In the Ciorlano survey area, two seismic soundings (T1 and T2 in Fig. 3) were carried out using a Tromino® digital tromograph (by MoHo s.r.l.) equipped with three short-period (proper frequency equal to 4.5 Hz) electrodynamic velocimeters oriented N–S, E–W and Up–Down, respectively. At each site, the data were collected in a time interval of 22 min with a sampling frequency of 128 Hz and then processed by the Grilla software (Micromed s.p.a.), which applies guidelines for processing ambient seismic noise records according to the SESAME recommendations (SESAME 2004). Specifically, the HVSR data processing involved a preliminary analysis by subdividing the signal acquired for each component into non-overlapping time windows of 20 s, which were analysed to separate the transient portion of the signal from its most stationary parts and to obtain the individual spectral components of the signal by the fast Fourier transform. The spectra were smoothed by a triangular smoothing function with a width equal to 10% of the central frequency. Afterward, for each frequency, the averages of the single-component spectra (i.e. N–S, E–W and Up–Down) were computed for all the analysed time windows, and the HVSRs of each window were then averaged to obtain the final HVSR mean curve. It is worth noting that, in some cases, the HVSR curves show anomalous peaks due to local cultural noise sources (e.g. industrial activities, water pumps, vehicular traffic), which could be characterized by frequencies comparable to those under study and therefore require to be removed. Consequently, the time stability of the spectral ratios and their directionality (i.e. projection of the HVSR along different directions, from 0° to 180°) were analysed with angular step of 10°. Finally, the international criteria for stability of results (i.e. reliability of the HVSR curve and clarity of the HVSR peak) (SESAME 2004) were verified before to interpret the HVSR curves in Fig. 8 coming from the above-described data processing applied to the two seismic soundings carried out in the Ciorlano survey area (T1 and T2 in Fig. 3). In particular, Fig. 8a,e exhibits the mean HVSR curves along with the 95% confidence interval, where the occurrence of more than one frequency peak suggests the possible presence of different seismic surfaces overlapping the rock basement, which would support the results of the 3D ERT survey. To evaluate the “quality” of the two single noise measurements, and consequently their reliability, the following criteria were applied:

Fig. 8

Results of the HVSR (H/V) data analysis obtained for the T1 and T2 soundings in Fig. 3.(a, e) Mean H/V curves (red lines) and standard deviations (black lines); (b, f) single-component spectra (i.e. N–S, E–W and Up–Down); (c, g) stationary values; (d, h) signal directions

1. recording duration long enough to produce “robust” estimates of the average field of ambient vibrations;

2. temporal stationarity of the spectral ratio;

3. isotropy of the signal in terms of spectral ratios;

4. absence of electromagnetic noise;

5. shape of the H/V curve in the analysed frequency range.

Following the above criteria, we note that the T1 site does not show any directivity (Fig. 8d) of the peaks present in the HVSR curve. There are no indications of electromagnetic noise (absence of disturbances), and the overall trend of the HVSR curve in the analysed frequency range is stationary for at least 50% of the recording duration. This means that there is a substantial stationarity (Fig. 8c) within the analysed single time window. Furthermore, the maxima of the HVSR ratio are characterized by a localized decrease in the amplitude of the vertical component spectrum (purple curve in Fig. 8b) with respect to the two horizontal components. On the other hand, for the T2 site we observe a strong directivity (Fig. 8h), particularly in the 0.4–2 Hz frequency band, centred on N80°. The overall trend of the curve does not remain stationary for at least 30% of the measurement duration (Fig. 8g). Consequently, there is no stationarity (Fig. 8g) within the analysed single time window, especially for frequency lower than 1 Hz, while there are no indications of electromagnetic noise. Finally, in correspondence with the maximums of the HVSR ratio, the decrease in the amplitude of the vertical component spectrum with respect to the two horizontal ones is not as evident as for the T1 site. Therefore, the seismic measurement carried out at site T1 shows high reliability, while the one performed at site T2 shows medium reliability.

To assess the Vs profiles, the code provided by the Grilla software (Micromed s.p.a) was applied to the HVSR curves in Fig. 8a,e. The code, which is based on the simulation of the surface wavefield (Rayleigh and Love) in multi-layer systems characterized by plane-parallel discontinuities (e.g., Castellaro 2016), provides a synthetic HVSR curve, to be compared with the experimental one, and the corresponding Vs profile. The modelling procedure is based on a priori constrain (generally depth or Vs of the first layer) coming from direct tests, such as penetrometer or well logs, available for the study area. Indeed, as well know, the lack of any constraint would generate infinite models (i.e. infinite Vs–H combinations) that could satisfy the observed HVSR curve. In particular, to interpret the two experimental HVSR curves, a stratigraphical log from a borehole located in a neighbouring area of the Ciorlano test site was used as initial guess for the HVSR data fit that, according to a trial-and-error approach, provided the synthetic HVSR curves and the relative Vs profiles shown in Fig. 9. As given in Table 1, which synthesizes the results of the trial-and-error process, four-layer 1D models provided the best fit with respect to the experimental HVSR curves. The blue points of the synthetic HVSR curves in Fig. 9a,c correspond to the three seismic contact surfaces reconstructed through the above forward modelling: The shallower the contact surface, the higher the frequency peak. We note that the synthetic HVSR curve related to the T1 site well fits the experimental curve (Fig. 9a), while for the T2 site a good match between experimental and synthetic curve (Fig. 9c) is observed for the highest frequency peak (30–40 Hz). The retrieved models are more appropriate than single-layer models consisting of a low-velocity layer over a half-space corresponding to the basement rock, because for a complex system, such as the Ciorlano survey area, continuous vertical variation of velocity and stratifications with appreciable velocity contrasts must be considered. Indeed, despite the seismic sounding at the T2 site did not show high reliability for the lower-frequency peak, the retrieved 1D models provided Vs and depth values that are consistent with the results from the 3D ERT survey (Fig. 6d). Specifically, the HVSR analysis confirms the presence of a bedrock-like layer at a depth of about 70 m from the ground level in the eastern part of the investigated area where the T1 seismic sounding is located, which rises in its central–western part where the T2 sounding is sited, with values of Vs that correlate well with a fractured carbonate basement.

Fig. 9

a, c Comparison between the experimental H/V curve (red lines) and the synthetic H/V curve (blue lines) provided by the trial-and-error modelling for the T1 and T2 sites in Fig. 3. b, d Vs profiles for the T1 and T2 sites

Table 1 Stratigraphical features and Vs values of the subsurface model used for modelling the two HVSR soundings (T1 and T2 in Fig. 3). The lithology of the layers derives from a stratigraphical well-log located close to the T1 and T2 sites

4.2.3 Self-Potential Profile

The SP method aims to retrieve the distribution of buried anomalous electric charge concentrations by measurements of potential differences occurring naturally on the ground surface (e.g., Parasnis 1997; Sharma 2002; Revil and Jardani 2013; and the references therein). The main SP sources are ascribable to redox reactions around ore deposits or metallic bodies, electrochemical processes taking place between regions of different ionic concentrations, electrokinetic phenomena generated by underground fluid and heat fluxes (e.g., Sen 1991; Birch 1998; Naudet et al. 2004; Revil and Leroy 2004; Castermant et al. 2008; Heinze et al. 2019). In particular, the streaming potential, due to the fluid circulation along pores and/or permeable fractures systems, proved a key parameter in the geothermal and hydrogeological research (e.g., Ishido et al. 1983; Fournier 1989; Revil et al. 1999, 2006; Fagerlund and Heinson 2003), in seismic and volcanic activity forecasting (e.g., Mitzutani et al. 1976; Di Maio and Patella 1991; Di Maio et al. 2000; Zlotnicki and Nishida 2003; Grobbe and Barde-Cabusson 2019), and, recently, in detection and monitoring CO₂ degassing processes along active fault zones due to the direct link observed between the SP anomalous pattern and the carbon dioxide release in porous and fractured rock systems (see Sect. 2.2).

The SP measurements in the Ciorlano study area were performed along the profile ERT 3 in Fig. 3. The data were acquired through the gradient technique using two unpolarizable electrodes and a voltmeter with high input impedance. An electrode spacing equal to that used for the ERT profiles (i.e. 10 m) was chosen, which allowed the acquisition of the measurement points shown in Fig. 10a with black dots. It is interesting to note that a strong negative SP anomaly is found precisely in correspondence with the Colle Sponeta fault, which, as shown in Fig. 7, is centred at 235 m along the SP profile. This occurrence confirms the interpretative hypothesis coming from the ERT survey, which attributes to the shallower conductive portion of the central part of the investigated area (Fig. 7) a permeable fracture zone that acts as a preferential path for CO₂ degassing.

Fig. 10

a Observed SP data (dots) and curve (red line) calculated by assuming a 2D inclined sheet b as anomaly source model. The calculated curve has been obtained for the following model parameters: x₀ = 237.41 m, z₀ = 96.12 m, a = 10.11 m, θ = 121.52°

To obtain quantitative information on the observed SP anomaly source, a global optimization method based on a Genetic-Price hybrid Algorithm (GPA) (Di Maio et al. 2016b) has been applied to the collected SP data. The inversion algorithm combines features of controlled random search and genetic algorithms and requires minimization of the following misfit function:

$$\sigma^{2} = \frac{1}{{\left( {N - m} \right)}}\mathop \sum \limits_{i = 1}^{N} \left( {V_{i,obs} - V_{i,c} } \right)^{2}$$

(2)

where N is the number of acquired data, m is the number of parameters and V_i,_obs and V_i,c are, respectively, the measured SP data and the computed SP data obtained from the assumed (forward) model. In particular, a 2D inclined sheet source has been assumed to model the hypothesized dipping planar fault damage zone. The general expression of the SP anomaly, $V_{i,c = sheet} ,$ generated by a such source (Fig. 10b) at any point of coordinate x_i on a profile perpendicular to the strike (e.g., Sundararajan et al. 1998; Gobashy et al. 2020), is given by the following equation:

$$V_{i,sheet} = kln\left\{ {\frac{{\left[ {\left( {x_{i} - x_{0} } \right) - acos\theta } \right]^{2} + \left( {z_{0} - asin\theta } \right)^{2} }}{{\left[ {\left( {x_{i} - x_{0} } \right) + acos\theta } \right]^{2} + \left( {z_{0} + asin\theta } \right)^{2} }}} \right\}$$

(3)

where k is the electric dipole moment, x₀ and z₀ are the horizontal position and the depth of the sheet centre, respectively, a is the half-width of the sheet and θ is the sheet inclination with respect to the horizontal axis (Fig. 10b). The synthetic SP curve (red line in Fig. 10a) that best matches the main feature of the experimental SP data is obtained for inversion parameters k = 643.47 ± 131.95 mV, x₀ = 237.41 ± 2.62 m, z₀ = 96.12 ± 1.93 m, a = 10.11 ± 3.94 m and θ = 121.52° ± 0.02°, which characterize the dipping planar fault damage zone and correlate well with the results of the ERT survey. We note that the above inversion values represent the optimal solutions for the source parameters estimated as average values over a series of 100 runs of the GPA optimization procedure, each of which involved a maximum of 5·10⁵ iterations and a convergence limit of 10⁻⁶. Significantly, the value of the sheet angle returns a fault dip angle of π - θ ≅ 59° (see dashed black line in Fig. 6a), which perfectly falls in the range of observed and predicted dips for normal faults (e.g., Jacskon and White, 1989). However, as expected, the inversion curve does not properly fit the experimental curve in its western and eastern parts, where the presence of other buried SP anomaly sources is reasonable to assume. Indeed, looking at Fig. 6d, the structural high detected by the ERT measurements in the western part of the survey area and attributed to the top of the carbonate basement, as well as the presence of N–S-trending flowing streams east of the Colle Sponeta fault that cross the survey area, could originate electrokinetic and junction potentials (e.g., Maineult et al. 2006; Schutze et al. 2015) responsible of the SP trend observed at both ends of the experimental curve for which more detailed modelling is required that is beyond the scope of the present work.

4.3 Petrophysical Model

Based on the results of the geophysical investigations carried out in the study area (see Sect. 4.2), a synthetic geological model 470 m long, 60 m width and 120 m thick was considered (Di Maio et al. 2019). Specifically, from the inversion of the 3D ERT data (Fig. 6a), the model was approximated by a sequence of three layers (Fig. 11) representative of the lithostratigraphic units of the study area, i.e.

• 1 shallow and impermeable clay layer with a thickness of 10 m.

• 2 water-saturated layer of alluvial deposits (silt and clays) with a thickness ranging from 50 m to 60 m b.g.l. proceeding from west to east.

• 3 carbonate basement from a depth of 58 m to 70 m b.g.l. proceeding from west to east.

Fig. 11

3D geological model used to simulate CO₂ emissions along the Colle Sponeta fault. The latter is described by the darker cells in the central part of the model (see text for details)

The chosen volume was then discretized into a grid consisted of 4230 regular-shaped elementary cells and referenced according to a XYZ coordinate system with the X-axis in the WE direction, the Y-axis in the NS direction and the Z-axis-oriented downward (Fig. 11). The discretization along the three axes is as follows: 47 cells 10 m long in the X-direction, 6 cells 10 m long in the Y-direction and 15 layers in the Z-direction, consisting of cells with a thickness of 2 m, 12 m and 10 m for the first, second and third layer, respectively. As shown in Fig. 11, the discretized geological model includes the Colle Sponeta fault damage zone (brown cells) as a homogeneous and isotropic region between two parallel sub-vertical planes equidistant from the progressive 235 m along the X-axis. The fault dip angle was fixed at 70° based on the inclination of the central resistive–conductive contact in Fig. 6a (dashed black line) attributed to the damage zone, which was also confirmed by the results of the forward modelling of the SP data (Fig. 10). The fault is assumed to extend along the Y- and Z-axes by 60 m and 120 m, respectively (Fig. 11), with normal kinematics and a displacement of about 10 m, as reported in Ascione et al. (2018).

The petrophysical parameters and the initial and boundary conditions required for the numerical simulation were retrieved from geological and geochemical literature data (Celico 1990; Colleselli and Colombo 1996; Eppelbaum et al. 1996; Mielke et al. 2017; Ascione et al. 2018) and are summarized in Tables 2 and 3. As shown in Table 2, the second layer is assumed to be the most permeable between the shallow clay layer and the underlying fractured carbonate basement. A decrease in porosity values with depth is also supposed. However, the simulations were performed by varying the values of the petrophysical parameters in different ranges to verify the stability of the presented solutions. Moreover, to model the fluid flow processes, a two-component (H₂O + CO₂) biphasic flow was assumed, while to define the initial thermodynamic state of the analysed system pressure, gas saturation degree and CO₂ partial pressure were chosen as primary variables (Table 3). Specifically, to each cell of the grid in Fig. 11 has been assigned a pressure corresponding to the lithostatic load, a gas saturation value of 70% and a CO₂ partial pressure given by 70% of the lithostatic load (Ascione et al. 2018). Finally, as concerns the boundary conditions, very large volumes (V = 10⁵⁰ m³) have been assumed for the outermost and surface blocks of the grid, which is equivalent to imposing that the thermodynamic conditions must not change due to mass flows.

Table 2 Parameters used for characterizing the discretized model in Fig. 11

Table 3 Parameters used for describing the initial conditions of the discretized model in Fig. 11

4.4 Numerical Simulations

The 3D numerical modelling has been realized by using the PetraSim software, which is the graphical interface for the TOUGH2 simulator family by an interactive 3D environment that includes mesh generation, parameter definition and result visualization. To define a source cell, injection rate (or flux) and enthalpy must be fixed as input parameters. The details about the model setup, the actual equations that are solved and the solution scheme are given in Appendices A and B. In the following, the results of numerical simulations performed by assuming three different geometries of a deep source supplying endogenous CO₂ are presented and discussed. Numerical simulation of the temporal evolution of CO₂ degassing in such three source models was aimed at defining the system that better than the others is able to reproduce the values of CO₂ flow and gas saturation measured in the survey area (Ascione et al. 2018). Specifically, for the characterization of the source cells a constant CO₂ flow rate of 5.78 kg/s (Ascione et al. 2018) and an enthalpy value of 505 J/kg were assigned to all the source cells of the grid, while a slightly higher rate (6.0 kg/s) and an enthalpy value of 205 J/kg were assumed for the water–steam mixture. Finally, the gas saturation of the cells was assumed equal to 98%, according to Ascione et al. (2018).

In Fig. 12a, the first investigated model is shown, where the CO₂ source cells (in yellow) are located in the hanging wall of the fault and, specifically, at the top of the carbonate basement. This model assumes that the permeable damage zone (fractured rocks) hosted within the carbonates is sealed by the overlying clays and silts and that the preferential migration path of water and CO₂ comes from the east, specifically, from the Matese Mts (Ascione et al. 2018). Based on this hypothesis, numerical simulations of water and CO₂ upward migration dynamics were performed to estimate CO₂ flow and gas saturation values within the investigated underground volume. The time required to reach stationary conditions was around 60 years. Figure 12b, c shows the results after 10 years and 60 years of simulated dynamics in terms of distribution of CO₂ flux vectors and gas saturation iso-surfaces, respectively. It is interesting to note the formation of convective circulation of the CO₂ flow in the hanging wall of the fault in the package composed of silts and clays. Indeed, because of buoyancy forces, CO₂ is expected to rise at the top of the permeable layer. However, a fraction of CO₂ dissolves in the aqueous phase, it increases its density and the associated negative buoyancy force can induce convective circulation that will carry dissolved CO₂ downward, while causing additional dissolution of CO₂ into upflowing waters that are poor in CO₂ (Pruess et al. 1999).

Fig. 12

I source system: a 3D geological model of Fig. 11 with location of the source cells (in yellow). CO₂ flux distribution vectors and gas saturation iso-surfaces (GS) after 10 years b and 60 years c of simulated dynamics

The second hypothesis proposed for the source system is illustrated in Fig. 13a, where the CO₂ source cells (in yellow) are positioned along the fault plane up to the top of the carbonate basement; therefore, it is assumed that the fault core acts as a preferential conduit for the CO₂ upward migration. Also for this model, the numerical simulation was performed until to reach the stationary conditions, which were obtained after about 100 years of simulation. As examples, Figs. 13b,c shows the results after 10 years and 100 years, respectively, of simulated dynamics. In this case, it is worth noting that after 10 years (Fig. 13b), CO₂ storage bags form in the silt–clay package, 60 m thick, at the easternmost part of the investigated volume, i.e. in the hanging wall block of the fault. The observed convective cells, indeed, seem to reproduce the sectors with higher resistivity identified within the shallow conductive layer of the 3D electrical resistivity model (green volumes in Fig. 6b) and attribute to CO₂ migration pathways.

Fig. 13

II source system: a 3D geological model of Fig. 11 with location of the source cells (in yellow). CO₂ flux distribution vectors and gas saturation iso-surfaces (GS) after 10 years b and 100 years c of simulated dynamics

Figure 14a illustrates the last proposed source model. It was formulated on the basis of the high-resolution 3D electrical resistivity tomography which allowed an extremely detailed characterization of the investigated subsurface volume, highlighting the possible fault plane geometry. Indeed, from the slices of Fig. 7, the damage zone could be identified as a large resistive zone at a depth of 84 m b.g.l. (Fig. 7a). Therefore, for this third system, the source cells (in yellow) extend for about 70 m along the X-axis up to the top of the carbonate basement (Fig. 14a), and it is assumed that the entire degraded zone at the fault boundaries acts as a preferential path for CO₂ migration. Once again, the model was used to simulate the upward migration dynamics of fluids and CO₂ and to estimate CO₂ flow and gas saturation values within the studied subsurface volume. The system evolution was simulated in a time interval of about 100 years. Figure 14b,c shows the simulation results after 10 years and 100 years, respectively. As for the first two systems, CO₂ flow convective circulation is well evident at the hanging wall block in the silt and clay package characterized by greater permeability.

Fig. 14

III source system: a 3D geological model of Fig. 11 with location of the source cells (in yellow). CO₂ flux distribution vectors and gas saturation iso-surfaces (GS) after 10 years b and 100 years c of simulated dynamics

To evaluate the differences in the CO₂ flow and gas saturation estimates retrieved from the numerical simulations for the three different CO₂ source geometries, the results obtained for the same grid cell placed on the Earth’s surface near the fault (identified with the cell number 3690) have been analysed. In particular, for this cell gas saturation and CO₂ flow values have been represented as a function of time in Fig. 15a,b, respectively, thus permitting the comparison between the estimates provided by the numerical simulations and the values measured in situ, which for the Ciorlano area are of the order of 33 g∙d⁻¹ m⁻² and 76% (v/v), respectively (Ascione et al. 2018). From the comparison shown in Fig. 15a, it turns out that the first two source systems return gas saturation values comparable with those measured in situ after about 60 years and 90 years, respectively, whereas the third hypothesized source system provides, in the simulated time interval, much lower gas saturation values, not comparable with those measured in situ. A reasonable explanation for this behaviour can be found in the size of the highly permeable area; in fact, the greater the area fractured and permeated by fluids, the longer the time for the silt and clay package to become saturated. By comparing the results of the three source models in terms of CO₂ flow estimates with the field values (Fig. 15b), it is possible to observe that the first and third hypothesized source systems return much larger and much lower flow estimates, respectively. On the contrary, the second source system appears to be able to reproduce flow values comparable with those observed in situ. Interestingly, a CO₂ flow of 33 g∙d⁻¹ m⁻² is found after 10 years of simulation. For longer simulation times, the estimated CO₂ flow does not appear to undergo significant variations, and after 100 years, it reaches an asymptotic value, which is indicative of stationary conditions.

Fig. 15

a Gas saturation and b CO₂ flow values estimated for the cell n. 3690 at different time steps for the three hypothesized source systems. The red dashed line marks the value measured in the survey area

In summary, the simulation results at different time steps for the three analysed source models suggest that a relatively narrow fractured area along the Colle Sponeta fault, i.e. the second model (Fig. 13), could well represent the source of CO₂ degassing in the survey area. Interestingly, for this model convective cells of CO₂ flow form in the eastern part of the investigated volume that reduce to a single large convective cell from 10 to 100 years, likely due to an increase in the saturation degree. On the western side, a large amount of gas under pressure is envisaged, which well correlates with the high resistivity anomalies observed in the shallow sedimentary succession in the western side of the investigated area.

5 Conclusions

A methodological approach based on high-resolution 3D ERT surveys and numerical modelling, which integrates geological and geochemical data is proposed to study the diffuse CO₂ degassing along fault zones. The suggested procedure has been applied to the Ciorlano area (Matese Mts., southern Apennines, Italy) at the aim to localize and model preferential non-volcanic CO₂ migration pathways near to a gas vent detected in the study area.

An innovative high-resolution 3D ERT technique based on non-conventional electrode configurations, generally applied to the civil engineering field, has been employed and optimized to investigate the chosen complex geological environment. The obtained well-resolved subsurface model of the test area, validated by self-potential and passive seismic surveys, provided a detailed geostructural and physical characterization of the investigated subsurface volumes, which allowed to identify and model, in terms of resistivity contrasts, the Colle Sponeta fault zone and the complex migration path of the non-volcanic CO₂ flow. In particular, the 3D resistivity model obtained for the Ciorlano area suggests that the relatively high resistivity (about 60–70 Ωm) patterns observed in the shallow clay layer at the easternmost part of the surveyed area, at a depth of about 60 m b.g.l., could be associated with possible CO₂ migration paths, which would cause an increase in the resistivity values. Furthermore, the very high resistivities (1000–2000 Ωm) observed in the central portion of the studied volume up to the maximum exploration depth (about 100 m b.g.l.) are very likely attributable to CO₂ volumes in the underlying fractured carbonate basement. Finally, the 3D geophysical model, suitably integrated with recent and accurate geological and geochemical analyses, drove the numerical simulation of the CO₂ uprising along the Colle Sponeta fault zone. Three different scenarios have been hypothesized corresponding to three different CO₂ source models. In particular, a model that hypothesizes the CO₂ source system located along the fault plane at the depth of the carbonate basement was found to be the best candidate to represent the surveyed area. Indeed, the performed numerical simulations allow to well reproduce the gas saturation and CO₂ flow values measured in situ.

It is worth noting that the proposed methodology can be helpful in many geological research fields, such as geothermal exploration and gas hazard assessment, due to the possibility of defining with high accuracy geometry and physical characteristics of fault systems in terms of resistivity contrasts, whose knowledge is fundamental to identify possible preferential paths for gas movements, as well as gas storage reservoirs. Furthermore, future perspectives of the present research provide for the use of new wireless systems for electrical tomography, which would allow to investigate very large areas and to reach high exploration depths in relatively short times, thus allowing 4D ERT surveys for monitoring carbon dioxide degassing. Such monitoring could provide useful constraints on the numerical modelling for predicting future hazard scenarios.

Appendices

Appendix A: The TOUGH2 numerical simulator

TOUGH2 (Transport Of Unsaturated Groundwater and Heat) numerical simulator solves the mass and energy balance equations that describe flows of heat and multiphase, multicomponent fluid mixtures in porous and fractured media. It is particularly effective for describing water properties and phase transition considering the relative permeability of the geological formations crossed by fluids and the effects of capillarity pressure. The code, developed by Lawrence Livermore Laboratory Berkeley (Pruess et al. 1999), uses a modular approach through interchangeable modules, called EOS (equation of state), which define components, phases and thermodynamic properties of the studied system, such as viscosity, enthalpy and fluids density. Using the EOS modules, TOUGH2 has both the ability to simulate different thermodynamic conditions and the flexibility to be applied to different areas of interest. In Fig. 

Fig. 16

Pruess et al. 1999)

Scheme of the structure and architecture of the TOUGH2 numerical simulator (modified from

16, a scheme of the TOUGH2 structure is shown. As it can be seen, the code is structured on two main one-dimensional arrays which contain, respectively, the primary thermodynamic variables and the secondary thermophysical parameters required to assemble the equation governing the flow (i.e. density, viscosity, enthalpy, etc.). The set of the primary variables depends on the choice of the EOS module, e.g. pressure, temperature and CO₂ partial pressure for the EOS2 module used in the present work.

The basic mass- and energy-balance equations solved by TOUGH2 can be written in the general form:

$$\underbrace {{\frac{d}{dt}\mathop \smallint \limits_{{V_{n} }} M^{k} dV_{n} }}_{accumulation} = \underbrace {{\mathop \smallint \limits_{{S_{n} }} \left( {\overrightarrow {{F^{k} }}\cdot\vec{n}} \right) dS_{n} }}_{advection} + \underbrace {{\mathop \smallint \limits_{{V_{n} }} q^{k} dV_{n} }}_{{\begin{array}{*{20}c} {production} \\ {/injection} \\ \end{array} }} ,$$

(4)

where Vn is an arbitrary sub-domain of the analysed flow system, which is bounded by the closed surface Sn, with $\vec{n}$ the normal vector to the surface element dS_n, M is the mass (or energy) per volume, with k = 1, …, NK indicating the mass components, and k = NK + 1 the heat component, $\vec{F}$ is the heat [kg/(m²s)] or mass [J/(m²s)] flow and q denotes the production (q < 0) or injection (q > 0) of fluids and heat.

In the mass balance of a system characterized by more than one component occurring in several phases β (e.g. liquid, gas, etc.), the mass accumulation term M^k in Eq. (4) takes the form:

$$M^{k} = \varphi \mathop \sum \limits_{\beta } S_{\beta } \rho_{\beta } X_{\beta }^{k} ,$$

(5)

where the porosity, φ, is multiplied for the sum of each phase of the k-component, S_β and ρ_β are, respectively, the saturation and the density of the phase β and X_β^k is the mass fraction of the k-component present in the β phase. Similarly, the heat accumulation, $M^{KN + 1} ,$ is given by two contributions:

$$M^{KN + 1} = \left( {1 - \varphi } \right)\rho_{R} C_{R} T + \varphi \mathop \sum \limits_{\beta } S_{\beta } \rho_{\beta } \mu_{\beta } .$$

(6)

The first contribution takes into account the matrix heat provision, where ρ_R and C_R are, respectively, the density and the specific heat of the rock and T is the temperature. The second contribution stands for the heat of each phase, where μ_β is the specific internal energy of the phase β.

The advective mass flux $\overrightarrow {{F^{k} }}$ in Eq. (4) is equal to the sum over all the individual phase fluxes, $\overrightarrow {{F_{\beta } }} ,$ weighted by the mass fraction, X_β^k:

$$\overrightarrow {{F^{k} }} = \mathop \sum \limits_{\beta } X_{\beta }^{k} \overrightarrow {{F_{\beta } }} .$$

(7)

The flux $\vec{F}_{\beta }$ is computed using the Darcy’s law:

$$\overrightarrow {{F_{\beta } }} = \rho_{\beta } \overrightarrow {{v_{\beta } }} = - K\frac{{K_{r\beta } \rho_{\beta } }}{{\mu_{\beta } }}\left( {\nabla P_{\beta } - \rho_{\beta } \vec{g}} \right),$$

(8)

where $\overrightarrow {{v_{\beta } }}$ is the Darcy’s velocity in the phase β, K is the absolute permeability, K_rβ is the relative permeability to phase β, μ_β is the viscosity coefficient, P_β is the fluid pressure related to the phase β and $\overrightarrow {{g{ }}}$ is the gravity acceleration vector. In this case, the sink and source contributions represent a mass rate per volume. As concerns the heat flux $\overrightarrow {{F^{NK + 1} }}$, it includes conductive (Fourier’s law) and convective components:

$$\overrightarrow {{F^{NK + 1} }} = - \lambda \nabla T + \mathop \sum \limits_{\beta } h_{\beta } \overrightarrow {{F_{\beta } }} ,$$

(9)

where λ is the thermal conductivity and h_β is the specific enthalpy in the phase β.

Generally, the continuum Eqs. (4) are solved numerically by a space/time discretization procedure to obtain an equivalent set of linear algebraic equations resolved through direct or iterative approaches. In particular, TOUGH2 uses the integral finite difference (IFD) method (e.g., Edwards 1972; Narasimhan and Witherspoon 1976), which allows to create regular and non-regular three-dimensional grids. The discretization approach used in the IFD method and the definition of the geometric parameters are illustrated in Fig. 

Fig. 17

Spatial discretization and geometric data in the IFD method (from Pruess et al. 1999)

17.

Introducing appropriate volume averages, the accumulation term in Eq. (4) results:

$$\mathop \smallint \limits_{{V_{n} }} M dV = V_{n} M_{n} ,$$

(10)

where M is a volume-normalized extensive quantity and M_n is the average value of M over the domain V_n (Corapcioglu 1996). But the surface integral in Eq. (4) is approximated as a discrete sum of the average values of the (inward) normal component of the flux $\overrightarrow {{F^{k} }}$ over the surface segments A_nm in Fig. 14:

$$\mathop \smallint \limits_{{S_{n} }} \left( {\overrightarrow {{F^{k} }}\cdot \vec{n}} \right) dS = \mathop \sum \limits_{m} A_{nm} F_{nm} ,$$

(11)

where A_nm is the contact surface between the volume elements (or grid blocks) V_n and V_m and F_nm is the average value of the $\overrightarrow {{F^{k} }}$ component perpendicular to the A_nm surface, as shown in Fig. 17.

Finally, the production/injection term in Eq. (4), similarly to the accumulation term, becomes:

$$\mathop \smallint \limits_{{V_{n} }} q^{k} dV_{n} = V_{n} q_{n,}^{k}$$

(12)

where $q_{n}^{k}$ is an average mass rate.

Substituting eqs. (10), (11) and (12) in Eq. (4), a set of first-order ordinary differential equations with respect to time is obtained:

$$\frac{{dM_{n}^{k} }}{dt} = \frac{1}{{V_{n} }}\mathop \sum \limits_{m} A_{nm} F_{nm}^{k} + q_{n}^{k} .$$

(13)

The time is discretized as first-order finite differences, and the advection and production/injection terms (second member in Eq. (13)) are evaluated at the new time level t ^j+1 = t ^j + ∆t by a “fully implicit” procedure (e.g., Peaceman 1977) to obtain numerical stability needed for an efficient calculation of the multiphase flow (Pruess et al. 1999).

Thus, the time–space discretization generates the system of nonlinear algebraic equations of the type:

$$R_{n}^{k, j + 1} = M_{n}^{k, j + 1} - M_{n}^{k, j} - \frac{\Delta t}{{V_{n} }}\left[ {\mathop \sum \limits_{m} A_{nm} F_{nm}^{k, j + 1} + V_{n} q_{n}^{k, j + 1} } \right] = 0,$$

(14)

where $R_{n}^{k, j + 1}$ is the residual corresponding to the kth equation (k = 1, 2, …, NEQ; NEQ = NK + 1 with number of fluid components) related to the nth volume element (or grid block), V_n, at the t ^j+1 time step.

Since for each volume V_n there are NEQ equations, Eq. (14) represents, for a flow system characterized by NEL grid blocks, a set of NEL x NEQ coupled nonlinear equations, with NEL x NEQ unknown independent primary variables which define the state of the flow system at the time step t ^j+1. The set of eqs. (14) is solved by the Newton–Raphson iterative method, which stops at the iteration that provides residues satisfying the following predetermined relative convergence ($\varepsilon$) condition:

$$\left| {\frac{{R_{n,p + 1}^{k,j + 1} }}{{M_{n,p + 1}^{k,j + 1} }}} \right| \le \varepsilon ,$$

(15)

where p indicates the iteration index. The default (relative) convergence criterion used by TOUGH2 is $\varepsilon$ = 10^–5.

To solve eqs. (14), the knowledge of the thermophysical properties of fluid mixtures is required; such properties are provided by appropriate equation of state which are chosen according to the characteristics of the flow system to be studied.

Appendix B: EOS2 Module

As specified in Sect. 4.4 and Appendix A, the thermophysical properties of fluid mixtures needed for assembling the equations that govern mass and energy balance are provided by the EOS modules. The fluid-phase conditions are recognized from the numerical values of the primary variables. Each EOS module fulfils three additional important functions:

• the phase conditions pertaining to a given set of primary variables are identified for all the volume elements (grid blocks);

• the appearance (or disappearance) of phases is recognized as primary variables change during the Newton–Raphson iteration process;

• primary variables are switched and properly re-initialized in response to a change of phase.

The meaning of the primary variables/secondary parameters as defined in TOUGH2 essentially eliminates any direct link between the choice of the primary variables and the secondary parameters used to set up the flow equations. This provides maximum flexibility and convenience in choosing primary variables, as only secondary parameters are used in the flow equations.

For the study presented here, the EOS2 module was used, because it solves the thermodynamic Eq. (4) for a mixture of water and CO₂ (k = 1, …, NK = 2) and accounts for the non-ideal behaviour of gaseous CO₂ and its dissolution in the aqueous phase. Therefore, in Eq. (5) β_H2O = 2 (liquid and vapour phases) and β_CO2 = 1 (gas phase) were used. The primary variables for such two-phase conditions are gas-phase pressure, gas saturation and CO₂ partial pressure.

According to the Henry’s law, the partial pressure of an uncondensed gas, P_NCG, in the gas phase is proportional to the molar fraction of NCG, $x_{aq}^{NCG}$, dissolved in the aqueous phase (O’Sullivan et al. 1985):

$$P_{NCG} = K_{h} x_{aq}^{NCG} ,$$

(16)

where Kh is the Henry’s law constant which is the ratio of a compound’s partial pressure in air to the concentration of the compound in water at a given temperature.

Due to the inaccuracy of the Henry subroutine used by the EOS2 module at low temperatures (T < 30 °C), this subroutine has been replaced with a new procedure based on the correlation developed by Battistelli et al. (1997), as it is more appropriate for a mixture of cold water and CO₂ (T ≤ 25 °C) considered for our study. Tables 2 and 3 (see text) report the variables/parameters used for characterizing the discretized model in Fig. 8 and describing its initial conditions.