The Paola Basin lies at the rear of the Calabrian Arc, in the upper plate of the northwestward dipping Tyrrhenian-Ionian subduction system, which includes the back-arc domain of the SE Tyrrhenian Sea, volcanic seamounts, the Aeolian Islands volcanic arc and the fore-arc region of the Calabrian accretionary wedge (Figure 1A). The basin formed since the late Miocene, as a result of compressional (Pepe et al., 2010) and strike-slip (Milia et al., 2009) tectonics, in an overall extensional regime.

The Paola Basin area is flanked to the north and west by active and inactive volcanoes. In the north, the area is characterized by the presence of magmatic intrusions that locally reach the seafloor forming the Diamante and Enotrio volcanic edifices and the flat-topped Ovidio Seamount (Würtz and Rovere, 2015, Figure 1B). Altogether, these form a volcanic-intrusive complex characterized by a deep-rooted, magma feeding system, formed in consequence of the ascent of subduction-induced mantle flow originated in the northwestern edge of the retreating ionian slab (De Ritis et al., 2019). The volcanoes are not presently active even if hydrothermal activity occurs, as the emplacement and cooling of the magma occurred sometime during the Brunhes Chron (De Ritis et al., 2019).

The Alcione Seamount is a ∼1,000 m-high conical volcano located on an almost flat seafloor on the lower slope of the Paola Basin (Marani and Gamberi, 2004) and is characterized by absence or very weak hydrothermal activity (Lupton et al., 2011). The Lametini Seamounts are two conical edifices on top of which Fe-Mn crusts characterized by a very low Fe/Mn ratio and high Cu-content indicative of hydrothermal venting were dredged (Rossi et al., 1980) (Figure 1B).

A 60-km-long anticline (Paola Ridge, Figure 1A) bounds westward the Paola Basin syncline (Gamberi and Rovere, 2010). Here, cold seep-like structures were discovered with full-ocean depth multibeam and backscatter data at 500–800 m water depth. Higher resolution acquisitions and seafloor sampling were carried out in 2011 (Rovere et al., 2014). The most prominent structures in the area are represented by three mound structures (RMV, MMV and R1MV, Figure 1B) which make up the Paola Ridge and stand 200–400 m above the Propeller Basin located landward (Figure 1B). Active gas venting at the seafloor was detected via multibeam echosounder at RMV location (Rovere et al., 2014), while oxy-hydroxides (goethite), sulfides (pyrite, marcasite and sphalerite) and siderites were collected in the sub-seafloor of all mounds (Rovere et al., 2015). Siderites precipitated in burrows and hardgrounds in the shallower sedimentary succession showing enrichment in δ¹³C (3–10 V-PDB ‰) and δ¹⁸O (8–9 V-PDB ‰) isotopes, which are compatible with their precipitation in the methanogenic zone (Rovere et al., 2015) during periods of lower gas discharges under prolonged anoxic conditions and intermittent venting of deep-sourced CO₂ (Franchi et al., 2017).

NW-SE and NNW-SSE-oriented normal faults, which can be regarded as the marine prolongation of the fault zones that dissect the Calabrian Arc from SE to NW, were suggested to be the most probable mechanism for the emplacement of the mud mobile structures and primary conduits for upward fluid migration (Rovere et al., 2014). It was further noted that the mound structures rise in coincidence with extensional faults that offset the Messinian evaporites, suggesting that pre-Messinian source rocks, mobilized along discrete belts of active tectonic deformation, were controlling the seepage of fluids in the study area (Gamberi and Rovere, 2010).

The acoustic detection of the water column flares was performed using Kongsberg EM710 (100 kHz, 1° × 1°) and EM302 (30 kHz, 1° × 2°) multibeam systems (MBES) that were both hull-mounted on board the R/V Urania. Bathymetric data were also collected and integrated the seafloor coverage achieved during the previous campaign carried out in 2011 (Rovere et al., 2014). Sound velocity profiles were obtained with a Seabird-Scientific SBE 911 plus CTD unit mounted on a rosette carrying 24 Niskin bottles and were applied both in real-time acquisition and post-processing. Bathymetric and seafloor reflectivity data were post-processed using the suite CARIS HIPS & SIPS.

Water column data were extracted mainly from EM302 data after quality check and processed with the QPS Fledermaus (version 7.8, including the FMMidwater module) software package which allows for manual flare identification and geo-picking of flare sources. A normalization filter was applied to identify background noise values (dB) and apply signal-to-noise level thresholds (see Rovere et al., 2020 for a complete description of the processing routine). Water column plume features (flares) were manually extracted and exported as ASCII points, including longitude, latitude, depth of the feature and signal amplitude (dB) corrected with sound velocity profiles. The georeferenced seep data were imported into the Fledermaus software to explore virtual 3D scenes containing digital elevation models of the seafloor and sub-seafloor (Figure 1C). Geographical visualization and statistical analysis were performed using ESRI ArcMap 10.8. The relative height of the flares was obtained using the spatial analyst tool “extract multi value points” applied to the raster obtained from the bathymetric data acquired simultaneously with the water column.

Bottom waters were collected with a box-corer from the two most energetic gas discharging locations (R1MV and RMV), one of which already identified and imaged (Rovere et al., 2014). Water samples were taken immediately and bubble free after retrieval using a silicon tube. Samples were directly transferred into 250 ml pre-evacuated glass flasks and were acidified with 2 ml 37% HCl. The bottles were sealed with a Teflon-coated butyl rubber seal and were closed with aluminum crimp caps. Five additional near-bottom water samples were collected by means of the rosette sampler alongside the two most energetic structures. Dissolved gas concentrations were determined applying headspace analysis and gas chromatography techniques described in Franchi et al. (2017).

Multichannel seismic reflection profiles (MCS, Figure 2A) were acquired by a single 60 in³ mini Sercel GI-gun set in Harmonic Mode (30 Generator +30 Injector) with a shot interval of 9.375 m at 2000 psi. The data were collected by a 300 m long, 96 channels Geometrics GeoEel digital streamer with a trace distance of 3.125 m. Both the gun and the streamer were towed at 1.5 m below sea level to minimize the ghost effect on the emitted spectrum (centered well above 200 Hz), thus preserving the high frequencies needed to better resolve the shallower targets. During the acquisition four Digicourse birds were used to keep the streamer at a constant depth. The acquisition setting allowed a compromise between penetration in 700–800 m water depth and resolution, which can be regarded as metric along the vertical axes, if assuming a λ/4 criterion, with an actual 1.56 m lateral distance between the traces in the stack section. With this configuration, resulting in an effective horizontal sampling of 1.56 m in the stacked section, the maximum attainable fold coverage was 16 traces/CDP. The polarity of the seismic data was recorded as minimum-phase reverse SEG Standard. A conventional processing was applied by means of the Schlumberger Vista package; the sequence consisted of shot gather editing, large band pass filtering, velocity analyses, velocity based spherical divergence compensation, NMO correction, stack in common midpoint domain. The stack sections were eventually time migrated.

MCS profiles were further treated for image visualization with the Geosuite AllWorks software in which an A.G.C. window length of 250 ms was added to enhance the deeper reflections. The Petrel E&P Software Platform was used to extract the Root Mean Square (RMS) amplitude (iterative) attribute and to detect bright spots.

Vintage sparker profiles were made available by the Institute of Marine Sciences as paper profiles that were scanned to high resolution raster image (TIFF) and converted into georeferenced SEGY format using the Matlab routine IMAGE2SEGY (Farran, 2008).

The reflector corresponding to the top of the Messinian units (MES) was mapped on available MCS and sparker profiles and interpolated with spherical kriging with 1 km radius; contours were extracted as isobaths expressed in ms TWT (Figure 2B).

The seismic data are presented in TWT, as an accurate depth conversion is not feasible due to the lack of boreholes and a proper velocity model. However, some features are described in distance-depth and the thickness of sediment and volcanic deposits were derived from time to depth conversion using velocities of 1,500 and 2,100 m/s for respectively the water column and the Plio-Quaternary deposits, for consistency with De Ritis et al. (2019).

Sub-bottom echosounder profiles were used to image shallow sedimentary structures and map gas indicators with a Teledyne Benthos CHIRP-III system. The system is comprised of a 16 hull-mounted transducer array with a 2–20 kHz sweep-modulated bandwidth and 4 kW power per-channel, which allows a vertical resolution of about 50 cm. Data were post-processed for trace equalization and band-pass filtering with the software Geosuite All Works. Calculations for converting from two-way-traveltime to depth used a sound velocity of 1,500 m/s.

In seismic reflection data, bright spots and high-amplitude anomalies, seismic attenuation, velocity pull-down events, flat spots and chimneys are indicative of the presence of gas in the sediment and fluid contact (Müller et al., 2018). In our seismic data, gas-induced bright spots have the characteristic “peak-over-trough” pattern, in which a high-amplitude trough follows a high amplitude peak (a decrease in acoustic impedance, negative amplitude anomaly). On the contrary, the trough-over-peak reflections or positive amplitude anomalies represent a transition to a harder material (increase in acoustic impedance) generally associated with the occurrence of denser sediments, hardgrounds, precipitated authigenic carbonates, gas hydrates or the seabed.

Furthermore, in seep environments pull-up events may be caused by the presence of rocks with higher seismic velocity than the strata above, such as authigenic carbonates (Madof, 2018).

Magmatic sills can be interpreted based on their high amplitude character, showing as localized brightening of positive amplitude reflections similar to the seabed reflector (Omosanya et al., 2018). In terms of geometry and lateral continuity within the host strata, sills are imaged in seismic reflection data as saucer-shaped or sheet-like, high-amplitude reflections with abrupt lateral terminations within the host-rock strata (Roelofse et al., 2021).

Hydrothermal vent complexes are imaged on seismic sections as pipe-like, vertical zones of low amplitude and chaotic reflections in the conduits, terminating as dome, eye-shaped or crater morphologies at their summits (Omosanya et al., 2018). The interior of the vents vary between chaotic seismic reflections and clear reflections that terminate within the vents as a result of the disruption of the originally stratified sedimentary rock during fluid expulsion (Roelofse et al., 2021).

In side-scan sonar and multibeam backscatter data, the recognition of seep sites on the seafloor is often favoured by the detection of anomalously high acoustic seafloor reflectivity, caused by enhanced acoustic impedance contrast related to: sharp changes in the seafloor morphology (pockmarks, mud volcanoes); precipitation of authigenic carbonates at the seafloor; bubbles or gas hydrates in the sediment (Naudts et al., 2008).

The seismic stratigraphic analysis allowed the reconstruction of the depositional architecture of the stratigraphic units of the late Miocene and Pliocene-Quaternary, as well as the identification of magmatic intrusions. However, due to the lack of boreholes in the study area and adjoining sectors, a precise temporal framework for our seismostratigraphic reconstruction is not achievable. Some major seismic reflections can be assessed with confidence by correlation with CROP line M36, which passes above the volcanic structure E5c (after De Ritis et al., 2019) and crosses several MCS profiles presented in this study (Figure 2A).

In addition, a regional reflector associated with the top of the evaporites deposited during the late Messinian Salinity Crisis (MSC) is recognized by analogy with basins in the Western and Central Mediterranean, where it coincides with the top of the trilogy of units of the MSC (Dal Cin et al., 2016; Camerlenghi et al., 2020). The Lower Unit (LU) is inferred to be composed of gypsum and clastics, but is not visible in our data due to the poor energy of our source. The Mobile Unit (MU) is mostly a seismically transparent unit characterized by an upper folded boundary corresponding to evaporites filling the basins and onlapping into the lower continental slopes. The Upper Unit (UU) forms a distinct seismic facies, consisting of a package of high-amplitude parallel and continuous reflections, interpreted as an alternation of anhydrite or gypsum and marl layers. The erosion of these latter two units can be marked by the Margin Erosion Surface (MES), which is widespread on the continental slope of the Western Sardinian margin, where Messinian units are generally absent (Geletti et al., 2014). The MES can correspond with the regional M reflector recognized as the top of the Messinian units or regional unconformity in the Tyrrhenian Sea (e.g., De Ritis et al., 2019). This reflector is often located below the maximum penetration of our MCS data (e.g., Figure 3) and its depth surface has been mapped mostly with vintage sparker data (Figure 2B).

The overlaying Plio-Quaternary sedimentary succession is characterized by layered, mostly continuous, medium-to high-frequency and moderate to high-amplitude reflections. The Lower Pliocene is generally characterized by semi-transparent reflections (Figure 3).

In the northern study area a seismic facies characterized, internally, by high frequency and variable amplitude reflectors with sigmoid external geometry interpreted as “Lowstand Infralittoral Prograding Wedge” (LIPW) formed during pre-Last Glacial Maximum relative sea level lowstands has been identified by analogy with De Ritis et al. (2019) (Figures 3, 4A).

Volcanic seismic features exhibit a mound-shaped or an eye-shaped geometry, with a flat base and convex top, discontinuous and chaotic internal reflectors, which can be interpreted as laccoliths and magmatic intrusions reaching the sea floor. Volcanic chimneys and seepage of fluids along vertical faults are widespread in the shallower sediments (Figure 4A). Similarly to the laccoliths, they are confined in the northernmost part of the Paola Basin, indicating that the ascent of magma wanes moving away from the northern volcanic-bearing sector of the study area.

The interpretation of magmatic sills is a problematic task in our data. In particular, the acoustic imaging of sills and the distinction between repeated sills and multiples beneath the first high amplitude sill reflection event (Figure 3) is challenging. However, the presence of sills is envisaged at this site in the lower-resolution CROP M36 line, on which, the same feature is interpreted as a magmatic sill associated with a conduit (see Figure 4 in De Ritis et al., 2019). It is however, noteworthy to say that the seismic facies and striking folding of the reflectors corresponding with the magmatic sill are reminiscent of the acoustic response of the Messinian UU. In the Western Sardinia margin, the Messinian UU can be intercalated by (salt) lenses characterized by a more transparent facies (Geletti et al., 2014), with an overall seismic facies similar to what we observe in our case study (Figure 3). However, the presence of the Messinian UU in the upper and mid slope of the study area, where only marginal evaporites consisting of gypsum are believed to be present (Bertoni and Cartwright, 2015), is here considered unlikely. Additional high-resolution seismic data would be necessary to investigate this possible occurrence. Other sheet-like reflections characterized by high-amplitude that cross-cut the host-rock strata may be interpreted as sill remnants, severely disrupted or eroded by the observed diapirism of Messinian or pre-Messinian units (Figure 4B). Chaotic reflections and seismic attenuation stem from the sill and form a large vertical conduit below a 1.5 km wide eye-shaped feature characterized by an inner chaotic seismic facies with an external wall, which cuts surrounding clear reflections (Figure 4A). Bright spots are visible downslope from the conical edifice of Ov1, that overall can be interpreted as a volcanic feature (Figure 4C).

Proceeding farther south, in the non-volcanic sector of the study area, doming and diapirism of sediment below the MES reflector of the Messinian units remain significant (Figure 2B). In correspondence of the mound structure R1MV, the diapir gives rise to a prominent conduit which reaches 50 m below the seabed (Figure 5A), where bright spots and chimneys are connected to the presence of free gas in the shallow sediment (Figure 5C). The small chimney ascending to the seabed is barely visible also in the CHIRP profile (Figure 5D). Mass-transport deposits (MTDs) are located both at the seabed (Figure 5D) and 130 m below the seafloor along the western slope of the mound structure (Figure 5B) suggesting a repeating pattern of discharge. A chimney is visible ascending from the top of the MTD, suggesting that gas-charged sediment is present within the MTD, although the latter appears to be sealed by positive amplitude anomalies, indicative of harder material (Figure 5B). Plio-Quaternary reflections pinch out against the diapir walls but are disrupted by crater-like reflections followed upward by chimney-like reflections (Figure 5A). Locally, the diapirism and associated conduits can reach to the seabed where the high amplitude reflectors are directly connected to the seabed (Figure 6A), and where acoustic flares have been detected in the MBES water column data. Sill remnants are possibly present also at these locations (Figures 5A, 6B), while bright spots are quite evident (Figure 6C), also in the unfiltered shot gathers (Figure 6D).

Farther south, the depth of MES is always below the maximum penetration of our MCS data, however, the isobaths map highlights that the MES is shallower (Figure 2B) in correspondence of two mounded structures of very low relief (Figure 7A). Seismic attenuation is pervasive below scattered through-over-peak reflections that correspond to an impedance acoustic contrast, evidence of a harder layer located 150–160 m below the seabed (Figure 7B). Where these positive anomalies are more continuous, pockmarks are not present at the seafloor. The mound structures are not connected with gas flares in the water column, but are characterized by pockmarks and a crater on the seafloor (Figure 7B). Bright spots are present, but are buried below 25–50 m of sediment (Figure 7C). A saucer-shape reflection, that may resemble a magmatic sill, is present in the lower part of the Plio-Quaternary (Figure 7A).

In correspondence of the mound structure RMV, where the highest gas flare has been recorded, acoustic blanking is relevant (Figure 8A), while below the MMV structure a diapiric-like feature is imaged especially by RMS amplitude (Figure 8B). Acoustic blanking predominates where the acoustic impedance contrast likely indicates the presence of layers comprised of harder material (Figure 8B). Here, the observed scattered positive anomalies (Figure 8D) and pull-up events (Figure 8A) are probably connected with the methane-derived carbonates that were actually collected by sediment coring. Bright spots are located near the seafloor where the gas flares are present (Figures 8C,E); in the mound structure where no flares were detected (same of Figure 7A) the bright spots are just 10 m below the seafloor (Figure 8D). In correspondence of MMV, faults radiate at the summit of the diapiric conduit outlined by seismic attenuation (Figure 8C).

Due to increased impedance contrasts, resulting from enriched free gas content in the pore-space, gas in the subsurface becomes visible in CHIRP profiles as enhanced reflectors resulting in high amplitudes of the seabed reflector and few meters below and acoustic blanking underneath (Figure 5D). A map of these areas of enhanced high-amplitude in the shallow sediment succession (1–10 m below the seabed) is shown in Figures 9B, 11B. Mapping of high-amplitude reflections indicates that their occurrence is concentrated in the southern sector, around MMV and RMV, suggesting that these areas are influenced by higher gas concentrations. Subsurface acoustic blanking in CHIRP profiles is observed also as few-meters wide vertical chimneys reaching within 3–4 m from the seafloor (Figure 5D). Since they correspond with gas flares in the water column, or their projection at distances <500 m, and correlate with similar features in MCS data (Figure 5B), they are interpreted as gas pockets. Another type of high amplitude reflections, consisting of narrow vertical lineations, a few meters long, is frequently observed close to the seafloor and connects gas chimneys with flares in the water column (vertical bright spots) (Figures 9A, 11A). Vertical bright spots mostly corresponds with areas of pervasive landslide headwall scars (Figure 9A) in the northern sector of the study area, or pockmarks in the southern sector (Figure 11A).

Gas emissions were detected and identified as flares in water column echograms of the MBES. Due to the limited coverage of the swath needed for water column observations, the total area covered for water column flare imaging by swath coverage was restricted to 2,476 km² (Figure 2A). These data were loaded in the QPS Fledermaus Midwater tool using the stacked view as initial guide for possible venting, and the fan-view for detailed detection of gas emissions (Figures 9D, 11D). When a flare was identified, its acoustic signal was traced to the seafloor through subsequent fan-images and its geographic location was picked at the central point. The MBES data allow also defining flare locations when offset from the central ship track. However, in this case, identifying the outlet at the seafloor can be problematic due to background noise, especially when flares are not standing out clearly above the seafloor and weather conditions are not ideal. Depending on water depths, bubble size distribution, bubble rise rates, and possible turbulent flows, some frequencies may be better suited for detecting gas in the water column, and EM302 30 kHz acoustic data were the most efficient to detect the gas flares.

In total, 15 water column anomalies were recorded in the study area (Table 1). This represents a minimum value, since some distinct flares may appear clustered because their spacing is smaller than half of the MBES footprint used to detect the flares (Figure 9D). Flare observations were classified according to their appearance being certainly caused by gas bubble emissions. Relatively weak anomalies or with shapes, deviant from a continuous linear feature connected to the seafloor, were therefore discarded; however, their number was negligible compared to the columnar shapes. So far manual identification and extraction, with quality check on, is the common methodology used (Römer et al., 2021), although some attempts for automation have been achieved to distinguish gas flares from fish schools (Minelli et al., 2021). However, the study area is not renowned for fishing grounds, therefore the occurrence of fish schools was unlikely or very limited. As some areas were studied multiple times, flare observations were partly repeated and flare numbers have been corrected for double counting. Double detection during different survey times suggests that most flares are spatially and temporally stable, this is further constrained by the observation of the same flare, characterized by similar height, recorded at RMV in 2011 (Rovere et al., 2014) and 2014 (this study).

Flares were detected at seafloor depths of 548–862 m. Flares appear in heights from less than 20 m in the outlets surrounding the main flares and end mostly within the water column, the only exception being the two flares identified on top of RMV that reach very shallow water depths. Flare height determination is generally limited by the swath geometry, and the upper parts of the highest flares detected at RMV are probably cut off about 15–25 m below the sea surface (Figure 11D). However, no natural ebullition at the sea surface was observed during the surveys and there is no reason to suggest that gas exchange with the atmosphere is active at least in calm weather conditions and during the summer season, when both the 2011 and the 2014 surveys have been conducted (August and June).

The flares are always located directly above bright spots and stacked bright spots connected to the seabed or very close to it in the MCS profiles (Figures 6, 8).

Bathymetric mapping of the study area revealed that flare locations are systematically related to specific morphological seafloor elements (e.g., mounds, pockmarks, and linear scarps) (Figures 9A,C, 11A,C) or seafloor backscatter anomalies (methane-derived carbonate precipitation and fluidized mud at the seafloor) that are usually associated with gas seepage (Figures 9B, 11B). The flares align almost N−S along the seaward slope of the Paola Ridge (Figure 2A), the only exception being the flares in correspondence of the RMV. This mound structure is located more landward and is confined by a NE-dipping normal fault, which uplifts the mound with respects to the adjoining Propeller Basin (Figure 8A) and is particularly evident to the north (Figure 13).

The northernmost isolated flare rises above a pockmark located on the southern steepest flank (7–10° gradient slope) of Ov1, a 100-m-high almost conical mound which may be interpreted as an HTVC (Table 1, Figure 9B1). A sediment core collected on top of Ov1 recovered less than 1 meter of sediment, consisting of conglomerates and carbonate hardgrounds. On the western flank of Ov1, pockmarks are present but limited in number (Figure 9A), the largest being 100 m wide, and are probably related to fluid seepage through fractures from the buried MTDs located downslope (Figure 4A). The presence of MTDs characterized by bright spots and chimneys, indicative of the presence of free-gas in the sediment and seepage, may be related to past hydrothermal activity and gas discharge.

Other flares develop along a prominent scarp that connects the Ov1 HTVC with R1MV (Table 1), and forms a N−S oriented, 12-km-long feature consisting of a series of structural terraces originating a topographically complex slope sector (Figures 9A,C). The tallest flares (223–271 m) develop on a small depression (Figure 10), downslope from the scarp that makes up the western flank of R1MV (Figures 9C, B3). Here, flares are not obviously associated with pockmarks that are <5 m deep, and averagely 100 m in diameter (Figure 10). Flares develop above the smaller pockmarks, a depression and the slope (Figure 10).

Gas flares are connected with areas of highest seafloor backscatter and high-amplitude shallow reflections in CHIRP profiles (Figure 9B); however, along the northern sector of the study area, the high seafloor backscatter is also indicative of the presence of several landslide headwall scars and erosion at the seafloor that contribute to the character of the acoustic response. Therefore, in this case, the seafloor backscatter is not an exhaustive tool to distinguish areas characterized by fluid seepage and presence of gas on or close to the seabed. Furthermore, a further contribution to the high backscatter can come from the widespread shallow-seated MTDs, possibly related to the discharge of fluidized mud and mudflows from R1MV, as testified by CHIRP profiles (Figure 5D). This setting is similar to that of MTDs located downslope of the HTVC of Ov1, where amplitude anomalies are associated with fluidized sediment and gases (Figure 4C), and may correspond to eruptive phases, which have now ceased. In addition, the presence of carbonate hardgrounds, in general not associated with active venting (e.g., Svensen et al., 2003), points to the extinction of the main hydrothermal vent activity, on top of Ov1. However, there might be still residual thermal effect reflected by the release of gas through fractures and along the flanks of the edifice.

Moving southward, MMV and RMV sit on a large almost rounded (13.5 × 11.5 km) mound 300 m high. There are no significant morphological differences between MMV and RMV, in terms of height and size, however RMV appears to be fault-controlled. Above MMV, several clusters of gas flares (Table 1) are associated with evenly distributed pockmarks with diameter of about 50 m and average depth of a few meters (Figures 11A, B2). The western flank of MMV is characterized by a belt of landslide headwall scars, which might coincide with shallow-seated mudflows (Figure 11C). MMV and RMV are morphologically separated by the largest depression observed in the study area. The latter is a 250 m wide and 10 m deep, crater hosted in a subtle saddle (< 40 m deep) between MMV and RMV, not associated with any gas flares (Figures 11A, C, 12). On top of RMV the two highest gas flares detected in the study area (700 m against a water column of 736 m) develop above pockmarks and a depression, which has dimensions similar to the crater (Figure 12).

Also in the southern sector of the study area, the gas flares always coincide with high seafloor backscatter areas, which are however larger than those affected by gas discharge in the water column. Although slope gradient at the mounds flanks may enhance seafloor reflectivity, the very high backscatter of the crater may be a further evidence of the presence of gas in the first 25–50 m of the sub-seabed. This is further supported by bright spots in MCS data (Figure 7C), that may be related with a past outburst that formed the crater.

In other seep areas, steeper dipping and higher mounds are indicative of recent or prolonged growth of structures, whereas, lower-relief and gentler topography are associated with older features (e.g., Serié et al., 2012). In the study area, the differences in the morphology of the mounds are not so significant, and therefore may not be indicative of their different development stages and seepage activity. Furthermore, the seafloor reflector does not show significant variations in seismic amplitude with depth and topography, suggesting that the seafloor is everywhere covered by the same type of sediment, with the only exception coinciding with the erosional seafloor downslope R1MV.

Gas samples from two box-cores collected on RMV (flare B1) and R1MV (flare B3) revealed a chemical composition dominated by CO₂ (up to 98.73% by vol.) and subordinately by N₂ (up to 1.26% by vol.) and methane (< 0.06% by vol.). The carbon and oxygen isotopic ratios of CO₂ are −1.8 and −1.1 (V-PDB ‰) and −2.4 and −4.4 (V-PDB ‰), respectively. Data collected with the Niskin bottles carry negligible CH₄ and other hydrocarbon content, probably due to the difficulty in controlling the lowering of the rosette sampler inside the flares and reaching close to the bottom. The chemical composition of the Niskin samples is similar to present-day seawater, whereas the alkalinity is slightly higher (Franchi et al., 2017).

All flares detected are found in correspondence of subsurface of diapirs apparently originating from upward folding of the top Messinian units and diapirism connected to gas conduits directed to the seabed. The highest abundance of flares occurs around the diapir observed below R1MV (Figure 6A), while the highest gas flares were detected in connection with acoustic blanking connected with up warping of the Messinian units (not imaged by our MCS data) below RMV and doming and diapirism below MMV (Figure 8A).

Gas chimneys or conduits, that represent the path of vertically migrating gas into shallow reservoir sediments below bright spots, and that are located above salt domes has been detected in various geological contexts, such as: the German North Sea (Müller et al., 2018); the upper continental slope of the Santos Basin offshore Brazil (de Mahiques et al., 2017), where Miocene diapirs are associated with fields of pockmarks on the seabed; the Angola offshore (Serié et al., 2012), where a cluster of gas hydrate pingoes is attributed to a gas migration system along the flanks of a salt diapir. In the Mediterranean Basin, Messinian evaporites deposited in thick sedimentary successions in the deepest sub-basin depocenters as a consequence of the MSC. Although salt diapirs and domes are normally implicated in shaping the thermal and dynamic fluid circulation and seepage, this straightforward relationships is not always observed in the Mediterranean. In the Nile deep-sea fan, not all mud volcanoes are linked to the presence of thick Messinian evaporites (Bertoni and Cartwright, 2015). Even in the Levantine Basin, diapirism of Messinian salt and focused fluid flow reaching the seabed in the form of chimneys, which cross the thick basinal Messinian evaporites and originate in the pre-Messinian units, are not directly linked with mud volcanoes (Hübscher et al., 2009; Bowman, 2011). In the Western Mediterranean, the largest Messinian diapiric structures form in the Plio-Quaternary depocentres, where the sediment over-burden is larger (Geletti et al., 2014) and the buoyancy difference between indurated gypsum and overlying soft deep-sea muds favours the initiation of diapirism (Dal Cin et al., 2016). Less frequently, in these deep basins, salt diapirs cause intense deformation in the upper Plio-Quaternary, breach through the sea bottom and give rise to prominent surface fluid flow. On the contrary, in the slope context of the study area, where only gypsum deposits are present, which makes quite a different lithological and dissolution pattern compared to halite, diapiric structures concentrate where the thickness of the Plio-Quaternary is reduced in connection with up warping of the MES reflector (Figures 5A, 6A).

Furthermore, in the Western Mediterranean, salt deformation was intense during the Lower Pliocene, whereas in the Paola Basin, a more prolonged phase of up-rising and deformation of the diapiric structures, not necessarily connected with passive diapirism of evaporites, is indicated by pinch-out terminations of the upper Quaternary units (Figure 5A). This may be explained with a tectonic control over the emplacement of the diapirism in the Paola Basin, which prevail over sedimentary processes. On the other hand, the presence of massive free gas in the very shallow subsurface and vigorous and steady gas flares in the water column, as above R1MV and RMV, requires high fluid flux which is usually associated with a supply of thermogenic fluids by a deep-rooted plumbing system (e.g., Plaza-Faverola et al., 2015). In our case study, this plumbing system likely exploits normal faults, which are not obviously imaged in our MCS data because of the intense acoustic masking. This notwithstanding, a normal fault along the western boundary of the Propeller Basin is imaged (Figure 13). Above the conduits, radial faults are often present, especially on top of MMV (Figure 8C). The conduits below R1MV and MMV, that are associated with gas flares in the water column, align along the westward slope of the Paola Ridge which is tectonically controlled by normal faults (Gamberi and Rovere, 2010) and experienced subsidence during Pliocene to Quaternary times (Pepe et al., 2010).

In addition, polygonal faults are observed in the upper Plio-Quaternary to the east of the gas flares area of R1MV, they are related to unfocused fluid flow probably by the early expulsion of pore waters in the fine-grained sediments of the Plio-Quaternary. Only few of these faults have a subtle manifestation at the seafloor, indicating recent activity (Figure 5A). In the Western Mediterranean, polygonal faulting directly above the thick stratified evaporite sequence may have formed in response to water release from the underlying evaporites, due to gypsum-anhydrite conversion. The gypsum-anhydrite diagenetic reaction is temperature and depth-dependent (Bertoni and Cartwright, 2015). Therefore, given the relatively shallow-burial depth in the R1MV area, the likely presence of only gypsum in the Messinian units, a high geothermal gradient/heat flow should be invoked for addressing the diffuse fluid seepage along the polygonal faulting area east of R1MV.

The presence of magmatic sills in the Plio-Quaternary succession of the study area (Figures 4A, 5A) may represent the explanation for the additional thermal gradient needed for the observed intense fluid flow and gas venting in connection with up warping of thin Messinian evaporites in the study area. The upward folding of the Messinian units and seal bypass can be therefore associated with fracturing caused by the hot igneous intrusions, that are interconnected or lie just above the sealing succession (Figures 4A, 5A), as in other contexts such as the Rockall Basin in the NE Atlantic (Cartwright et al., 2007). The heating of the colder and wetter Plio-Quaternary sediments may have resulted in boiling of pore water, rapid maturation of organic matter in the sill aureole and overpressure around the sill with the formation of the fracture network along which gaseous hydrocarbons, water and magmatic material escaped to the seabed and formed eye-shaped or domed HTVCs (Figure 4A), where material is extruded into the seafloor and may resemble a mud volcano system (e.g., Iyer et al., 2013; Roelofse et al., 2021).

The assessment of the presence of sills in our MCS data remains problematic, although subtle evidence of magmatic sills may be detected in the lower part of the Plio-Quaternary (Figures 5A, 6B, 7A). The subtle evidence of magmatic sills may be due to: lack of seismic signal penetration; acoustic masking caused by the transparent facies of the diapirs and gas conduits (Figure 8A); fractures that may form for either thermal contraction during cooling, metamorphism in the contact aureole and rapid loss of hydrothermal fluids from the surrounding sediments (e.g., Cartwright et al., 2007; Figure 4A).

Seismic acoustic masking has been interpreted as gas chimneys or gas-charged fluid migration conduits also in areas not in connection with salt domes (Plaza-Faverola et al., 2015; Waage et al., 2019). A minor fraction of free gas in the pore space of sediments is sufficient to produce signal blanking in form of chimneys that can be thus regarded as gas-charged fluid migration pathways. On the other hand, acoustic masking and seismic attenuation can be due to the shadow effect by a sharp impedance contrast when short offsets are used in the acquisition of the seismic data, such as in the system used for this study. The acoustic masking zones underneath RMV and the other mounds in the southern sector of the study area crosscut the entire Plio-Quaternary succession (Figure 7A). The masking zones are located below scattered high-amplitude positive anomalies located 150–170 m below the seafloor, intercalated with bright spots and seismic pull-up events that may be associated with authigenic minerals (Figures 6A,B). This is confirmed by the fact that, where no high-impedance contrast is visible, the conduit rising to the seabed is well imaged by RMS amplitude, such as below MMV (Figure 8B). The scattered character of bright spots in the mounds, where no gas flares have been recorded (Figures 7C, 8D), may further suggest that the lower-impedance contrast is caused by gas within the pore space and fractures rather than gas bubbling at the seabed surface. On the contrary, the positive anomalies may be related to the presence of authigenic carbonates or other types of deeper mineralizations that hamper the seepage towards the seabed. The isotopic signature of the authigenic siderites (δ¹³C 3–10 V-PDB ‰, δ¹⁸O 8–9 V-PDB ‰) collected in the sub-seafloor of R1MV, RMV and MMV (< 10 m) indicate a precipitation in the methanogenic zone (Rovere et al., 2015). Therefore, the seep-carbonate horizons may be fed by the migration of thermogenic methane and hydrocarbons from deeper levels below the Messinian units, or by maturation of the organic matter favoured by the thermal anomaly connected to the magmatic sill as, for example, in the Guaymas Basin (e.g., Simoneit et al., 1988). In the area of the crater separating MMV and RMV, where gas flares are absent, a bowl-shaped deep reflection can be interpreted as a magmatic sill connected with the development of HTVCs that are not presently active. Here, fluids migrated through the lower part of the hydrothermal vent complex, which acted as a vertical zone of high-permeability, enhanced fracturing of the host sediment and circulation of fluids. However, the fractures that develop during the initial stages of hydrothermal activity within the metamorphic contact aureole are highly subject to mineralization and clogging by cementation (Cartwright et al., 2007). Therefore, it is here suggested that the majority of the HTVCs in the southern sector of the study area formed in consequence of an initial gas-dominated explosive phase, but soon became inactive due to their pervasive mineralization. The thermal effect continued to promote the maturation of organic matter and deposition of methane-derived siderites that in turn progressively acted as major barrier for fluid seepage to the seafloor (Figure 7A). Where gas venting have been observed in the water column, either the hydrothermal activity is still prevalent or precipitation of methane-derived carbonates is not sufficient to prevent prolonged gas migration.

Magmatic activity in sedimentary basins has a critical impact on fluid migration as hot intrusive rocks elevate the local geothermal gradient. The main temperature effect of magmatic sills is believed to last only few thousand years (Roelofse et al., 2021). However, normal faulting, which is quite common in hybrid sedimentary-volcanic systems, causes downward displacement of hanging wall sediments that are initially colder. This process results in the thermal instability of the basin, which may regain its steady state only after a long period, in the order of a few millions years (e.g., Sydnes et al., 2019). There are cases, such as the offshore of eastern Australia, where an Eocene phase of magmatism feeding hydrothermal vent complexes has been linked to the formation of buried mud volcanoes and is believed to still influence the occurrence of seep-related features at the seafloor (Rollet et al., 2012). Therefore, we interpret the Paola Basin as a site where magmatic intrusions, which remain challenging to image in our data but are likely present based of numerous linked observations, have provided the necessary thermal gradient to maintain an effective hydrothermal system.

In hybrid sedimentary-volcanic systems, inorganic CO₂ concentrations generally exceed 50 vol%, whereas CH₄ concentrations are generally lower (roughly > 1–2 vol%). Their surface manifestations, like muddy craters or bubbling pools, may be similar to, and thus may be confused with, pure sedimentary gas manifestations (hydrocarbon seeps), such as mud volcanoes (Procesi et al., 2019). According to a global review (Procesi et al., 2019), the δ¹³C-CO₂ carbon isotope signature, related to decarbonation reactions (thermo-metamorphism) and magma-mantle degassing at high temperatures (> 250°C) varies from −8‰ to +2‰. These values are in agreement with those recorded at the two gas venting sites sampled above R1MV and RMV (−1.8 and −1.1 V-PDB ‰), implying that there must be a contribution from deep magmatic and metamorphic systems to the gas released at the seabed in the Paola Basin, although a certain contamination from seawater CO₂ cannot be excluded in our isotopic signal, due to the sampling of gas bubbles from the box core. It must be added that the presence of an evaporitic sequence, although thin, such as the top Messinian unit, modifies pressure and temperature regime, adding a forcing to the sediment burial and fluid history of the basin, including complex feedbacks of diagenesis and dissolution. Furthermore, the diapirs originating from the Messinian units are connected through extensional faults to the deep-rooted plumbing systems. This setting enhances fluid migration of both abiotic CH₄, formed by deep magmatic and post-magmatic high-temperature reactions and biotic thermogenic CH₄ related to thermal reactions on organic matter, which is reflected in the precipitation of isotopically heavy carbonates.

The timing of the onset of fluid seepage in the study area must be better constrained. In the Mediterranean area, fluid flow systems formed around the Messinian time, and were sourced from the pre-MSC under saturated fluids and overpressured sediments, below the thick gypsum-halite sequences (Bertoni and Cartwright, 2015). The absence of paleo-pockmarks or limited presence of pockmarks buried by a thin (> 5°m) veneer of sediment (Figures 7B,D–G) and evidence of pinch-out terminations against the diapiric structures in the upper Plio-Quaternary sediment (Figure 5) point to a post-MSC activity of the seepage in the study area. In general, fluid flow dated after the MSC seems more dependent on the local geodynamic and depositional settings (Bertoni and Cartwright, 2015); therefore, a Pleistocene age would be more coherent with the establishment of magmatism in the area.

As a matter of fact, the geodynamic history of the Tyrrhenian Sea started at 10 Ma with the opening of the basin. The eastward migration and intensification of the rifting, between 10 and 5 Ma as a consequence of the slab roll-back, triggered seafloor spreading in the Vavilov back-arc basin and the formation of subduction-related volcanism. Between 5 and 0.6 Ma, the rifting migrated further to the southeast and seafloor spreading started in the Marsili back-arc basin, with the emplacemeant of oceanic crust (see Malinverno, 2012 and references therein). Around the same time (0.6 Ma), STEP (Subduction-Transform-Edge-Propagator) faults, associated with lateral lithospheric tearing, are believed to have formed (Cocchi et al., 2017; De Ritis et al., 2019) in the eastern Tyrrhenian domain. The STEP faults acted as the main conduit for upwelling of the isotherms due to mantle melts upraising from the edge of the slab along the lithospheric tear. Therefore, these inferred STEP faults would have controlled the magma uprising that fed the Palinuro volcanic complex (west of Glabro, outside of the area depicted in Figure 1B) and the Diamante-Enotrio-Ovidio volcanic-intrusive complex (Figure 1B). At present, the upwelling of subduction-induced mantle flow has stopped in the Diamante-Enotrio-Ovidio volcanic-intrusive complex, while the subduction-related volcanism continues in the Aeolian Arc (De Ritis et al., 2019). This reconstruction seems compatible with magmatism of Pleistocene age in the area located south of the Ovidio complex, which may have enhanced the fluid circulation in the Paola Basin, in a complex feedback system, including the top of Messinian units.

The hydrothermal vent complexes that develop along the eastern side of the Tyrrhenian domain deserve better understanding. There are indeed striking morphological and geochemical similarities with volcanic-hydrothermal systems located in shallow waters along the margin. For example, the morphology of seabed and nature of the gas is similar to active seabed doming and gas discharge in the Gulf of Naples (Figure 1A), where He and CO₂ are sourced from mantle melts and decarbonation reactions of crustal rocks (Passaro et al., 2016). Other large fluid-escape depressions and seafloor mounds have been observed on the continental shelf of the Pontine Archipelago (Figure 1A), where hydrothermal sulfides collected on the seabed and sub-seabed point to the possible degassing of magma similar to the one feeding the latest volcanic activity occurred on the islands in the Middle Pleistocene (Conte et al., 2020). Hydrothermal fluids discharged at this location are CO₂-rich and show a mantle-derived signature indicating that cooling magmas are still releasing enough thermal energy to feed an active hydrothermal system (Italiano et al., 2019). An example in water depths closer to our case study, is the fault-controlled system of deep-hydrothermal circulation discovered on the Cape Vaticano Ridge (Figure 1B). Here, δ³He enriched bottom waters are related to magmatic intrusions generated from mantle-wedge partial melting above the Tyrrhenian-Ionian subducting slab and melt upward migration (Loreto et al., 2019). Therefore, submarine hydrothermal systems and vent complexes along the eastern side of the Tyrrhenian domain remain overlooked especially in deeper waters, where geophysical and geochemical data are sparse. The study of surface expressions of fluid circulation, such as near-seabed seepage and fluid discharge in the water column, would enable to better understanding their relationships with the deep processes associated with slab-subduction, slab tearing and mantle-wedge partial melting.