3D subsoil reconstruction of a mud volcano in central Sicily by means of geophysical surveys

The upwelling of fluids, subject to overpressure, along with discontinuities in the subsoil, causes the formation of geological structures known as mud volcanoes. These structures, very widespread in the world and in some cases located near inhabited centers, can represent a considerable risk factor for the population, as they can give rise to paroxysmal eruptions, even very violent. The assessment of the characteristics and structure of the subsoil of the areas affected by this phenomenon can prove to be a useful tool for risk mitigation. Non-invasive geophysical surveys were carried out in the area of the active cone of the Santa Barbara mud volcano in order to obtain a 3D characterization of the subsoil. Through the processing and integration of data derived from active and passive seismic surveys was obtained, a 3D model showing the seismostratigraphic subsoil structure. The electrical resistivity tomography surveys provided results comparable to those obtained from seismic surveys and supplied information on the perimetral areas of the mud volcano. The information obtained is useful to study the boundary conditions that influence short-scale activity. Furthermore, this study seeks to evaluate the possibility of using the proposed methodological approach for monitoring the variations that occur over time in the areas affected by mud volcanoes phenomenon.


Introduction
Mud volcanoes are geological structures that arise from the emission on the surface land or on the seabed of significant quantities of fluids, such as mud, water and hydrocarbons, in a variable percentage depending on their origin. The distribution of mud volcanoes on the earth's surface is controlled by geological and tectonic factors. In fact, these structures are mainly located in geological contexts characterized by active compression Le Pichon et al. 1990;Dimitrov 2002;Etiope et al. 2002;Martinelli and Judd 2004;Milkov 2005) or with high sedimentation rates (Graue 2000;Mascle et al. 2001).
Overpressure and under-consolidation are the key features that trigger the upward migration of water and gas-rich mud, often along fault planes and fracture lines (Brown 1990;Milkov 2000;Dimitrov 2002;Kopf 2002;Niemann and Boetius 2010). Mud volcanoes are positive morphologies, which have a central main emission duct and surface elements linked to the emission of fluids. Often, minor cones called griffins are present on the flanks of the volcanic building.
Although most of the mud volcanoes have the appearance of a conical hill, they can take on a great variety of morphologies. The morphology is mainly a function of duct size, pressure of interstitial fluids, frequency and type of emission, as well as viscosity of the emitted material. The density of the emitted material controls the maximum height that the mud volcano can reach.
These structures have aroused a great scientific interest especially for the strong impact of the paroxysmal and non-paroxysmal events that take place in many areas of the Earth (Mazzini et al. 2007;Davies et al. 2008;Madonia et al. 2011). The sudden release of large volumes of mud from structures located near populated areas can have harmful effects on the communities residing within their areas of influence, representing a source of natural hazard. An emblematic case is the powerful eruption of the Lusi mud volcano in Indonesia, starting in 2006, causing the displacement of some thousands of people and huge economic damage (Manga and Bonini 2012). Moreover, as a result of eruptive events deep cracks extending for several kilometres are often observed in surrounding areas (Baloglanov et al. 2018), causing significant damages to buildings, roads and pipelines (Gattuso et al. 2021).
Several scientific studies are focused on the application of geophysical surveys in mud volcanoes areas. In many of these researches, the survey methods used to image the near-surface structure of a mud volcano mainly consist of electrical resistivity surveys (Keller and Rapolla 1974;Ping-Yu et al. 2011;Zeyen et al. 2011), in some cases used together with active seismic surveys (Accaino et al. 2007;Rainone et al. 2015). Recently, Imposa et al. (2016) have applied passive seismic surveys in the mud volcano context. This study highlights how the passive seismic surveys can provide valuable information in defining the underground structure of a mud volcano. This work describes the results of some geophysical field surveys carried out in an area affected by the presence of a mud volcano, known as "Santa Barbara" mud volcano.
The mud volcano is located near the Santa Barbara village, hence its name, at an altitude of about 520 m above sea level (m asl), within the area of Caltanissetta municipality, in central Sicily (Fig. 1a-c). Although the area has always been known for the mud volcano phenomenon, over the years, it was turned into a residential area with buildings and public works. The activity in the area has been documented since 1800, by several scholars that have studied this area and described its paroxysmal events (Li Volsi 1826;Baldacci 1886;Roda 1968;Mezzadri 1989).
The last significant paroxysmal eruption occurred on 11 August 2008. During this event, the mud was thrown upward at about 30 m. The mud deposit created at the end of the event was about 5 m-thick with a diameter of about 200 m. In August 2009, another associated minor event occurred. Both of these events produced ground deformations in the surrounding area, modifications in the mud volcano morphology and changes in the number of the fluid emission points (Madonia et al. 2011).
The area of Santa Barbara mud volcano has never been investigated with applied geophysics surveys aimed at obtaining a 3D subsoil model. In this study, integrating the results of different geophysical techniques, a 3D reconstruction of the area affected by the presence of the emission points was obtained. This model shows the likely ascent pathway, and the shallower accumulation areas of the muddy material emitted at the surface. Moreover, an attempt to monitor the subsoil conditions over time was carried out in the area immediately surrounding the main emission point, comparing the 3D models obtained from two different field surveys, performed three months apart.

Geodynamic and geologic framework
The central Sicily area must be framed within the broader geodynamic context that involves the entire central Mediterranean region (Scandone 1979). It is the result of a complex deformation history linked to subduction and collision mechanisms, which have occurred starting from the Tertiary era, between the African and European margins (Malinverno and Ryan 1986;Dewey et al. 1989;Boccaletti et al. 1990;Lentini et al. 1994Lentini et al. , 1996De Guidi et al. 2014). The Caltanissetta basin, falls within the foreland fold and thrust belt domain, characterized by crustal structures with brittle deformation and on the surface by multiplicative structure with ductile deformation.
From a geological point of view, the area of Santa Barbara mud volcano is characterized by clay deposits (Tortorici et al. 2014) belonging to the Argille Brecciate formation (Early Pliocene) (Ogniben 1954) that widely outcrop in the area (Fig. 1d). This formation covers and/or intercalates to the Trubi formation, consisting of marls and marly limestones of Low Pliocene. The Trubi formation is in stratigraphic continuity but in geometric discordance above the Gessoso Solfifera formation; the latter consists of the evaporites deposited during the salinity crisis in the Mediterranean Sea, which occurred during the Upper Miocene (Messinian). The Middle-Upper Miocene clays (Argille Brecciate) (Ogniben 1954) are at the bottom of this geological succession.
The area has a hilly morphology, interrupted by mountain ridges, it has a considerable variety of shapes linked to the different lithologies outcropping. The hilly morphologies that mainly affect the clayey successions are characterized by modest mammellonary reliefs and by forms of ravine erosion. The flat portions are very limited, and they are linked to the alluvial belts of the main rivers.

Survey design
Overpressure and under-consolidation are the triggers that allow the upward migration of water and gas-rich mud, in mud volcanoes systems. The rising material has different physical and mechanical properties than the medium within which it moves, with particular reference to its density. Therefore, between these two media an impedance contrast should be present. With the aim to detect and localize these contrasts, consequently identifying the probable ascent paths 1 3 and accumulation areas of the emitted material, a geophysical field survey was carried out in the area where the emission points are present (Fig. 2a, b). Many scientific papers show the aptitude of the passive seismic single station technique to detect the presence of impedance contrasts in the subsoil (Yamanaka et al. 1994;Duval et al. 1995;Zhao et al. 2000;Panou et al. 2005;Imposa et al. 2016Imposa et al. , 2017aPappalardo et al. 2016Pappalardo et al. , 2018Pappalardo et al. , 2020Grassi et al. 2019Grassi et al. , 2021. In order to localize these contrasts, it is necessary to know the trend of shear waves velocity with depth, this data can be retrieved by performing active seismic geophysical surveys. Therefore, during the first field survey Horizontal to Vertical Spectral Ratio (HVSR) and Multichannel Analysis of Surface Waves (MASW) seismic prospections were carried out.
Furthermore, two 2D electrical resistivity tomography (ERT) surveys were carried out in order to identify in the subsoil the presence of the areas affected by emitted material through the association with the distribution of resistivity variations (Accaino et al. 2007;Ping-Yu et al. 2011;Zeyen et al. 2011;Rainone et al. 2015).
A second field survey was performed, after three months, carrying out passive seismic single-station samplings in the area where the main emission points are located, in order to attempt to evaluate the variations of impedance contrasts.

HVSR
The HVSR or H/V technique allows to identify the frequencies at which the ground motion is amplified by stratigraphic resonance computing the spectral ratio between the average of the horizontal components and the vertical one (Nogoshi and Igarashi 1970;Nakamura 1989). This method can be used to process seismic ambient noise (Lachet and Bard 1994;Tokimatsu 1997;Yamamoto 2000;Okada 2003;Molnar et al. 2018;Xu and Wang 2021), acquired in three components seismic single station surveys in order to obtain information on the seismic structure of the subsoil. The presence of an impedance contrast between different lithotypes or between rock portions with different mechanical features produces a peak in the H/V spectrum. If a velocity inversion is present, the H/V amplitude value drops below one (Castellaro and Mulargia 2009a). This phenomenon could be linked to the presence of water and/or materials with poor mechanical characteristics (Wang and Hao 2002;Yang 2006). In the presence of such conditions, the horizontal component spectra undergoes a deamplification, while the vertical component is generally less affected (Castellaro and Mulargia 2009b).
The presence of dominant sources and other 2D effect could make the HVSR results invalid (Bonnefoy-Claudet et al. 2006;Cara et al. 2010;Benkaci et al. 2021), consequently it is necessary to verify that such sources are not present when the ambient noise is acquired. By evaluating the polarization effects connected with spectral peaks, it is possible to verify stationarity noise and recognize the possible presence of dominant sources. These sources produces marked azimuthal heterogeneity in the characteristics of recorded noise field.
During the first field survey, carried out on 28 August 2018, 156 ambient noise samplings were acquired at the nodes of a regular grid (Fig. 2c) arranged on the area affected by the paroxysmal event of 2008, which still has active emission points. The acquisition grid, with dimension of 110 × 120 m, has an interdistance between measuring points of 10 m. The grid nodes were geo-referenced through a Stop-and-Go GPS survey, using a dual frequency (L1 and L2) GNSS (global navigation system) Topcon HiperSR receiver. The coordinates have an average standard deviation of about 0.025 m for the East component, 0.026 m for the North and 0.064 m for the Up. At each node, the ambient noise was acquired with a sampling frequency of 128 Hz for 16 min.
The ambient noise is characterized by a temporal variability (Cara et al. 2003;Okada 2003;Bonnefoy-Claudet et al. 2006;D'Alessandro et al. 2021) linked to different phenomena depending on the frequencies considered. The noise temporal variability at frequencies below 1 Hz is linked to the variation of meteorological phenomena, at higher frequencies instead, it is influenced by anthropic activity. Fortunately, the latter phenomenon does not seem to affect the area where the samplings were carried out, although it is located not far from the Santa Barbara village. This small neighborhood, located 2 km east of the populous city of Caltanissetta, is characterized only by the presence of some houses, with no industrial or commercial activities in the surroundings.
To ensure that the recorded noise was stationary over time, ten 3D digital portable seismometer were used for the acquisitions which allowed us to acquire all the samplings on the same day during stable weather conditions. In fact, it is known that the amplitude of the noise seems variable for sampling carried out on different days, but it does not vary particularly in a few hours (Okada 2003).
The second field survey was performed after 3 months, on 30 October 2018. In order to acquire an ambient noise that could have characteristics as similar as possible to that acquired previously, it was performed on a day with weather conditions similar to that of the first field survey. The survey was performed sampling the ambient noise on 25 nodes of the acquisition grid located in the central area, where the main emission points are present (Fig. 2c). The acquisition parameters are the same used for the previously field survey.
The records acquired during the two field surveys were processed using the same elaboration steps. The trace related to each acquisition was divided in 20 s time windows, to each window was applied a triangular smoothing function of 10% of the central frequency. The H/V spectrum was calculated in the frequency range 0-64 Hz. The check of SESAME (2005) parameters was performed for all the records. The three parameters for a reliable H/V curve are always satisfied, but not always all six parameters related to clarity of the peak are verified. In particular, the spectral ratio amplitude at peak frequency (f 0 ) does not always have a value greater than 2, the frequency on the H/V curves corresponding to mean + and -one standard deviation is not always located in the frequency range f 0 ± 5%. Finally, in some cases the standard deviation of peak frequency does not stay below a threshold dependent on f 0 , tabulated in the SESAME (2005) guidelines.
Moreover the spectral ratio was calculated along various directions, turning the horizontal components of 10° azimuthal intervals, proceeding from north (0°) to south (180°). This technique allows highlighting the possible presence of dominant sources responsible for polarization effects, which could influence the trend of the H/V spectra.

MASW and joint fit
The MASW method is an active seismic geophysical survey that provides information on shear wave velocity distribution through the study of the dispersion of Rayleigh waves 1 3 (Park et al. 1999;Richart et al. 1970;Achenbach 2000). Two MASW surveys (Fig. 2d) were performed during the first field survey, on 29 August 2018. The survey MSW1 was performed in an area close to that affected by the paroxysmal event of 2008 and which currently has no emission points. This survey was carried out in order to verify the presence of muddy material in the subsoil, outside the area affected by recent deposits. The survey was performed using an array of 25 vertical geophones with natural frequency of 4.5 Hz that were spaced 3 m-apart, for a total length of 72 m.
The survey MSW2 was carried out in the area currently affected by the emission points and also in this case was used an array of 25 geophones, but with a spacing of 2 m, for a total array length of 48 m. For both surveys, the shots were 6 m outside by the first and the last geophone. An 8 kg-hammer was used as a source. To improve the signal to noise ratio (SNR), five-shot vertical stacking was applied. The signal was acquired for 3 s length with a sampling frequency of 256 Hz.
To identify the energy associated with Rayleigh waves, the time series were processed in the frequency-phase velocity domain by applying the Fourier transform. Subsequently, in order to better constraining the 1D shear wave velocity profile and reducing the associated error a joint fit was performed between the phase-velocity dispersion curve of the Rayleigh wave obtained, and the H/V spectrum related to one sampling performed near the MASW alignment (Castellaro and Mulargia 2009b).
Starting from a subsoil model, with assigned values of layer thickness, seismic waves velocity (V p and V s ), Poisson coefficient and density, based on site geological information, a theoretical dispersion curve and a synthetic H/V spectrum are calculated and superimposed to the respective experimental curves. The code, which generates the synthetic H/V curve, uses a forward modeling approach (Castellaro 2016), based on the simulation of the surface-wave (Rayleigh and Love) field in plane-parallel multi-layered systems, according to the theory described in Aki (1964) and Ben-Menahem and Singh (1981). The starting model is then repeatedly modified to obtain the best overlap between theoretical and experimental curves. The model that provides the best fit can be considered the most representative of the seismostratigraphic conditions of the investigated subsoil.

Seismic impedance contrast sections
The integration of the data obtained from the H/V spectra with those relating to the distribution of the shear wave velocity allows to obtain seismo-stratigraphic representations know as impedance contrast sections, which show the distribution of H/V amplitude values in the subsoil.
In an area were no significant 2D variations occur, it is reasonable to hypothesize that the shear waves velocity increases with depth due to the lithostatic load increase, this trend can be modeled according to the following relationship (Ibs-von Seht and Wohlenberg 1999): where V s (z) is shear wave velocity versus depth; V 0 is initial velocity; Z is depth; α is velocity-depth proportionality coefficient.
This function can be used to fit the 1D V s -depth profile in order to obtain the V 0 and α parameters for which there is the minimum misfit between this function and the velocity profile.
The parameters obtained substituted in the Ibs-von Seht and Wohlenberg (1999) formula, allow transforming the H/V spectra frequency values (f) into depth values (h), to which the respective H/V amplitude value is associated: The acquisition of ambient noise samplings at the nodes of a regular grid allows performing an interpolation of the data related to adjacent H/V spectra and consequently to obtain the impedance contrast sections. The interpolation of all the grid data allowed to obtain an impedance contrast 3D model.

2D ERT
Electrical resistivity tomography surveys are intended to determine the distribution of resistivity values below ground level. Resistivity is influenced by various geological parameters, such as mineral and fluid content, porosity and saturation degree of rocks (Keller and Rapolla 1974).
The quality and characteristics of the data obtained depend not only on the properties of the crossed medium, but also by the configuration used for the acquisition.
The resistivity measurements are carried out by injecting an electrical current into the ground through two current electrodes and measuring the resulting potential difference in two potential electrodes. From the values of current (I) and voltage (V), the apparent resistivity (ρ a ) can be calculated using the following equation: where k is the geometric factor that depends on the electrodes configuration.
The best geometry to use is chosen depending on the subsoil structure, instrument sensitivity and degree of background noise present at the investigated site (Loke 2000; (1) V s (z) = V 0 (1 + Z) , Martorana et al. 2009Martorana et al. , 2017Carrara et al. 2015;Corrao and Coco 2021;Patti et al, 2021). After acquiring the data, in order to determine the real resistivity and obtain the electrical resistivity section it is necessary to make an inversion of the apparent resistivity measurements by means of the smoothless-constrained least-squares method (deGroot-Hedlin and ConsTable 1990;Sasaki 1992).
In the study area, two electrical resistivity tomography surveys were performed orthogonal to each other (Fig. 2e), with the crossing point located in correspondence to the area affected by the main emission points. Both arrays were deployed by placing on the ground 48 electrodes, spaced 5 m-apart, for a total length of 235 m.
Since the purpose of the survey was to identify areas characterized by the presence of the emitted material, a Wenner-Alfa configuration was used for data acquisition, which is suitable for highlighting vertical variations and is characterized by a good signal to noise ratio.
The Schwarz-Christoffel transformation method was used to calculate the distortion in the surface levels, linked to the topography. Data inversion was performed with the smoothless-constrained least-squares method with quasi-Newton optimization (Loke and Barker 1996). Five iterations were performed for processing the ERT1 and eight for the ERT2, with an associate RMS of 1.8% and 2.7%, respectively.

HVSR and MASW
The processing of ambient noise samplings, with the HVSR technique, has no highlighted significant resonance peaks (SESAME 2005). However, a seismo-stratigraphic interpretation can be made anyway by observing the single component Fourier spectra. Indeed, it is possible to notice local minima in the vertical component, which produce the typical 'eye-shape', indicating the presence of a seismo-stratigraphic transition between materials that have different mechanical characteristics (Castellaro and Mulargia 2009b;Molnar et al. 2018;Xu and Wang 2021). In some spectra, resonance peaks with greater amplitude (2-4) can be observed at higher frequencies (> 20 Hz). These higher frequencies peaks suggest existence of an impedance contrast at shallow depth.
The obtained HVSR spectra (Fig. 3a) have highlighted, in almost all acquisitions made during the two field surveys, a deamplification under the amplitude value 1 for a wide range of frequencies, on average between 1 and 20 Hz. The related single component spectra (Fig. 3b) show indeed the decrease of the horizontal components below the vertical one. This phenomenon attributable to the presence of a velocity inversion in the subsoil was carefully studied to identify the areas that could be characterized by the presence of fluid-rich muddy material emitted on the surface. Previous studies have observed that the saturation increase in clayey sediments causes a reduction in shear wave velocity (Xu and White 1996;Carrière et al. 2018). Therefore, it is reasonable to assume that this material, due to its more fluid-rich composition than the surrounding lithotype, has poorer physical-mechanical characteristics, leading to a consequent reduction in shear wave velocity within it.
The directional analysis performed did not show the presence of strong directional effects for both field surveys. The results are shown through polar diagrams where H/V amplitude values are plotted for different directions (Fig. 3c). The polar diagrams highlight no significant deviations from the isotropy condition for both field surveys. It is, therefore, possible to exclude the presence of dominant sources that may influence the results of ambient noise records. The spectra characteristics are ascribable only to the actual subsoil conditions.
The 1D V s -depth profiles obtained, through MASW surveys processing, show in both cases low shear wave velocity values according to the outcrop lithotype. In particular, the 1D V s -depth profile, related to the MSW1 survey (Fig. 4a,  b), performed in an area which currently has no emission points, it is characterized up to about 2.5 m of depth by V s values that are around 125-129 m/s. Under this layer, the V s value tend to increase slightly, reaching a value equal to 186 m/s in the last seismic layer identified at depths greater than 21 m. The small increase in the shear wave velocity does not indicate the presence of stratigraphic variations, it is attributable to increased compaction of the material as a result of increasing lithostatic loading.
The 1D V s -depth profile related to MSW2 alignment, carried out within the area that has active emission points, was obtained through the HVSR-MASW joint fit technique (Fig. 4c-e). This technique allows to overcome the problem linked with the limited exploration depth of the MASW survey. Since the HVSR surveys analyse seismic waves with larger wavelengths, which enables investigating greater depths than those achieved with only MASW surveys, albeit with less accuracy. Therefore, the H/V inversion extends the V s profiles to greater depths. The 1D V s -depth profile shows near the surface a low velocity layer with a thickness of 4 m, probably due to the presence of very saturated clays. Below this layer, the shear wave velocity tend to increase slowly reaching a value of about 160 m/s at depths greater than 42 m.
Since the subsoil of the study area is characterized by a single lithotype, the lower velocities values observed in the 1D V s -depth profile relating to the MSW2 alignment confirm that the survey involves an area where muddy material is present in the subsoil (Fig. 4f). Given the small velocity inversion (the velocity is reduced by only 5 m/s) that affects a limited thickness (only 4 m), it is still possible to consider a velocity average trend that increases with depth.
A fit of the latter V s -depth profile was carried out using the function (1) in order to obtaining the V 0 and α parameters necessary to convert the frequency values of H/V spectra into depth values (Fig. 4e). The fit was performed on the V s -depth profile of MSW2, because it is representative of the subsoil conditions in the area where the ambient noise samplings were carried out. The minimum misfit between 1 3

Impedance contrast sections and 3D model
The parameters obtained by fitting the Vs-depth profile with the function (1) were used in the Eq. (2) to convert all the frequency values of the H/V spectra, from both field surveys, to depth values. In this way, each H/V amplitude value was associated with a depth, in order to show the distribution of seismic impedance contrasts in the subsoil. It is reasonable to believe that these parameters are not subject to significant changes between the two field surveys. In fact, no significant paroxysmal events were recorded between the two field surveys. In addition, low rainfall cannot affect the subsurface distribution of the areas affected by the presence of muddy material, since the lithotype on site is impermeable. It must be considered that the 1D V s -depth profile, obtained from a MASW survey, is the average result of the V s velocity distribution beneath the array (Park et al. 1999). Consequently, if no significant changes occur in the subsoil, this profile should not changes significantly in a short time like that elapsed between the two field surveys.
Interpolating the data related to samplings carried out along the same alignment with regular interdistance, the impedance contrast sections were reconstructed. These sections show on the x-axis the progressive distance of the measurement points, on the y-axis the depth, and the plotted value corresponds to H/V amplitude; variations in this value are associated with a color scale.
From the 156 ambient noise samplings, performed during the first field survey, 25 impedance contrast sections were obtained, 12 in SSW-NNE direction and 13 in WNW-ESE direction (Fig. 5a, b).
All the impedance contrast sections are characterized by low H/V amplitude values (< 1) in their shallow parts up to a depth of about 5 m from the ground level. Another area with same characteristics can be observed at a depth of approximately 40 m in the northern portion of the acquisition grid, were the main emission point are located. This area, with an average thickness of about 20 m, decreases in thickness in the eastern portion of the acquisition grid until it almost disappears. In particular, in the sections oriented WNW-ESE this area is observable, only in those located further north (G-H-I-L-M-N-O), from 0 to 90 m to the origin of the profiles. In the SSW-NNE sections, this area is observable from 60 m to the end of the profiles. It tends to thin in the sections located further east (10-11-12).
By interpolating the data of the impedance contrast sections, it was possible to reconstruct a 3D model, showing the volumetric extension of the areas characterized by different H/V amplitude values (Fig. 6a, b). The model highlights in the northern portion a connection zone between the two areas characterized by H/V amplitude values lower than 1. The area located at shallower depth seems to have a greater extension in the horizontal direction than the one covered by the sampling points. The second area, deeper, appears to be confined to the southern sector of the acquisition grid and extends beyond the area investigated in the northern, western and eastern sectors.
Also from the measurements carried out during the second field survey, 10 impedance contrast sections were obtained, 5 in SSW-NNE direction and 5 in WNW-ESE direction. These sections confirm the results of the previous field survey, as also in this case the impedance contrast sections are characterized by two areas with H/V amplitude value lower than 1. The first just below the surface with an average thickness of about 5 m and the second deeper at about 40 m of depth with greater thickness. Although having on average the same characteristics as the areas observed previously, if a comparison is made between homologous sections it is possible to find differences in the shape, size and position of these areas. The 3D model obtained by interpolating the data from the impedance contrast sections (Fig. 7a-c) confirms the presence of a connection zone between the two areas, located in the NW corner of the model. This area is located approximately in the same position as that identified in the model obtained from the previous field survey.

2D ERT
The electrical resistivity tomographies obtained are characterized by low resistivity values due to the presence of clayey lithotype. However, it is possible to observe areas characterized by particularly low values that may be associated with the stationing areas or the ascent path of emitted material.
Both the electrical resistivity tomographies have reached an investigation depth of about 40 m.
In the ERT1 section (Fig. 8a), the lower resistivity values (< 0.5 Ω* m) are observable in the center of the alignment, between 100 and 140 m from the origin. This low resistivity area extends from 10 m of depth up to the maximum investigation depth. Two other shallow areas with lower resistivity, characterized by a consistent lateral extension but little vertical development, can be observed south and north of the previously identified area. Furthermore, the section shows two areas with higher resistivity, one located at the roof of the area with low resistivity, placed in the center of the profile, and another at the end of the section. These areas located near the surface have a moderate lateral extension and a very small vertical extension.
Even the ERT2 section (Fig. 8b) shows an area characterized by lower resistivity values at the center of the profile, which extends up to the maximum investigated depth. Another significant area with lower resistivity value is observable at the end of the section. Three areas with greater resistivity values, characterized by a limited horizontal and vertical extension, are observable at the beginning, in the middle and at the end of the profile.
In order to better correlate and interpret the ERT results, the sections obtained were visualized in a three-dimensional environment, and a Digital Elevation Model (DEM) of the area was superimposed on them (Fig. 8c). It was, thus, possible to note immediately that the area characterized by low resistivity values, which extends to the maximum investigated depth, is located at the crossing point of the two ERT surveys and is placed in correspondence with the area characterized by the presence of the main emission points.

Discussion
The processing, interpretation and integration of data obtained from different geophysical surveys allowed us to obtain information on the likely ascend path and accumulation zones of the emitted material in the Santa Barbara mud volcano area.
The spectra obtained from the processing of ambient noise samplings do not show significant resonance peaks (Fig. 3a). Considering the geological framework, this data are in accordance with the subsoil characteristics. In fact, in Fig. 9 Comparison between the isosurfaces that enclose areas characterized by H/V amplitude values lower than 1 related to the 3D models obtained from the two field surveys: perspective (a) and lateral (b, c) views; the box at the top left shows the point of view the area there is no a lithologic succession, which can produce resonance peaks, as the subsoil is characterized by the presence of a single lithotype, the clays. However, observing the single components Fourier spectra (Fig. 3b), it is possible to identify the presence of resonance peaks, even if with H/V amplitude values are lower than 2. These peaks are probably linked to the seismostratigraphic transition between portions of the same lithotype with different mechanical characteristics. In order to exclude the possibility that the data obtained from ambient noise samplings could be influenced by the presence of dominant sources that a study on the possible connection of polarization effects with spectral peaks (Fig. 3c) was carried out. The maintenance of the noise isotropy conditions in both field surveys, allowed us to exclude the influence of dominant sources on the acquired data. This condition have encouraged and allowed a comparison between the results of the two field surveys in an attempt to monitor the variations that occur over time.
As shown in the V s -depth profiles (Fig. 4b, e), the study area is characterized by very low shear wave velocities. By integrating the data of the H/V spectra with the information concerning the trend of the V s values with depth, the impedance contrast sections were obtained (Fig. 5).
The first meters of depth of all impedance contrast sections are characterized by low H/V amplitude values (< 1), these areas could be associated with the presence of muddy material. A valid explanation on the origin of these areas could be the overlap of mudflow not fully consolidated yet. Another area with the same H/V amplitude value is observable at a depth of about 40 m in the north-west sector of the acquisition grid. This area, with an average thickness of about 20 m, highlights the possible presence of another zone affected by muddy material. The parts of sections characterized by major H/V amplitude values can be related with the transition between the muddy material and the surrounding lithotype. The 3D seismostratigraphic model (Fig. 6) allowed us to highlight the volumetric extensions of the two area affected by low H/V amplitude value (< 1). The 3D model highlighted the presence of a connection zone between the two areas that were located in the north-west sector. This may be associated with the likely ascent path of emitted material. The presence of this connection area is also highlighted, with approximately the same location and extension, in the 3D model obtained from the data acquired in the second field survey (Fig. 7). Even in this case, two areas characterized by H/V amplitude value lower than 1 were detected. These zones may be associated with those previously identified because they have on average the same location. However, comparing the surfaces that circumscribe these volumes (Fig. 9a-c, Online Resource 1) it can be possible to note a variation in shape and size. In particular, in the 3D model related to the second field survey a reduction in the volume of these areas is observable and an apparent migration of the same toward slightly lower depths. The lack of recharge of emitted material from higher depths could be, for example, the cause of such changes.
In addition, two ERT surveys were performed in the Santa Barbara mud volcano area (Fig. 8) in order to compare and validate the results obtained. The acquisition of these data in the considered context presents some difficulties linked to the characteristics of the areas affected by emission points, the most significant of which is the high imbibition with saline fluids of the surface levels. This could cause some current leakage phenomena and produce unreliable results. Therefore, the data obtained were evaluated with the appropriate caution and compared with those obtained from surface observation (location of emission points) and other geophysical surveys (Fig. 10 a, b).
The ERT alignments, orthogonal to each other, have a greater extension than the area affected by the deposits of the last paroxysmal event. The crossing point between the two profiles is located in correspondence with the area where the main emission points are present.
In this sector, both the geoelectric sections identify an area characterized by particularly low resistivity values that extends up to the maximum depth investigated. Its location corresponds to the connection area previously identified in the 3D models that can be associated with the likely ascent path of the emitted material (Fig. 10 a, b). The fluid-rich muddy material that rises within the surrounding lithotype should be the cause of the very low resistivity values observed in this area. Areas of higher resistivity can be observed in the profiles edges, near the surface, these identify the mud volcano peripheral areas currently covered by vegetation and not affected by recent mudflows. Other low resistivity areas are observable in the ERTs edges, which could be interpreted with the presence of accumulation areas of the emitted material, located at low depth. Their feeding area is not visible in the tomographic sections. These areas could be fed by subsurface channels within which the fluid-rich muddy materials migrate to different areas of the volcanic system, acting as feeding zone for potential emissions located in peripheral areas.

Conclusions
The activity in the Santa Barbara mud volcano area was documented since 1800. The last intense paroxystic event, occurred in 2008, has caused damage to edifices and infrastructures located in the surrounding area. When these geological structures are located near inhabited centers, it is essential to carry out studies aimed at a more complete characterization of the processes that cause the paroxysmal events in order to safeguard the population. In this work, the applicability of some survey techniques rarely used up to now in the study of mud volcanoes has been evaluated.
The study of mud volcanoes through the application of geophysical techniques may prove useful for the evaluation of the characteristics and the detailed structure of the emission system in the first tens of meters below the surface. This aspect is often difficult to determine with survey tools generally applied to the phenomenon. The acquired information can be useful for assessing the surrounding conditions that govern the short scale activity of a site and consequently contribute to the mitigation of the risk associated with them.
The integration and joint interpretation of the acquired data made it possible to mutually validate the results and obtain a model that most likely corresponds to the real conditions that characterize the investigated site.
The variations identified by comparing the 3D models, obtained by the two field surveys carried out three months apart suggest, as indeed it is reasonable to expect that the investigated area is subject to continuous changes even after short periods of time. The possibility to monitor these changes over time could prove to be a useful tool for risk mitigation in the areas affected by the presence of mud volcanoes.
The results obtained suggest the aptitude of the methodological approach employed for use in monitoring changes in the subsoil characteristics in areas affected by these phenomena. However, we believe that further tests and verifications are necessary to fully ascertain the validity of the method in monitoring studies.
From this study, it emerges that the integration of the results obtained from the applied geophysical methodologies can have good potential in the study of the mud volcanoes. Indeed, it was possible define the migration pathway and accumulation zones of fluid-rich muddy material in the first few tens of meters of depth. The detection of buried accumulation zones may be a clue for predicting the possible development of new emission points, or for assessing changes in emissions from existing ones.
Author contribution SG and SI conceived and designed the experiments; SG, GP, GDeG, FB, FC and SI performed the experiments and analyzed the data; SG, GP and SI wrote the paper.

Conflict of interest
The authors have no financial or proprietary interests in any material discussed in this article.
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 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/.