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

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 CO2. 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 CO2 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 CO2 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 CO2 degassing scenarios could forecast possible critical variations in the amount of CO2 emitted near the fault.


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 2 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 2 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 2 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 2 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 2 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 2 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 2 degassing for different choices of gas sources, with the purpose to reproduce the observed CO 2 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.

State of Art
The natural degassing of CO 2 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 1996Sibson , 2000Guerriero 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. 2009Byrdina et al. , 2014Pettinelli 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 2 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 3 A New Approach for Studying CO 2 Degassing Based on Geophysical Surveys and Numerical Simulations
The proposed approach aims to reproduce non-volcanic CO 2 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 2 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 2 (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. 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 2 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 2 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 2 and accounts for the non-ideal behaviour of gaseous CO 2 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 2 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 2 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.

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 2 ever measured on Earth (Ascione et al. 2018).

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 2 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 2 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.

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

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.
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 3 ) shown in Fig. 5 after highly noisy data filtering.
The 3D inversion of the whole apparent resistivity dataset was performed by the ERT-Lab64 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 1 3 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 2 , 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 2 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 2 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 2 , the latter coming from fractures of the carbonate basement.
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 resistiveconductive 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 2 .
To support the interpretative hypotheses suggested by the ERT model, ambient seismic noise and self-potential measurements along profile were performed.

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) (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 shortperiod (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  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: 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 planeparallel 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 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 1 3 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-anderror 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.

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. 1999Revil et al. , 2006Fagerlund 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 2 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 2 degassing.
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: 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: where k is the electric dipole moment, x 0 and z 0 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 0 = 237.41 ± 2.62 m, z 0 = 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 5 iterations and a convergence limit of 10 −6 . 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.

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 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.
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 2 O + CO 2 ) biphasic flow was assumed, while to define the initial thermodynamic state of the analysed system pressure, gas saturation degree and CO 2 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 2 partial pressure given by 70% of the lithostatic load (Ascione et al. 2018). Finally, as concerns the boundary conditions, very large volumes (V = 10 50 m 3 ) 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.

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 1 3 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 2 are presented and discussed. Numerical simulation of the temporal evolution of CO 2 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 2 flow and gas saturation measured in the survey area (Ascione et al. 2018). Specifically, for the characterization of the source cells a constant CO 2 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 2 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 2 comes from the east, specifically, from the Matese Mts (Ascione et al. 2018). Based on this hypothesis, numerical simulations of water and CO 2 upward migration dynamics were performed to estimate CO 2 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 2 flux vectors and gas saturation isosurfaces, respectively. It is interesting to note the formation of convective circulation of the CO 2 flow in the hanging wall of the fault in the package composed of silts and clays. Indeed, because of buoyancy forces, CO 2 is expected to rise at the top of the permeable layer. However, a fraction of CO 2 dissolves in the aqueous phase, it increases its density and the associated negative buoyancy force can induce convective circulation that will carry dissolved CO 2 downward, while causing additional dissolution of CO 2 into upflowing waters that are poor in CO 2 (Pruess et al. 1999).
The second hypothesis proposed for the source system is illustrated in Fig. 13a, where the CO 2 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 2 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 2 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 2 migration pathways. 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 2 migration. Once again, the model was used to simulate the upward migration dynamics of fluids and CO 2 and to estimate CO 2 flow and gas saturation values within the studied subsurface volume. The system evolution was 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 Layer number  Figure 14b,c shows the simulation results after 10 years and 100 years, respectively. As for the first two systems, CO 2 flow convective circulation is well evident at the hanging wall block in the silt and clay package characterized by greater permeability.
To evaluate the differences in the CO 2 flow and gas saturation estimates retrieved from the numerical simulations for the three different CO 2 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 2 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 −1 m −2 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 2 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 2 flow of 33 g•d −1 m −2 is found after 10 years of simulation. For longer simulation times, the estimated CO 2 flow does not appear to undergo significant variations, and after 100 years, it reaches an asymptotic value, which is indicative of stationary conditions.
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 2 degassing in the survey area. Interestingly, for this model convective cells of CO 2 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.

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 2 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 2 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 wellresolved 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 2 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 2 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 2 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 2 uprising along the Colle Sponeta fault zone. Three different scenarios have been hypothesized corresponding to three different CO 2 source models. In particular, a model that hypothesizes the CO 2 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 2 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.
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. 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 2 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: where Vn is an arbitrary sub-domain of the analysed flow system, which is bounded by the closed surface Sn, with ⃗ 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, ⃗ F is the heat [kg/(m 2 s)] or mass [J/(m 2 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: 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: 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 ��� ⃗ F k in Eq. (4) is equal to the sum over all the individual phase fluxes, ��� ⃗ F , weighted by the mass fraction, X β k : The flux ⃗ F is computed using the Darcy's law: where ��� ⃗ v 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 � ⃗ 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 ���������� ⃗ F NK+1 , it includes conductive (Fourier's law) and convective components: 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. 17.
Introducing appropriate volume averages, the accumulation term in Eq. (4) results: 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 ��� ⃗ F k over the surface segments A nm in Fig. 14: 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 ��� ⃗ 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: where q k n 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: 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:  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 2 (k = 1, …, NK = 2) and accounts for the non-ideal behaviour of gaseous CO 2 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 2 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 NCG aq , dissolved in the aqueous phase (O'Sullivan et al. 1985): (16) P NCG = K h x NCG aq , 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 2 (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. Acknowledgements The authors thank Dr. Sascha Brune and an anonymous reviewer for their insightful comments and suggestions that helped to improve the manuscript. Financial support to the present research from Italian National Agency for the Evaluation of Universities and Research Institutes (ANVUR), via Funding for Basic Activities Related to Research (Grant no. E61I18001660005), is gratefully acknowledged.
Funding Open access funding provided by Università degli Studi di Napoli Federico II within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons   Fig. 16 Scheme of the structure and architecture of the TOUGH2 numerical simulator (modified from Pruess et al. 1999) Fig. 17 Spatial discretization and geometric data in the IFD method (from Pruess et al. 1999) licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.