The pre-processing phase (Ismail et al., 2020a, Figure 3) involves the preparation processes that are being used to sweep up both major and relatively minor issues in the raw data, such as the core processing flow (i.e., velocity analysis, frequency filtering, post-stack migration, etc). Even though the pre-processing processes appear to be common and even regular, each step is crucial to the final results (Zhou, 2014). The techniques that resulted in the final seismic image include a standard methodology, as shown in the workflow in Figure 3.

Pre-processing steps include editing and muting, which were required to adapt the dataset to the main processing steps (Figure 3). During the pre-processing stage, the key challenge has been identifying the reflection phases, which are obscured by the very high amplitude of ground roll and very close in time to the refraction waves. By recovering the amplitude of the weaker reflection phases, Automatic Gain Control (AGC) scaling has solved the amplitude problem. In the L1 norm, a time gate of 250 ms was utilized for AGC scaling. The consistency and coherence of seismic reflectors have significantly improved the AGC output. Then, to reduce high- and low-frequency noise, the band-pass Ormsby filter was utilized.

The pre-stack stage is the core of the whole procedure (Figure 3); it starts converting the seismic data from the time to the frequency domain. Subsequently, the F-K filter is used to remove the ground roll noise that can be easily detected as it corresponds to the low frequencies in the filtered signal. According to Serra et al. (2016), the frequency-wavenumber spectrum has been analyzed to find a correlation between the frequency content and the velocity of ground roll phases. Then, the F-K filter was constructed to isolate and eliminate the ground roll phases. The reflected waves lost a high percentage of energy due to the diffraction, scattering, and absorption during the propagation. In particular, because of the complexity of the subsurface, the deeper reflectors from the two-way time (TWT) of 500 ms did not show any coherency in the semblance velocity analysis, therefore can not be picked through the Normal-Move-Out (NMO) correction. For all these reasons, just the first 500 ms of seismic data were used in the subsequent processing steps.

The arrival times to the receiver depend on the position of the source and their offset and on the velocity model, which is an unknown parameter. The evaluation of the velocity model is an inverse problem that can be solved using several methods. In this work, the NMO velocity method has been considered and the velocity spectrum has been computed to determine the stacking velocity of each reflector.

The main goal of the velocity analysis (Figure 4) is to provide a multi-1D (or pseudo-2D) velocity model (Figure 5) profile to be used for an optimal Common-Mid-Point (CMP) stacking section with the identification of all reflection events in the correct TWT positions which allows enhancing their amplitudes and signal-to-noise ratio in the stacking section. In Figure 5, the velocity model used for profile A obtained from the velocity analysis is shown. The velocity interval for the NMO analysis has been inferred from previous works (Gammaldi et al., 2018; Ismail et al., 2020a) and has been set to 200–2500 m/s (Supplementary Figure S1).

The previous processing and comparison steps are needed to help us improve and test the semblance plot to provide a good estimation of the final stacking velocity model. This is the main reason behind the good quality of the obtained Common-Depth-Point (CDP) stacking that lets us identify clear horizons displaying the disrupted beds and faults related to collapsed sediments in the post-migration profiles.

The time shifts related to the source-receiver position at the surface are evaluated using a combination of cross-correlation technique and stack power optimization in residual static corrections. By measuring the time shifts between each trace within a CMP gather and a “pilot” trace, the selected times are then decomposed to solve for the source and receiver statics in a surface-consistent approach. The main objective of stack power enhancement is to highlight the effectiveness of a stacked trace along the traces in a CDP gather. For a given source and receiver location, a variety of time shifts are tested, and the one that generates the highest stack output is selected as a surface-consistent residual static value for that site.

Finally, to recover events with conflicting dips during CMP stacking, Dip-Moveout correction (DMO) is required to adjust the dip impact affecting stacking velocities (Yilmaz, 2001). After NMO correction, a dip-dependent partial migration is applied to nonzero offset seismic profiles to convert them from non-hyperbolic offsets to zero offsets. This shift from non-hyperbolic to hyperbolic (zero) offset results in enhanced velocity estimations and lateral resolution, as well as a few other useful side effects including coherence noise attenuation (Bruno et al., 2002). CMP stacking is composed of three stages: CMP binning, NMO correction, and NMO adjusted traces stacking (Zhou, 2014). The maximum offset covered in each CMP gather, which is the length of the active spread, is nearly equivalent to 100 m, which would be the determined maximum depth in the multi-2D image. An interface’s seismic identification refers to whether or not there is a reflection from the interface that can be detected in seismic data. This is a far more challenging task with seismic resolution since it involves many additional variables, some of which are poorly constrained and challenging to be analyzed (Magee et al., 2018).

Seismic migration is a technique for removing distortions from reflection records by relocating events to their proper spatial position and reducing the energy of scattering backward toward the dispersing spots (Schneider 1978; Bleistein 1987; Gray et al., 2001; Yilmaz, 2001). Post-stack Kirchhoff migration is used as the final processing step in the processing flowchart after relocating the seismic events geometrically in both time and space (Gray et al., 2001; Ismail et al., 2020a) and before interpreting all features in the seismic profile A.

The subsurface of the volcanic area is characterized by complex geology, such as faults and collapsed areas as already shown in the SV (Ismail et al., 2020a). The Kirchhoff migration approach was employed in both time and space. After resolving diffraction patterns, Kirchhoff migration shifts dipping reflectors correctly to their corrected and real subsurface position (Figure 6). Kirchhoff migration has many benefits over other approaches, including its flexibility and ability to manage lateral velocity fluctuations well.

By converting the vertical axis of the stacked data from TWT to the depth domain via the use of the final stacking velocities image built during the CDP processing, the Time-Depth (T-D) conversion has indeed been crucial in obtaining the position of all features and faults beneath the surface in depth coordinate (Supplementary Figure S2).

The processing steps and signal enhancement have been performed to generate the 2D velocity profile (Figure 5) used for NMO correction as the main processing step to produce the angle and full-stack profiles (Supplementary Figures S1, S2). Moreover, the generated velocity profile played a significant role in relocating the reflection events to their correct time and space by using the post-stack Kirchhoff migration and producing the final 2D TWT WNW-ESE seismic profile (profile A) (Figure 7) and profile B (Supplementary Figure S4). Then, time to depth conversion is performed to the final 2D profile by using the generated 2D velocity file to obtain a better geometrical representation for all seismic features including seismic horizons and faults (Figure 7B) where the profile extends around 400 m as a total distance and around 150 m as total depth. Several classes of post-stack seismic attributes (Figure 7) are extracted to delineate more information from the seismic profiles. The time gain attribute (Figure 7B) extracted for the 2D depth profile to improve the display, enhances the amplitude of deep horizons and gives a clear image of the structural frame in the profile which shows the near-surface layer from the surface to 50 m depth. Furthermore, the near-surface depth interval (0–50 m) shows horizons discontinuities in the profile, with most faults extending from the profile’s highest depth to the extensively deformed near-surface horizons, particularly between CMPs 22 and 42 (Figure 7).

The maximum values of amplitude are exhibited between 30 and 50 m in depth. Therefore, both profiles A and B show seismic reflections on the east side between 30-and 50 m. Hence, the WNW-ESE seismic profile (A) describes three main depth intervals: in the near-surface area at the depth range from 0 to 30 m bending reflectors with a damaged zone in the central part are shown; going deeper into the profiles, the depth ranges from 30 to 50 m where the area is characterized by the presence of bright anomalies (Figure 7B). Finally, at the depth range between 50 and 150 m, there are poorly organized reflections affected by discontinuities associated with more faults extending between CDP 20 and CDP 60 which play a significant role in the fluids uprising to the near-surface depths and are characterized by poorly organized pattern reflections with amplitude attenuation in the faulted zones.

The main target of using unconventional seismic techniques is to improve the imaging of different geo-bodies and features such as gas-bearing zones, bright spots, faults, and acoustic blanking in both WNW-ESE (Profile A) and NNE-SSW (Profile B) seismic profiles.

Seismic attributes analysis started by extracting different classes of post-stack seismic attributes such as time gain, envelope, RMS, and energy attributes to obtain the best results much more easily and quickly. Due to the abrupt change in reflectivity over the shallower part of the profiles, the detected bright anomalies in the original seismic profile are also characterized by RMS, high Energy, and envelope values (Figures 7A,C,D). Furthermore, in Figure 7C, the energy attribute displays and distinguishes the maximum amplitude values throughout a depth range of 20–50, with the maximum energy range depicted in red. Finally, we present poorly organized reflections impacted by discontinuities associated with more faults located between CDP 20 and CDP 60, as evaluated after extracting RMS, time gain, and envelope attributes at depths between 50 and 150 m (Figures 7A–C) which play a significant role in the fluids uprising to the near-surface depths.

AVO analysis provides more information from the pre-stack seismic dataset. It is used to produce a new vision for identifying gas zones, using pre-stack AVO attributes based on the sPr changes. AVO analysis is performed by obtaining specific AVO attributes, which include intercept (A) and gradient (B) (Supplementary Figure S3). This step lets us implement the sPr change (Figure 8) and generate cross plots to distinguish the gas-bearing zone. The performed cross plots of the intercept and gradient for all time samples throughout all trace positions within the same setting window provide a more precise solution (Supplementary Figures S4, S5). Cross plots discriminate between gas zones (white color) and background zone (red color). Therefore, the gas zone is separated and distinguished for profile B in Figure 8 at TWT range of 80–120 ms between CDPs 21 to 55 and in the central part of profile A between CDPs 45 to 75 and TWT range of 80–140 ms (Figure 8).

To achieve a final interpretation of the shallower structures and fluids interactions using all the images, several multi-2D views are performed.

Merging profiles A and B, the first multi-2D image obtained using 3DMove, in Figure 9A shows the seismic sections and the interpretation obtained by means of the seismic attributes (Supplementary Figure S4).

To interpret the main bright spots and the structures, together with the AVO analysis, a prior multi-2D interpretation of the AVO section is needed (Supplementary Figure S7). According to the sPr shown in Figure 8, superimposed in Figure 9B, it is possible to note that the fluids accumulation zone is mostly located in the north-east area. Finally, the interpretation of the sPr, according to the electrical resistivity model beneath the fumarolic area in Gresse et al. (2018) and their interpretation of the gas reservoir (a mixture of CO₂ and H₂O) and condensate flow, shows the contact between the GDF and LDF in the first 50 m of the SV subsoil structure.

In the images, the structural evidence performed by seismic attributes in the depth domain is plotted as red dotted lines (interpreted as faults) and black dotted lines (interpreted as chaotic zones) (see Figures 9A,B). The time evidence of the fluids is now superimposed in Figure 9B in terms of continuities (GDF and LDF) with the black arrows pointing to the GDF and LDF continuities contact. The joint multi-2D image shows the results obtained from different analyses which merge at the precision of meters. Moreover, the low-velocity anomalies and the interpretation in Figure 2 of the water and gas-saturated layers are compatible with the faults and the GDF and LDF extracted in this work.

These anomalies, together with the previous research performed in the SV, strongly suggest that the faulted zone represents the pathway for the uprising gases coming from the hydrothermal plume in the Solfatara Volcano testified in the shallower images by the chaotic zones. Then, in the shallower area of about 40–50 m, the gases appear to be trapped. As shown in Figure 9B, the continuities are well enhanced by the AVO analysis because of the contrast between the liquids located in the area. In addition, in this zone, liquids have been already interpreted and are very well correlated with the Vp velocity by Gammaldi et al. (2018), shown in Figure 2, and the reflections displayed in Figure 9A for profile A and Ismail et al. (2020a).

The updated multi-2D image demonstrates that the fault structures here identified are probably working as a pathway for fluids from the deep hydrothermal plume to the surface. Meanwhile, for the first time, AVO analysis highlights and distinguishes the shallower interplay between gas and condensate flow in the north-east area of the crater. This contact probably causes a geochemical alteration in the shallow unconsolidated tephra (Figure 2). According to Mayer et al. (2016), the interaction shown in Figure 9 with a resolution of meters may be caused by the presence of the Alunite in the formation due to the hydrothermal alteration between GDF and the LDF. Therefore, the alunite represents the cap-rock where the fluids, from the deeper hydrothermal plume through the faults and the chaotic zones identified by the seismic, are trapped. It is reasonable to conclude that, over time, the continuous charging of the hydrothermal plume was able to form the cap-rock, which is constituted by the alunite, using the hydrothermal alteration of superficial rocks between the condensate flow and the uprising gases.