Gas migration and slope instability in the Danube Fan: insights from integrated OBS-MCS seismic analysis

Gas hydrates and deltaic deep-sea fans are main features in continental margin systems, influencing slope stability, fluid migration, and carbon cycling. In the northwestern Black Sea, the Danube Fan remains poorly constrained with respect to subsurface structure, sediment strength, and hydrate dynamics. Here, we present high-resolution multichannel seismic (MCS) and ocean-bottom seismometer (OBS) data to characterise sedimentary structure and fluid-related features. Two integrated OBS–MCS profiles reveal under-consolidated, clay-rich levee deposits interspersed with mass-transport units, chaotic facies, and gas-related anomalies. Derived P- and S-wave velocity models indicate low shear strength and high $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratios in shallow units, consistent with soft, water-saturated sediments. Deeper layers display compaction-driven velocity increases but remain mechanically weak, rendering the slope prone to failure. Our findings suggest that vertical gas migration is widespread, expressed by seismic chimneys, polarity reversals, and velocity pull-downs, with free gas confined below bottom simulating reflectors and in stratigraphic traps. Hydrates likely occur as sparse, patchy pore-filling accumulations, and the lack of S-wave velocity anomalies suggests they do not act as cementing phases, implying little direct influence on sediment strength or slope stability. The hydrate system appears hydrate-poor, possibly reflecting post-glacial re-equilibration. Overall, lithology, gas migration pathways, and high sedimentation rates emerge as primary controls on hydrate formation and slope instability in the Danube Fan.

1 Introduction

Submarine slope stability and gas hydrate systems are critical components of marine geohazards and carbon cycling along continental margins. In regions such as the Danube Fan, located on the northwestern Black Sea margin, rapid sedimentation, overpressure development, and complex fluid migration contribute to slope instabilities and localised gas accumulations (Popescu et al., 2006, 2004). Understanding these processes is essential not only for assessing landslide risks but also for improving our knowledge of subsurface gas dynamics, which has important implications for marine resource management, carbon storage, and climate change feedbacks.

Gas hydrates, ice-like crystalline compounds of water and gas, typically methane, are stable under high-pressure and low-temperature conditions in marine sediments (Sloan and Koh, 2008). Their stability is sensitive to environmental changes, such as warming bottom waters or sea-level fluctuations, which may trigger dissociation and migration of free gas in the sediment and into the water column (Ruppel and Kessler, 2017). These dynamics are particularly relevant in continental margin settings like the Danube Fan, where gas migration pathways, overpressured sediments, and stratigraphic traps generate favourable conditions for both gas hydrate formation and slope failure (Riedel et al., 2021; Bialas et al., 2020a; Hillman et al., 2018; Zander et al., 2018).

Combined P- and S-wave analysis, provides a valuable tool for characterising sediment properties and further allow the estimation of $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratios that allow the identification of lithological changes, detect free gas and gas hydrates in the subsurface. Together with high-resolution multichannel seismic (MCS) imaging and ground-truth data, fluid dynamics can be studied intensively.

This study presents two high-resolution seismic profiles across a slump scar in the Danube Fan to investigate the subsurface architecture, sediment consolidation state, and gas hydrate occurrence. By integrating velocity modelling with stratigraphic interpretation and drilling data (MeBo200), we aim to constrain the mechanical and fluid-related processes governing slope failure.

Beyond addressing geohazards, this study has the potential to significantly advance our understanding of fluid migration and gas accumulation in under-consolidated marine sediments. Our insights are crucial for optimising site selection and enhancing monitoring strategies in subsurface projects involving fluid-injection, leakage control, and long-term storage, particularly in high-stakes applications, such as carbon capture and storage (CCS). The observed gas trapping and migration mechanisms are to a certain degree analogous to those required for long-term CO₂ storage in marine settings (e.g., Yip et al., 2013).

2 Geological Setting

The Danube Fan is a fine-grained turbidite system on a particularly wide shelf (∼120 km) of the northwestern Black Sea. It developed over the past ∼900 ka through a succession of stacked channel-levee systems interbedded with mass-transport deposits (MTDs) (Wong et al., 1997; Winguth et al., 2000; Popescu et al., 2001). The fan system extends from the shelf edge at ∼100 m water depth to the abyssal plain, where water depths exceed 2200 m (Wong et al., 1997).

Erosional processes dominate the upper slope, primarily within the channel-levee systems, while depositional and aggradational processes prevail along the middle and lower slopes (Popescu et al., 2001). Bathymetric data reveal the presence of several submarine canyons and associated channels, with a predominant NW-SE orientation (Fig. 1).

For consistency with previous studies, we adopt the canyon nomenclature established by Wong et al. (1997) and Winguth (1998), where the prefix “S” refers to the southern Danube Fan complex. Canyon numbers ascend from the oldest to the youngest, corresponding to successive activity phases and sedimentary sequences. These canyons represent the most recent depositional phase of the Danube Fan, which developed between ∼75 to 22 ka ago, during the sea-level lowstand associated with the Last Glacial Maximum (LGM; Late Pleistocene Neoeuxinian Stage) (Winguth et al., 2000; Lericolais et al., 2013). During this period, sea-level was ∼100–150 m lower than present, resulting in an isolation of the Black Sea from the global oceans (Ryan et al., 1997; Lericolais et al., 2009).

Figure 1 Bathymetric overview of the Danube deep-sea fan complex, acquired during MSM-34 (Wintersteller and Bialas, 2016), showing the location of the seismic datasets collected on the western flank of the S2 canyon. Red diamonds indicate sites where gas hydrates were sampled during GHASS (Gas Hydrates, fluid Activities and Sediment deformations in the western Black Sea) cruise (Ker and Riboulot, 2015). Inset: Location of the study area in the north-western Black Sea (ETOPO-1; Amante and Eakins, 2009). Country abbreviations: BG: Bulgaria, RO: Romania, MD: Moldavia, UA: Ukraine, RU: Russia, GE: Georgia, TR: Turkey.

Several studies suggest, that only one canyon was active during each sea level low-stand, with the Viteaz canyon (S5) representing the recently active system (Wong et al., 1997; Winguth et al., 2000; Popescu et al., 2001, 2004). During active phases, the canyons were directly connected to the Danube River and transported large volumes of sediment via channelised turbidity currents towards the deep sea (Popescu et al., 2001). Compared to the more distinctive S5 canyon, S2 exhibits lower sinuosity (Fig. 1). The surrounding seafloor contains extensive MTDs, levee deposits, and fine-grained overbank sediments, including channel overflows (Wong et al., 1997; Badhani, 2016). The isolation of the Black Sea during the LGM and the associated fan sedimentation ceased with postglacial sea-level rise and global reconnection at the end of the Quaternary glaciation, approximately 7500 years ago (Lericolais et al., 2009). Hydrological conditions in the Black Sea have fluctuated in response to glacio-eustatic sea-level changes, alternating between low salinity, oxygenated lacustrine environments and strongly stratified, anoxic marine conditions in deeper waters (Deuser, 1974; Özsoy and Ünlüata, 1997). The onset of deep-water anoxia is closely linked to the formation of sapropels, which result from enhanced biogenic productivity during transitions from limnic to marine settings (Popescu et al., 2007). Today the Black Sea is the world's largest anoxic basin, providing favourable conditions for hydrocarbon formation (Vassilev and Dimitrov, 2002; Jørgensen et al., 2004).

The estimated upper limit of the gas hydrate stability zone (GHSZ), assuming an average bottom-water temperature of 9 °C (Ross and Degens, 1974), lies at ∼720 m water depth (Zander et al., 2017; Riedel et al., 2020), though this depth also depends on the pore water salinity. Pore water salinity in the Black Sea ranges from 22.3 practical salinity units (psu) under marine conditions (Özsoy and Ünlüata, 1997) to as low as 3 psu in limnic conditions (Svitoch et al., 2000).

Within the study area, gas emissions (gas flares) are observed at water depths shallower than 655 m, the upper limit of the GHSZ (Hillman et al., 2018). Gas samples from the Danube Fan show a primarily biogenic origin, with concentrations up to 99 % methane, indicating in-situ formation (Vassilev and Dimitrov, 2002).

3 Data Acquisition, Processing and Modelling

As a part of the multidisciplinary research project SUGAR (SUbmarine GAs Hydrate Reservoirs) and MIDAS (Managing Impacts of Deep-SeA ReSource Exploitation), the German research cruise MSM-34 acquired 2-D multichannel seismic (MCS) and 4-component ocean-bottom-seismometer (OBS) data (Dannowski and Bialas, 2025), along with multibeam bathymetry (Bialas et al., 2014). Bathymetric data were collected continuously throughout the cruise using a hull-mounted EM 122 echosounder operating at 12 kHz. The data were processed onboard using generic multibeam processing software (Caress et al., 2008), with sound velocity corrections from regular profiles, manual artefact removal and cleaning, resulting in a gridded dataset with a spatial resolution of 25 m × 25 m.

Seismic data acquisition employed a 0.74 L $/$ 0.74 L Generator-Injector (GI) airgun, towed at 2 m depth and fired at 3 s intervals. At an average vessel speed of ∼3.5 kn, this configuration yielded a shot spacing of ∼5.4 m. The GI airgun had a centre frequency of ∼130 Hz.

3.1 OBS Data and Processing

Two OBS profiles, designated A and B, were acquired on the western flank of the S2 canyon within the Danube Fan. Each line included four OBS deployments (Fig. 2). The instruments were equipped with 3-component 4.5 Hz geophones and hydrophones, and data were recorded continuously using a Marine Broadband Seismic recorder (Bialas and Flueh, 1999) at a sampling rate of 1000 Hz.

OBS data processing included the following steps: (1) clock drift correction, (2) instrument relocation on the seafloor using direct wave arrivals, (3) rotation of the horizontal components into radial and transverse components for each shot (Maercklin, 2001), (4) amplitude balancing (max/min normalisation), and (5) bandpass or low-pass filtering to enhance phase identification during traveltime picking. A detailed description of the wavefield rotation procedure, including vector fidelity considerations and the basis for horizontal component rotation, is provided in Sect. S1 in the Supplement.

Overall, the OBS data quality was high. While simple plane-wave theory (e.g., Stein and Wysession, 2003; Yilmaz, 2001) predicts a smooth decay of amplitudes toward zero at normal incidence, the data shown in Fig. 3 deviate from this idealised behaviour. This is because the OBS receivers are located at water depths of ∼600–800 m, placing them in a transition zone between the near-field and far-field of the seismic source. In this regime, the recorded wavefield is strongly non-planar and influenced by near-field contributions, finite-source effects, and instrument-seafloor coupling, including potential tilt effects (e.g., Duennebier and Sutton, 1995; Essing et al., 2021). Non-planar and multimodal wavefield behaviour has been shown in numerical wave-equation simulations to deviate substantially from predictions based on simple plane-wave theory assumptions (e.g., Bai et al., 2014).

Figure 2 Detailed bathymetric map of the study area on the western flank of the S2 canyon, gridded at 25 m × 25 m resolution. Data gaps were filled with the global ETOPO-1 dataset (1 arcmin resolution; Amante and Eakins, 2009). The map shows the spatial extent of the velocity models and associated 2-D MCS profiles (red lines), ocean-bottom-seismometer (OBS) positions (blue dots), gas flares observed in the water column (white stars; Hillman et al., 2018), sampled gas hydrates (red diamonds; Ker and Riboulot, 2015), and MeBo200 drilling sites (yellow rectangles; Riedel et al., 2020). The dashed black line indicates the 720 m isobath, representing the estimated upper boundary of the GHSZ (Zander et al., 2017; Riedel et al., 2020).

Figure 3 (a) Hydrophone component of OBS 3005 and (b) identified P-wave phases used for forward modelling. Inset: Enlarged view of near-seafloor P-wave arrivals. (c) Radial component of OBS 3005 and (d) corresponding converted S-wave phases used in forward modelling. (e) Hydrophone component of OBS 3003 and (f) identified P-wave phases. Note that the horizontally stretched onsets around x=0 km result from a lateral offline position of OBS 3003. Inset: Enlarged view of near-seafloor P-wave arrivals. (g) Radial component of OBS 3003 and (h) corresponding converted S-wave phases used for forward modelling. Note that the S-wave panels are shown on a different timescale due to slower wave propagation velocities.

Such behaviour has been documented in previous OBS studies, which highlight significant deviations from plane-wave assumptions in complex geological settings (e.g., Bünz et al., 2005; Exley et al., 2010; Wang et al., 2018). Additional amplification of horizontal motion may arise from seafloor effects such as sediment compliance, impedance contrasts, ray-angle projections, and shallow conversions (e.g., Shearer and Orcutt, 1987), all of which contribute to the observed signal characteristics at normal incidence.

Multiple P-wave (hydrophone) and converted S-wave (geophone) reflections were clearly identified (Fig. 3). In general, S-wave reflections exhibited a higher signal-to-noise ratio than P-wave signals and were less affected by water-column multiples due to the different frequencies analysed (Fig. 4). To suppress high-frequency noise and emphasise the stable, low-frequency signal content, a bandpass filter of 20–40 Hz for P-waves and a lowpass of 20 Hz for S-waves was applied. Several stations drifted when they were deployed and are now about 30 m offline.

Figure 4 OBS 3005 hydrophone dataset (a) Alignment of the hydrophone with the 2-D high-resolution MCS data after applying a Butterworth bandpass filter (20–40 Hz) and equal amplitude scaling. The hydrophone data are adapted to the MCS section in terms of onset and offset timing. The MCS time axis is converted into one-way traveltime. (b) Bottom-simulating reflector (BSR) identified in the first arrival (left) and corresponding first multiple (right). The BSR is observed as a distinct reflection with reversed polarity relative to the seafloor reflection, consistent with a P-wave velocity decrease below the reflector (Fig. 6), likely due to free-gas accumulations. Such free-gas accumulations reverse polarity by reducing acoustic impedance relative to the overlying sediments.

One station (OBS 3002) was excluded from S-wave analysis due to insufficient seafloor coupling. The absence of clear refracted phases across the array restricts the velocity analysis to a localised zone directly underneath each instrument.

3.2 2-D High-resolution MCS Data and Processing

Multichannel seismic (MCS) data were acquired simultaneously with the OBS deployments along profiles A and B (Fig. 2), using a 144-channel streamer. The streamer configuration included a 25 m vibration isolation section, a 55 m tow cable, and 19 active sections of 12.5 m each, yielding a total active length of 237.5 m. Each active section contained 8 channels, with a hydrophone group spacing of 1.5625 m. Data were recorded at a sampling rate of 1000 Hz and binned with a common depth point (CDP) spacing of 1.5625 m.

The MCS data processing workflow comprised the following steps: (1) navigation correction for both streamer and shot positions, (2) crooked-line common mid-point (CMP) binning, (3) gain correction, (4) attenuation of anomalous amplitudes, (5) radon filtering, (6) CMP stacking using water velocity, necessitated by the short streamer length, (7) deconvolution with a prediction lag of 60 ms, averaged over 1000 traces, to suppress injector bubble effects, (8) water-velocity Stolt migration, followed by residual finite-difference migration using spatially and temporally varying velocities from regional profiles (Hillman et al., 2018). Final depth conversion was performed in generic seismic interpretation software, utilising the velocity models derived from the OBS data.

3.3 Phase Identification and OBS-MCS Correlation

The OBS records show both compressional (P-) and shear (S-) wave phases. Phase identification began with picking the P-wave reflections, which were then assigned to specific layers within the seismic velocity model for each station. Prominent reflectors observed in the MCS data served as reference horizons and were cross-checked against corresponding phases in the OBS records. Therefore, the OBS data were aligned to the MCS at zero offset and time-shifted so that the direct water arrivals matched the seafloor reflection in the MCS data, compensating for the two-way traveltime (Fig. 4a). This alignment allowed consistent phase identification across all OBS stations. Additionally, first water multiples were used to help discriminate overlapping arrivals to identify reflection polarity more clearly, especially where direct arrivals were masked by source-related noise (Fig. 4b). P-wave reflections were typically traceable to source-receiver offsets of up to ∼1.2 km (Fig. 3).

Once P-wave arrivals were identified, converted S-wave phases, recorded on the radial components (Fig. 3c and g), were assigned to the same reflectors within the velocity model. This step was the most critical in the analysis, relying on joint forward modelling of both P- and S-wave arrivals. A detailed description of the identification strategy for converted S-wave phases is provided in Sect. S2. Following Digranes et al. (1998), the conversion from P- to S-waves is assumed to occur at the reflector itself. Due to their lower velocities, S-wave reflections could only be traced over shorter offsets, with a maximum source-receiver distance of ∼0.9 km (Fig. 3), resulting in a smaller spatial coverage compared to the P-wave data.

3.4 Forward Modelling

Seismic phases were picked using PASTEUP (Fujie et al., 2008), a graphical user interface (GUI) for Seismic Unix (Stockwell, 1999). Forward modelling was performed with RAYINVR (Zelt and Smith, 1992) using the companion GUI MODELING (Fujie et al., 2008). The initial P-wave velocity structure was adapted from Hillman et al. (2018) and iteratively refined using a trial-and-error forward modelling approach. Interface geometries were locally adjusted by up to ±5 m over ∼200 m of profile length in the upper sections and by ±10 m over ∼500 m in deeper sections, particularly where OBS instruments were located offline.

RAYINVR computes traveltimes through user-defined velocity and depth grids using a 2-D ray-tracing and forward modelling algorithm (Zelt and Smith, 1992). This process yielded detailed final models extending to ∼1.2 km depth, well below the BSR.

Subsequently, the velocity structure and layer geometry were jointly refined for both P- and S-waves by modifying the $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratios within each layer. The plausibility of the P-S-phase correlations was evaluated by examining the S-wave residuals for systematic misfits or unrealistic $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratio values.

In addition, the S-wave parametrisation was intentionally kept less flexible compared to the P-wave model. Although the modelling software allows the definition of vertically varying $V_{\mathrm{p}}/V_{\mathrm{s}}$ values (top and bottom layer), this option was not used. Instead, constant $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratios were applied across layers. This choice was made to maintain model stability and interpretability within the constraints of the data.

Given the limited offset range and reduced ray coverage of converted S-wave phases, the sensitivity to vertical velocity gradient is low. Introducing additional degrees of freedom would increase the non-uniqueness in the S-wave model. The resulting V_s structure, therefore, reflects a balance between data resolution and model complexity.

4 Results

The results presented here integrate OBS-derived P- and S-wave velocity models with coincident MCS reflection data to characterise the subsurface of the S2 slope failure zone. Together, these datasets provide complementary constraints on sedimentary architectural, physical properties and fluid-related processes across profiles A and B. Both profiles reveal a consistent three-unit stratigraphic framework, while also highlighting local variations associated with mass-transport deposits, deformation structures, gas-related seismic anomalies, and potential fluid migration pathways.

4.1 Subsurface Structure of Profile A

Profile A traverses the central part of a slope failure zone west of the S2 canyon. OBSs 3001 to 3003 are located within a slump scar, with an average dip angle of 2.5°, while OBS 3004 lies toward the southeastern margin of the study area (Fig. 2). The northwestern headwall forms a ∼100 m high step in the seafloor (Fig. 5a). The profile ends within levee deposits of the Danube Fan.

Based on traveltime picking and forward modelling, three main units (Unit 1–3; Fig. 5) are identified, comprising ten model layers (A1–A10; details in Sect. S3). A1 represents the water column.

Unit 1. Unit 1 shows continuous, high-amplitude reflections with a thickness of ∼40–80 m and can be traced upslope within the slump scar (1.0–5.5 km along profile) and southeast of OBS 3004 (Figs. 5 and 6). Amplitude strength decreases slightly toward the southeast, while stratification remains continuous. A lenticular reflection body (M1), ∼ 20 m thick and ∼500 m wide, is observed between 1.6–2.0 km model offset (Figs. 5 and 8). This body displays semi-transparent, sub-parallel internal reflections.

Unit 2. Unit 2 is characterised by chaotic seismic facies (CSF) with patchy high amplitudes (Figs. 5 and 8). It extends from ∼0.9 to ∼3.6 km along the profile and is ∼150 m thick. Its base is defined by a prominent downward-dipping reflection (DDR; Figs. 5 and 8). Locally, southeast-tilted reflections are observed down to ∼800 m depth. Reflections between 3.5 and 3.8 km model offset dip toward the S2 canyon. Localised high-amplitude reflections occur above a zone of discontinuous and disturbed reflections at ∼4.4 km model offset at 800–1000 m depth. Below, a velocity pull-down anomaly is observed (labelled SC in Fig. 8).

Unit 3. Unit 3 consists of stratified, laterally continuous reflections with a slight northwestward dip, extending to ∼1.2 km depth. Reflections amplitudes decrease with depth and towards the northwest (Figs. 5 and 8), partly due to processing artefacts. Local polarity reversals are observed near a zone of velocity pull-downs (SC). Below ∼1.1 km depth, reflections remain continuous.

4.2 P- and S-wave Velocities of Profile A

The three seismic units identified in the reflection data on profile A are further characterised using the derived P- and S-wave velocity models. The derived V_p and V_s models for profile A are displayed in Fig. 5a and b and summarised in Table S1 in Sect. S3. V_p increases from 1510 m s⁻¹ at the seafloor to 1900 m s⁻¹ at ∼1.2 km depth, while V_s increases from 140 to 860 m s⁻¹ over the same interval (Figs. 5 and 6). The water column is assigned a constant velocity of 1483 m s⁻¹. Both V_p and V_s show minimal lateral variation, resulting in laterally stable $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratios within each layer.

Within Unit 1, V_p increases from 1510 m s⁻¹ beneath the seafloor to 1530 m s⁻¹ in ∼720 m depth, and V_s from 140 to 230 m s⁻¹, resulting in a $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratio decrease from 10.6 to 6.8. Within Unit 2, V_p increases from 1580 to 1760 m s⁻¹ at ∼880 m depth, and V_s from 320 to 460 m s⁻¹, with $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratios decreasing from 5 to 3.8. In Unit 3, V_p increases from 1760 to 1900 m s⁻¹ in ∼1.2 km depth, and V_s from 500 to 860 m s⁻¹, yielding in a $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratio decrease from 3.6 to 2.2. A local velocity decrease (∼80 m s⁻¹ in V_p) is observed at the base of Unit 2 (DDR) associated with a normal-polarity amplitude reflection (Fig. 6).

4.3 Subsurface Structure of Profile B

Profile B crosses a steeper terrain than profile A and intersects the slope failure zone near the NNE sidewall and a gully system entering the S2 canyon (Fig. 2). Based on traveltime picking and final modelling, the subsurface is divided again into three units (Unit 1–3), comprising eleven model layers (B1–B11), with B1 representing the water column (details in Sect. S3). This structure is consistent with profile A.

Unit 1. Similar to profile A, Unit 1 is characterised by continuous high-amplitude reflections, but with a stronger variation in thickness, ranging from ∼100 m northwest of the slump scar (0–1.8 km model offset) to ∼30–50 m within (Figs. 7 and 9). Two normal faults (F2 and F3) locally disrupt reflectivity between 1.8 and 2.7 km model offset and extend to ∼100 m below the seafloor.

Unit 2. As in profile A, Unit 2 consists of a chaotic seismic facies (CSF) confined between the faults F2 and F3 (∼1.8–2.7 km model offset) and extends down to ∼900 m depths (Figs. 7 and 9). The CSF is characterised by discontinuous and spatially variable reflection patterns.

Figure 5 Seismic velocity structure models of profile A. (a) P-wave velocity, (b) S-wave velocity, and the corresponding (c) $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratio models overlaying the 2-D MCS seismic data. The seismic section is vertically exaggerated by a factor of three. Inverted red triangles indicate the position of the OBSs along the profile. Seismic anomalies are pointed out: BSR: bottom-simulating reflector, CSF: chaotic seismic facies, DDR: downward dipping reflector, F1: fault, M1: mass transport deposits (MTD), and SC: seismic chimney structure. The dashed white lines mark the boundaries of Unit 2. The dotted red line marks the theoretical base of the present-day GHSZ, calculated using a bottom- water temperature of 9 °C and a geothermal gradient of 22 °C km⁻¹. The blue dotted line indicates the estimated GHSZ base during the Last Glacial Maximum (LGM), assuming a bottom temperature of 4 °C and the same thermal gradient (both after Badhani, 2016). Note: The theoretical BSR marks the upper limit of the GHSZ. At this shallow depth, conditions are not fully stable, so the BSR does not precisely follow the seafloor.

Figure 6 1-D velocity and ratio profiles for OBS stations along profiles A and B. (a) 1-D P-wave (solid lines) and S-wave (dashed lines) velocity profiles for OBSs 3001–3004 along profile A. (b) Corresponding 1-D $V_{\mathrm{p}}/V_{\mathrm{s}}$ and Poisson's ratios. Note: Poisson's ratio is non-linearly related to the $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratio. (c) 1-D P-wave (solid lines) and S-wave (dashed lines) velocity profiles for OBSs 3005–3008 along profile B. (d) Corresponding 1-D $V_{\mathrm{p}}/V_{\mathrm{s}}$ and Poisson's ratios. The horizontally dashed lines represent the theoretical base of the GHSZ (after Badhani, 2016) and the red triangle marks the observed BSR depth at OBS 3005. The different depth of the theoretical GHSZ are a direct response to the prevailing pressure and temperature conditions in the shallow setting.

Southeast, the CSF abuts a lens-shaped seismic body (M2), displaying semi-transparent internal reflections and local cross-cutting reflectors, some of which may represent processing artefacts. Adjacent reflections are locally deformed around M2.

Further southeast, two additional semi-transparent seismic bodies (M3 and M4) occur within otherwise sub-horizontal strata. M3 forms an elongated, lens-shaped body extending from ∼4 to 5.2 km model offset at ∼850 m depth, with an internal topographic step of ∼20–30 m. A laterally continuous reflector separates M3 from the underlying body M4 (Figs. 7 and 9).

M4 extends from 3.7–5.4 km model offset at depths of ∼850–950 m and thickens from ∼50 to ∼70 m toward the southeast. Internally, M4 shows deformed reflections and localised zones of elevated amplitude, while reflectivity gradually decreases toward the northwest.

Southeast of 5.4 km model offset, several V-shaped reflectors extend from the seafloor to depths of ∼1 km. Owing to their lack of lateral continuity and stratigraphic coherence, these features are interpreted as likely processing artefacts.

Unit 3. Unit 3 exhibits the same general character as observed in profile A, consisting of sub-horizontal and laterally continuous reflections (Figs. 7 and 9). A major difference is the presence of a strong upward-bending reflector (BSR1) forming the upper boundary of this unit. BSR1 crosscuts the surrounding strata, mimics the seafloor morphology, and exhibits reversed polarity. Below BSR1, several zones of enhanced-amplitude, reversed-polarity reflections are observed. A second, weaker, reversed-polarity reflector (BSR2) occurs at ∼1050 m depth (Figs. 7 and 9). Beneath both BSR horizons, reflections remain laterally continuous and sub-horizontal, with amplitudes gradually decreasing toward the northwest.

4.4 P- and S-wave Velocities of Profile B

The three seismic units identified on profile B are further characterised using the derived P- and S-wave velocity models. The derived V_p and V_s models for profile B are displayed in Fig. 7a and b and summarised in Table S2 in Sect. S3. V_p increases from 1510 m s⁻¹ at the seafloor to 1900 m s⁻¹ at ∼1.15 km depth, and V_s increases from 140 to 810 m s⁻¹. The water column (B1) is assigned a constant velocity of 1483 m s⁻¹.

Overall, the velocity structure of profile B is consistent with profile A, showing laterally stable V_p and V_s distributions and comparable vertical gradients within the three-unit framework. No major lateral velocity heterogeneities are observed that would require a deviation from the established stratigraphic subdivision.

Within Unit 1, V_p increases from 1510 m s⁻¹ at the seafloor to 1530 m s⁻¹ at 50–100 m depth, while V_s increases from 140 to 230 m s⁻¹, corresponding to a $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratios decrease from 10.6 to 6.8.

Within Unit 2, V_p increases from 1600 to 1780 m s⁻¹ at ∼1 km depth, and V_s increases from 280 to 480 m s⁻¹, resulting in a $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratio decrease from 5.8 to 3.7. This decrease is consistent with profile A, but with slightly higher absolute V_s values.

In Unit 3, V_p increases from 1800 to 1900 m s⁻¹ at ∼ 1.2 km depth, while V_s increases from 490 to 810 m s⁻¹, leading to a $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratios decrease from 3.7 to 2.4.

A key difference compared to profile A is observed at the position of BSR1, where a local decrease in V_p of ∼ 190 m s⁻¹ occurs. V_s values remain largely unchanged across this interval, resulting in a distinct reduction in the $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratio. At BSR2, velocity values are not constrained by ray coverage, and the reflector is therefore not resolved in the model.

Figure 7 Seismic velocity structure models of profile B. (a) P-wave velocity, (b) S-wave velocity, and the corresponding (c) $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratio models overlaying the 2-D MCS seismic data. The seismic section is vertically exaggerated by a factor of three. Inverted red triangles indicate the position of the OBSs along the profile. Seismic anomalies are pointed out: BSR: bottom-simulating-reflector, CSF: chaotic seismic facies, F2 and F3: fault, and M2 to M4: MTD (mass-transport deposit). The dashed white lines mark the boundaries of Unit 2. The dotted red line marks the theoretical base of the present-day GHSZ, calculated using a bottom-water temperature of 9 °C and a geothermal gradient of 22 °C km⁻¹. The blue dotted line indicates the estimated GHSZ base during the Last Glacial Maximum (LGM), assuming a bottom temperature of 4 °C and the same thermal gradient (both after Badhani, 2016). Where the upper limit of the GHSZ outcrops, BSRs may cross underlying units and do not perfectly follow the seafloor, reflecting local pressure and temperature conditions (e.g., Zander et al., 2017, 2018). (d) Simplified stratigraphic framework at the projected MeBo-17/19 borehole location (modified after Riedel et al., 2020) showing three main units: Unit A (silty clays with sand layers) intercalated MTDs representing chaotically reworked material, and Unit B (older fine-grained sediments overchannel-overspill turbidites).

4.5 Constraints and Model Uncertainties

Phase picking and forward modelling are inherently non-unique problems, as the solution depends on manual phase identification, model parametrisation, and the iterative fitting procedure (Diehl, 2008). To quantify and constrain model reliability, we distinguish between traveltime misfit, parameter uncertainty, and interpretational ambiguity.

Traveltime misfit is assessed through picking uncertainties, ray coverage (Figs. S10–S13 in Sect. S3), and statistical residuals (Zelt, 1999). The frequency-dependent effects of filtering (0–20 Hz for S-waves, and 20–40 Hz for P-waves) suppress high-frequency noise and emphasise the stable, low-frequency signal content that improves the reliability of traveltime picking. Pick uncertainties range from 5 to 24 ms (Tables S1 and S2 in Sect. S3) and increase with depth due to decreasing signal-to-noise ratio and reduced phase coherence. In total 8192 and 10 063 picks were used for analysis for profiles A and B, respectively. The resulting Root-Mean-Square (RMS) residuals (3 and 4 ms) and normalized χ²-values (0.27 and 0.33), for profile A and B, respectively, indicate an overall consistent fit between observed and synthetic arrivals. These metrics, however, describe data fit only and do not imply solution uniqueness.

Parameter uncertainty is primarily controlled by ray coverage (Figs. S10–S13 in Sect. S3), discretisation of the velocity model, and the trade-off between velocity and interface geometry inherent to forward modelling (Zelt and Smith, 1992). Ray density reflects both resolution and smoothness and is generally high across the profiles. Single-parameter sensitivity tests indicate allowable perturbations of approximately 1 %–7 % (10–100 m s⁻¹) in V_p, 2 %–14 % (8–30 m s⁻¹) in V_s, and 1 %–9 % (±0.031–0.707) in $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratios within the traveltime uncertainty bounds. A reduced perturbation can be reached by jointly constraining the model for P- and S-waves, resulting in 1 %–3 % (10–50 m s⁻¹). Perturbations were applied layer-wise by changing the top and bottom velocities and keeping the gradients constant; effects on adjacent layers were not included. Vertical resolution is the highest in the upper ∼100 m (10 m node spacing) and decreases with depth (up to 100 m node spacing), while lateral resolution is limited by velocity node spacing (100–500 m). Therefore, depth nodes have a sensitivity of 1–2 m for perturbation. Consequently, small-scale heterogeneities below these scales are not resolvable and are effectively absorbed into smooth velocity gradients or interface adjustments.

Interpretational uncertainty arises from the non-uniqueness of structural assignment. This is mitigated by independent a priori constraints from MCS stratigraphy, previously published velocity models of the Danube Fan (Dannowski et al., 2016; Popescu et al., 2006; Winguth, 1998; Zander et al., 2017), and MeBo200 drilling results (Riedel et al., 2020), which define robust bounds on lithological layering and sediment properties. These constraints reduce ambiguity in correlating velocity contrasts with geological units but do not eliminate equivalent model solutions that satisfy the traveltime data.

Differences in the resolution of P- and S-wave models reflect acquisition geometry and phase physics rather than methodological inconsistency. P-wave modelling benefits from higher signal-to-noise ratios and a broader offset range (∼1.2 km), resulting in denser ray coverage and smoother velocity fields. In contrast, S-wave constraints are limited to P-to-S converted phases with restricted aperture (∼0.9 km), reducing resolution and increasing spatial smoothing. Accordingly, a single $V_{\mathrm{p}}/V_{\mathrm{s}}$ value per layer is used to avoid overparameterization unsupported by the data. Since V_s values depend on V_p and $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratios, they were derived using error propagation.

Overall, the presented models represent a geologically constrained, data-consistent solution within a bounded ensemble of admissible models rather than a unique subsurface reconstruction.

5 Discussion

This study integrates OBS and MCS profiles from the central (profile A) and eastern (profile B) parts of the S2 slope failure zone. The data are interpreted in the context of velocity observations from the Danube Fan (e.g., Popescu et al., 2006; Zander et al., 2017; Bialas et al., 2020a) and MeBo200 drilling results (Riedel et al., 2020; Bohrmann et al., 2018). The combined dataset reveals a range of structural, lithological and fluid-related features including under-consolidated sediments, indications of focused fluid migration, gas-related seismic anomalies, and MTDs.

5.1 Sedimentary Architecture

Seismic velocities and the integrated OBS-MCS data reveal systematically low P-wave velocities and elevated $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratios in the shallow subsurface (Figs. 5 and 7, Tables S1 and S2 in Sect. S3), reflecting water-rich clays and silts with poor consolidation (Sava and Hardage, 2009; Hamilton, 1979a). Measured P-wave velocities remain below typical reference values for normally consolidated sediments at comparable burial depths (Hamilton, 1979a), indicating a persistently under-consolidated mechanical state throughout the upper sedimentary units. This suggests limited dewatering and incomplete compaction within rapidly accumulated fine-grained deposits.

Drilling results from the MeBo200 campaign (Bohrmann et al., 2018; Riedel et al., 2020, 2021) provide lithological ground truth for the seismic interpretation. The drilling campaign in the Danube deep-sea fan (cruise M142; Bohrmann et al., 2018; Riedel et al., 2020) was conducted to investigate gas hydrate occurrence and constrain sedimentary structures using in-situ coring and logging. The drilling activities were carried out in the Romanian sector of the Black Sea, where multiple sites along the S2 Channel and adjacent canyon flanks were targeted based on pre-existing seismic data. Boreholes reached depths of up to ∼144–147 m below seafloor, recovering Holocene to glacial sedimentary sequences that provide key lithological and petrophysical constraints for seismic interpretation.

Figure 8 Seismic section (upper panel) and interpretation (lower panel) of high-resolution MCS profile A. Gas-charged zones are shaded in light blue. Small-scale mass-transport deposits (MTDs), including M1 and chaotic seismic facies (CSF1), bottom-simulating reflector (BSR), as well as the lithological boundary DDR (downward-dipping reflector), are highlighted in yellow. A seismic chimney (SC) is visible southeast of OBS 3004. The dotted red line marks the theoretical base of the present-day GHSZ, calculated using a bottom-water temperature of 9 °C and a geothermal gradient of 22 °C km⁻¹. The orange dotted line indicates the estimated GHSZ base during the Last Glacial Maximum (LGM), assuming a bottom temperature of 4 °C and the same thermal gradient (both after Badhani, 2016). Boxes indicate enlarged sections shown in Fig. 10.

The MeBo-17/19 sites are located in proximity to OBS 3005 along profile B (Fig. 2), when projected onto profile B, with an approximate lateral offset of ∼125 m northwest of the OBS and ∼250 m southwest of it. Given these relatively small distances compared to the seismic resolution of the MCS and the laterally continuous character of the stratigraphy, the borehole data can be considered broadly representative of the seismic units, while small-scale heterogeneity remains unresolved. The drilled sections are dominated by hemipelagic clays and silts interbedded with sandy levee deposits, consistent with channel-levee and overbank environments of the Danube Fan system (Riedel et al., 2020, 2021).

Stratigraphic correlation indicates that the shallow drilled intervals (0–38 m; A in Fig. 7d) corresponding to Unit 1 are characterised by weakly under-consolidated, high-porosity sediments. Deeper intervals (42–142 m; B in Fig. 7d) corresponding to Unit 2 represent more compacted and stratified levee and channel-related deposits. The drilled MTD interval (MTD in Fig. 7d) is very thin and might correspond to the outcrop of the MTD below OBS 3005.

Measured physical properties from the boreholes, including high porosity, low shear strength, and variable gas content (Riedel et al., 2020, 2021) are consistent with the low P-wave velocities and elevated $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratios observed in the seismic data (Figs. 5 and 7). These agreements provide independent support for the seismic unit interpretation.

These constraints provide independent validation for the seismic interpretation of lithological units and reduce ambiguity in linking seismic velocity structure to sedimentary boundaries. However, due to spatial separation between boreholes and seismic profiles and inherent small-scale heterogeneity, the correlation remains locally constrained and does not resolve fine-scale variability beneath the seismic data resolution.

To quantitatively constrain these observations, the velocity and $V_{\mathrm{p}}/V_{\mathrm{s}}$ distributions are evaluated against global marine reference trends (Hamilton, 1979a, b). In terrigenous marine sediments, P-wave velocities typically increase from ∼1500–1600 m s⁻¹ near the seafloor to >2000 m s⁻¹ at greater burial depth, while $V_{\mathrm{p}}/V_{\mathrm{s}}$-ratios decrease from >10 in very soft sediments to ∼2.6 in fully compacted sequences at a depth of ∼1 km. Sandy sediments exhibit similarly high near-surface ratios (>6) but steep near-surface gradients, whereas carbonates and basalts define distinct lower-ratio endmembers (∼1.8–2.1).

Within this framework, the observed seismic units can be quantitatively positioned along the compaction trend. Unit 1 (V_p≈1510–1530 m s⁻¹; $V_{\mathrm{p}}/V_{\mathrm{s}}\approx \mathrm{10.6}$–6.8) lies within the upper range of consolidated marine sediments, consistent with high porosity and low effective stress. Unit 2 (V_p≈1580–1780 m s⁻¹; $V_{\mathrm{p}}/V_{\mathrm{s}}\approx \mathrm{5.8}$–3.7) represents an intermediate mechanical state consistent with partially compacted or reworked sediments, in agreement with the CSF. Unit 3 (V_p≈1760–1900 m s⁻¹; $V_{\mathrm{p}}/V_{\mathrm{s}}\approx \mathrm{3.7}$–2.2) approaches the lower end of fine-grained marine sediments but remains above the fully normally consolidated reference conditions. MeBo200 borehole measurements provide a higher-resolution constraint on physical properties at decimetre to meter scale compared to the seismic datasets. P-velocities range between ∼1400 and 1700 m s⁻¹, with the trend to increase within 140 m depth (Riedel et al., 2020). Locally elevated values greater than 1800 m s⁻¹ are observed within specific stratigraphic intervals. These values are consistent with the velocities observed in the upper V_p model (Fig. 7).

The upper sedimentary drape (Unit 1) shows well-stratified, high-amplitude reflections but relatively low inferred shear strength. Intercalated MTDs, such as the lens-shaped body M1 on profile A (Fig. 8), record past slope failures that infilled depressions near the headwall (e.g., Schnellmann et al., 2006). On profile B, disturbed high-amplitude reflections confined between faults F2 and F3 (Fig. 9) indicate downslope deformation adjacent to the sidewall.

Figure 9 Seismic section (upper panel) and interpretation of high-resolution MCS profile B. Gas-charged zones are highlighted in light blue, while small-scale MTDs, including M2, M3, and chaotic seismic facies (CSF2), are marked in yellow. Solid red lines indicate BSRs. The dotted red line marks the theoretical base of the present-day GHSZ, based on a bottom-water temperature of 9 °C and a geothermal gradient of 22 °C km⁻¹. The dotted orange line shows the inferred base of the GHSZ during the Last Glacial Maximum (LGM), assuming a bottom-water temperature of 4 °C (both after Badhani, 2016). White stars indicate gas flares observed during cruise MSM-34 (Hillman et al., 2018). Boxes indicate enlarged sections shown in Fig. 10.

The chaotic seismic facies in Unit 2 (Figs. 8 and 9) is characterised by discontinuous reflections and irregular geometry. Their confinement within the slump scar suggests a local failure origin associated with the S2 slope failure event, followed by deposition of Unit 1. It is interpreted as acoustic turbidity caused by sediment deformation linked to the S2 failure or even the remnants of this event (e.g., Vanneste et al., 2014; Chand and Minshull, 2003; Judd and Hovland, 1992). Semi-transparent MTDs with polarity changes between top and reversed at the base (e.g., M2, M3 on profile B; Fig. 9) indicate vertical density contrasts due to sediment mixing during material transport (Riedel et al., 2020; Mulder and Cochonat, 1996). Together, these deposits record a repeating cycle of failure, remobilisation, and redeposition.

The downward-dipping reflector (DDR) marks a laterally continuous seismic boundary that likely reflects a contrast in lithology, compaction-state, or pore-fluid content. While its exact physical origin cannot be uniquely constrained, its continuity and position suggest a stratigraphic or diagenetic interface separating mechanically distinct sediment packages.

Underlying sediments (Unit 3) represent older levee and overbank deposits formed during active canyon sedimentation. These are dominated by fine-grained clays and silts with minor sandy interbeds (Riedel et al., 2020; Bohrmann et al., 2018). Although compaction increases velocities with depth, the P-wave values remain below normal consolidation trends (Hamilton, 1979a), confirming continued under-consolidation.

This effect is consistent with high-water-content clay-rich sediments, which exhibit low shear strength, cohesion and high compressibility (Ivanik et al., 2025; Ballas et al., 2018; Bartetzko and Kopf, 2007; Keller, 1982), making them sensitive to external forcing and pore-pressure changes (Mitchell and Soga, 2005). Mechanical contrasts between silt- and sand-rich layers generate local stress concentrations at stratigraphic interfaces due to contrasts in stiffness, compaction behaviour, and pore-pressure response, which promotes strain localisation and potential failure planes (Ivanik et al., 2025). The stratigraphic architecture, levee deposits capped by hemipelagic deposits and overprinted MTDs, generates a heterogeneous depositional environment.

Normal faults (F1–F3) mapped in the upper ∼100 m (Figs. 8 and 9) penetrate only shallow units but represent structural discontinuities that may act as local zones of mechanical weakness and potential fluid migration pathways (e.g., Hillman et al., 2018), particularly where overpressure develops. Sand-rich interbeds and chaotic facies (CSF) introduce mechanical heterogeneity, reflecting both past instability and preconditioning for future failure under stress redistribution or pore pressure changes.

Overall, the sedimentary architecture of the S2 slope area is characterised by fine-grained, under-consolidated levee and overbank deposits with intercalated MTDs and sand-rich layers. High $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratios, low seismic velocities relative to normal compaction trends, and borehole-derived physical properties consistently indicate a slope system that remains mechanically weak and highly sensitive to external forcing.

Figure 10 Enlarged seismic sections highlighting key structural and fluid-related features within profiles A and B. Left panel: Structural deformation features showing (a) mass-transport deposits (MTDs) and chaotic seismic facies (CSF, yellow semi-transparent zone) from profile A; (b) faults and CSF together with observed gas-flares (Hillman et al., 2018) on profile B; and (c) MTDs along profile B. The faults are interpreted as zones of structural deformation and potential fluid pathways. Right panel: (d) Fluid-related features exhibiting the downward-dipping reflector (DDR; yellow arrows), the seismic chimney (SC, semi-transparent light blue area), and associated high-amplitude reflections beneath the DDR in profile A. The SC is marked by a velocity pull-down in the upper part of the seismic anomaly. (e) displays the bottom-simulating reflectors (BSRs; red arrows) and underlying high-amplitude anomalies in profile B. The high-amplitude reflections are interpreted as indicators for free gas accumulation. Boxes in Figs. 8 and 9 indicate the location of the enlarged sections shown here.

5.2 Gas Migration and Hydrate Dynamics

Gas migration and hydrate processes in the Danube Fan are closely linked to the fine-grained under-consolidated sedimentary architecture described above. Seismic observations reveal a dynamic interplay between free gas occurrence, focused fluid pathways, and transient gas hydrate stability.

Seismic anomalies, including polarity reversals, abrupt amplitude changes, and reflector geometries consistent with pull-down effects, indicate heterogeneity related to gas-bearing fluids. P-wave velocity variations further support the interpretation of gas-related fluid effects, with locally reduced velocities of up to ∼80–190 m s⁻¹ in zones associated with seismic anomalies. Such reductions are consistent with the presence of low gas saturations in porous sediments, which typically reduce P-wave velocities by ∼5 %–15 % in comparable marine systems (e.g., Lee, 2008; Sava and Hardage, 2009).

Beneath OBS 3004, the MCS section shows a local downward deflection of reflections (Fig. 8), which we describe as a velocity pull-down within the seismic imaging. This reflector distortion coincides with a vertically focused zone of disrupted and low-coherency reflections in Unit 2. We refer to this combined expression as a potential seismic chimney (SC). The pull-down is primarily identified in the MCS image; its presence in the forward-modelled velocity structure is less constrained due to the limited OBS offsets, sparse instrument spacing, and the locally restricted modelling approach beneath each station. Therefore, we cannot exclude that parts of the observed distortions arise from structural complexity or imaging limitations rather than a true low-velocity fluid pathway. Nonetheless, the spatial correlation between (i) the downward-deflected reflections, (ii) the vertically disrupted reflection zone, and (iii) the previously documented local gas flares in close proximity to OBS 3004 suggests that vertical fluid migration is a plausible explanation.

If this vertical anomaly represents a focused fluid-migration structure, its characteristics would be consistent with features commonly linked to local overpressure build-up (Judd and Hovland, 2007) or hydrate dissociation (Hillman et al., 2018). Its vertical extent exceeds the present-day base of the GHSZ, which is compatible with scenarios in which upward fluid migration was facilitated by hydrate dissociation following post-glacial warming since the Last Glacial Maximum (Riedel et al., 2020, 2021). The anomaly appears to be capped by a sealing layer, potentially due to insufficient gas pressure to propagate further. However, while these processes provide plausible mechanisms, the available data do not allow to uniquely resolve the origin or evolution of the observed feature.

Alternative explanations for the observed SC and associated velocity pull-down, such as imaging artefacts, local structure complexity, or stratigraphic heterogeneity, cannot be fully excluded given the limited acquisition geometry and the inherent smoothing of the forward modelling approach. However, several observations argue against a purely artefactual origin. Seismic chimneys of comparable geometry and expression have been widely documented in the Danube Fan and adjacent Black Sea region, where they are consistently associated with fluid migration indicators such as gas flares, polarity reversals, and disrupted reflectivity patterns (Hillman et al., 2018; Nasif et al., 2020; Özel et al., 2022; Riedel et al., 2021). In addition, the vertical continuity of the anomaly and its spatial association with disturbed reflectors and amplitude anomalies are characteristics of focused fluid-flow systems rather than isolated stratigraphic features. While local stratigraphic variations or structural heterogeneities may contribute to parts of the observed signal, the integrated seismic expression is most consistent with a vertically focused fluid-migration system (e.g., Judd and Hovland, 2007).

Additional seismic anomalies include polarity-reversed reflections beneath the downward-dipping reflector (DDR), and a drop of ∼80 m s⁻¹ in the V_p (Fig. 5a; Table S1 in Sect. S3). These polarity reversals are consistent with the presence of free gas accumulation in porous, sand-rich horizons sealed by low-permeability silts and clays (e.g., Riedel et al., 2020; Sava and Hardage, 2009). This is further supported by polarity-reversed, high-amplitude horizons on top of the S2 slope failure site that indicate previous gas charging. Such stratigraphically confined accumulations highlight the role of lithological layering and permeability contrasts in controlling gas distribution.

During cruise MSM-34, several active gas flares were detected near the faults F2 and F3, along the headwall, and the southeast sidewall, including sites only ∼40 m from OBS 3004 (Hillman et al., 2018; Bialas et al., 2014). Additional gas flares within the slump scar were revealed at water depths shallower than the predicted GHSZ (Fig. 2), and coincide with vertical fluid escape structures, such as seismic chimneys and disrupted, gas-bearing horizons (Hillman et al., 2018). Even though there is no clear seismic evidence along the profiles analysed confirming faults as active conduits, suggesting that gas escape is more likely episodic and localised similar to the situation at Hydrate Ridge off Oregon (Crutchley et al., 2013) and Formosa Ridge off Taiwan (Kunath et al., 2022). As gas accumulates in sand-rich horizons beneath low-permeability seals, the resulting gas overpressure can locally exceed hydrostatic pressure and capillary entry pressures, enabling episodic migration of gas towards shallower sediments (Løseth et al., 2009; Judd and Hovland, 2007). The microbial degradation of organic-rich sediments supplied by hemipelagic and turbiditic input is the most plausible gas source (Popescu et al., 2006; Jørgensen et al., 2004; Winguth et al., 2000). Thermogenic contributions from deeper reservoirs, common in other deltaic hydrate provinces (such as the Mississippi Fan; Madof et al., 2019; Sassen et al., 2001), are not supported within the Danube Fan (Starostenko et al., 2010). Gas hydrate dissociation could contribute to an episodic increase in gas concentrations.

BSRs with reversed polarity are present in profile B (Fig. 9) and can be interpreted as a direct indicator for the presence of gas hydrates (Boswell et al., 2016). However, their discontinuous character, combined with observed velocity reduction of up to −190 m s⁻¹ across the upper BSR (Figs. 7 and 9, Table S2 in Sect. S3), indicates that it primarily reflects free gas beneath the base of the GHSZ rather than significant hydrate accumulation (e.g., Andreassen et al., 2003).

Recent numerical modelling studies of the Danube paleo-delta hydrate system by Gupta et al. (2024) provide valuable context for the transient and patchy nature of gas hydrate stability observed here. Their simulations highlight the sensitivity of hydrate formation and dissociation to sedimentation sequences and methane phase transition kinetics, which control the dynamic behaviour of hydrate stability in this setting. This explains the observed discontinuous and stacked BSRs, attributed to narrow windows of hydrate formation and dissociation rates that limit hydrate accumulation and promote the coexistence of free gas beneath the GHSZ. Such dynamic processes are consistent with recent modelling emphasising how sedimentation rates and methane phase transitions drive transient hydrate systems (e.g., Iemelianov et al., 2025). Together, these findings support the interpretation of a hydrate-poor but gas-active system undergoing ongoing re-equilibration since the LGM, prone to episodic venting and overpressure.

Local phase reversals at the intersection of BSRs and stratigraphic reflectors (Fig. 9) suggest that free gas may migrate into sand-rich layers. Nevertheless, the limited extent of these accumulations implies that gas overpressure is locally confined, restricting gas migration to discrete pathways rather than widespread upflow (Judd and Hovland, 2007).

The thermally derived current base of the GHSZ lies ∼ 20 m shallower than the observed BSR, suggesting a state of transient disequilibrium in the hydrate system (Burwicz-Galerne et al., 2024; Riedel et al., 2021). This discrepancy suggests that formerly trapped gas has migrated upwards, recharging shallower sediments. The transient disequilibrium state indicated by the observed BSR depth being deeper than the thermally derived GHSZ base aligns well with Gupta et al. (2024), who highlight how gas migration upward, triggered by hydrate dissociation and sediment loading, can cause dynamic fluctuations in hydrate stability and gas saturation over geological timescales. Such transient behaviour likely drives episodic gas venting and overpressure build-up observed in seismic chimneys and fluid escape features, emphasizing the system's sensitivity to recent environmental changes. The presence of multiple former BSRs preserved in the stratigraphy supports a history of dynamic hydrate re-equilibration since the LGM (Riedel et al., 2021; Grevemeyer and Villinger, 2001) or rapid deposition of levee deposits next to the channel as reported elsewhere in the Danube Fan (Zander et al., 2017).

Controlled-source electromagnetic (CSEM) and drilling results indicate very low hydrate saturations, generally close to 0 % within the GHSZ and <10 % just at the interface of the BSR (Burwicz-Galerne et al., 2024; Bialas et al., 2020a; Schwalenberg et al., 2020). Although slightly elevated P-wave velocities were observed (up to 1970 m s⁻¹ above the BSR), the absence of corresponding V_s anomalies suggests only a small fraction of gas hydrates is present and disseminated within the sediment or patchy pore-filling (Priest et al., 2009; Yun et al., 2005). Importantly, within the resolution and uncertainties of our dataset, the lack of S-wave velocity anomalies indicates that hydrates are not acting as cementing phases at grain contacts. Instead, they occur as disseminated or patchy pore-filling accumulations that do not significantly alter bulk shear strength. This implies that hydrate presence does not directly stabilise the sediment framework, and hydrate dissociation has no intrinsic weakening effect. Rather, the impact of hydrates on slope stability in this setting is indirect, operating mainly through dissociation-driven gas release and pore-pressure increase. This form of gas hydrate distribution does not significantly alter the bulk elastic properties of the host sediment. Gas release during MeBo200 drilling near the BSR (Bohrmann et al., 2018) further indicates hydrate destabilisation at shallow levels.

Gas migration in the S2 slope of the Danube Fan is controlled by stratigraphic and structural discontinuities, with sand-rich layers acting as preferred migration pathways (Gupta et al., 2024). Vertical transport is indicated by the seismic chimney, but it cannot be determined which impact such features have due to the limited observation. Below the SC, seismic reflections remain continuous and undisturbed, suggesting that fluid migration is laterally confined and does not penetrate deeper stratigraphic levels. Free gas accumulations are locally trapped beneath low permeability seals and seem to periodically discharge through venting, while hydrate presence is minimal and dynamically re-equilibrating in response to post-glacial conditional changes. Despite the presence of BSRs with reversed polarity, their discontinuous nature and limited hydrate saturation (<10 %) suggest a predominantly hydrate-poor system where free gas accumulates beneath the base of the GHSZ in porous sand horizons sealed by low-permeability layers. This hydrate-poor but gas-active environment promotes localised overpressure and episodic venting, contributing to slope instability and complex methane migration pathways within the Danube Fan.

5.3 Geohazard and Slope Instability

The combination of under-consolidated sediments, focused gas migration, and weak mechanical layering creates a slope system prone to failure in the Danube Fan. Results from the OBS-MCS integrated analysis, MeBo200 drilling, and core analysis support this assumption.

The S2 slope failure zone is bounded by a prominent horse-shoe shaped headwall (Fig. 2), documenting a major mass-wasting event (Hillman et al., 2018). Within this zone, chaotic seismic facies and multiple MTDs record past slope failure episodes (Riedel et al., 2020; Vanneste et al., 2014). Although shallow faults (F1–F3) only penetrate the upper ∼100 m, they represent structural weaknesses that may act as local failure initiation points, particularly where overpressured gas or fluids accumulate (Hillman et al., 2018).

Potential triggers include (1) rapid sedimentation loading during phases of canyon activity, (2) earthquakes, (3) overpressure development through gas charging and hydrate dissociation, and (4) reduction in effective stress due to excess pore pressure in clays and silts (Gupta et al., 2024; Urlaub et al., 2013; Mitchell and Soga, 2005). These triggers range from gradual preconditioning by under-consolidation to episodic external forcing. The role of gas hydrate dissociation in promoting localized overpressure and slope instability has been increasingly recognized in recent numerical and thermodynamic models (Iemelianov et al., 2025). Our integrated seismic and drilling evidence supports this mechanism as a key driver in weakening the under-consolidated sedimentary sequence, thereby enhancing slope failure susceptibility.

Seismic velocity profiles consistently show lower-than-expected P-wave velocities compared to normally consolidated reference curves (Hamilton, 1979a), underscoring prevailing under-consolidation throughout Units 1–3. This condition is confirmed by drilling, which revealed high porosities and clay-rich levee deposits with low undrained shear strength (Riedel et al., 2020; Bohrmann et al., 2018). Such sediments are prone to strength reduction upon external forcing (Gupta et al., 2024; Ballas et al., 2018; Keller, 1982), amplifying slope failure instability.

Localised free gas accumulations, as evidenced by polarity reversals and velocity pull-downs (Riedel et al., 2020; Sava and Hardage, 2009), further contribute to elevated pore pressures. Chimneys and disrupted horizons demonstrate focused zones of fluid migration, redistributing overpressure vertically and laterally across the slope area (Løseth et al., 2009; Judd and Hovland, 2007). Overpressure conditions in such weak sediments reduce the effective stress and can lead to failure propagation across mechanically layered sequences.

Integrated OBS-MCS seismic data and MeBo200 ground-truth data provide complementary evidence for slope instability. OBS results show shallow velocity anomalies associated with chimneys and DDR-related gas accumulations, while cores confirm the presence of mousse-like, homogeneous silty clays in chaotic facies (Riedel et al., 2020; Bohrmann et al., 2018). These deposits lack internal layering and exhibit reduced shear strength, consistent with mass-transport deposits and disturbed slope material. This suggests that localised material and porosity contrasts, rather than gas accumulations, are the key factors regarding slope stability. Together, these data point to a slope system, where stress redistribution and under-consolidation are critical factors controlling failure dynamics.

The geohazard potential of the S2 slope is high, based on the combined effects of (1) weak, under-consolidated clays and silts; (2) stratigraphic discontinuities and shallow faults acting as planes of weakness; and (3) localised overpressure from gas accumulation and hydrate disequilibrium. These conditions suggest that even moderate external triggers, such as sediment loading or regional seismicity, may be sufficient to induce slope failure.

5.4 Comparative Fluid Systems in Continental Margin Fans

Deltaic deep-water fan systems elsewhere show closely comparable controls on under-consolidation and mechanical weakness. The Danube Fan is situated on a passive continental margin of the northwestern Black Sea (Fig. 1), shaped by Pleistocene fluvial influx and hemipelagic sedimentation under a tectonically quiescent regime, with no evidence of significant deformation (Popescu et al., 2001; Winguth et al., 2000). Gas-migration indicators and hydrate behaviour in other systems place the Danube Fan as an endmember case.

The Amazonas Fan, also on a passive margin, is driven by massive sediment supply from the Amazon River, producing thick channel–levee complexes and frequent MTDs (Maslin, 2009; Cobbold et al., 2004). Fine-grained levee deposits (mostly silts and clays) promote under-consolidation and reduced shear strength (Ketzer et al., 2018), while continuous discharge and channel migration generate lithologic heterogeneity and local instability (Souza et al., 2025; Maslin, 2009; Audet, 1998; Flood et al., 1991). Thermogenic-derived hydrocarbons increase the fluid overpressure (Cobbold et al., 2004). BSRs, chimney structures, and widespread venting from within the GHSZ show an active hydrate-free-gas system and transient hydrate stability, with direct hydrate recoveries reported (Praeg et al., 2022; Ketzer et al., 2018). Fault-controlled gas migration forms gas-charged zones with low seismic velocity and acoustic wipeouts, which correspond to weak layers (Ketzer et al., 2018). Gas hydrate concentrations can strengthen the levee deposits initially (e.g., Priest et al., 2009; Yun et al., 2005). Evidence for gas hydrates within the GHSZ is largely inferred from cores presenting different gas hydrate distributions at vent-proximal sites and levee settings (Ketzer et al., 2018). While at vent-proximal sites, nodular to massive types occur, gas hydrates within the levee deposits are rather disseminated, pointing to low-concentrations (Ketzer et al., 2018). But the dissociation of these low concentrated gas hydrates can reduce the shear strength and raises the pore pressure within these deposits, destabilising slopes (Praeg et al., 2022; Ketzer et al., 2018). Those zones are commonly associated with MTDs, indicating past slope failures. The combined effects of hydrate dissociation, pore pressure buildup, and gas-related weakening make the Amazonas Fan particularly prone to submarine landslides, even on gentle slopes.

The Krishna-Godavari Basin is located on a rifted passive margin, where Early Cretaceous rifting superimposed on older Gondwana-age graben systems created a horst-graben architecture (Radhakrishna et al., 2012; Rao, 2001). Thick, rapidly deposited sediments form sand-rich channel-levee deposits embedded in a mud-rich surrounding (Collett et al., 2008, 2019; Kumar et al., 2014). Rapid sediment accumulation and limited dewatering lead to low shear strength and localised overpressure zones, especially in the graben structures (Collett et al., 2019; Chatterjee et al., 2011; Riedel et al., 2010). Inside the Krishna-Godavari Basin, occasionally gas hydrates occur as fracture-filling hydrates in clay-silt-rich sections (Collett et al., 2014), but predominantly as laterally extensive, sand-hosted, pore-filling systems within the GHSZ (Collett et al., 2019). Inside the sandy levee deposits the gas hydrate concentration is locally very high (>50 %, occasionally ∼80 %), with free gas accumulations beneath (Collett et al., 2008, 2019; Riedel et al., 2010). These conditions create strong contrasts between the hydrate-bearing sands and the weak surrounding mud, enhancing slope failure potential. Large-scale slope failures, including an ∼11 km³ event, confirm the high geohazard risk (Dewangan et al., 2025). Together, these factors make the Krishna-Godavari Basin a geohazard-prone margin where sediment instability is closely linked to complex structural settings, rapid sedimentation, and dynamic gas hydrate systems.

The Mississippi Fan, located in the Gulf of Mexico, is strongly influenced by salt tectonics; evaporites deform extensively to form salt diapirs and structural pathways (Boswell et al., 2012; Sassen et al., 2001). The fan lobes are the major structural elements within the fan system (Twichell, 2009). They comprise a sandy channel base, silty-clayey overbanks, sheeted sands, with repeated levee failures linked to fluid migration and undrained loading of low-permeability muds (Damuth and Olson, 2015; Twichell, 2009; Weimer, 1990). The Gulf of Mexico Gas Hydrate Joint Industry Project (JIP) Leg II confirmed high-saturation pore-filling hydrates in thin-bedded to stacked turbidite sands (e.g., Haines et al., 2017). The hydrate saturations are locally up to ∼80 % of the pore space near the base of the GHSZ and are sourced largely by deep thermogenic systems and focused along salt-related faults (Madof et al., 2019; Boswell et al., 2012). In adjacent clay-rich sections, gas hydrates occur primarily as fracture-fill, mostly veins, rather than pore-filling (Lee and Collett, 2012). Dissociation or free-gas build-up beneath sealing layers can reduce effective stress reduction and destabilise slopes (e.g., Madof et al., 2019). These factors make the Mississippi Fan a complex and potentially unstable system where salt tectonics, fluid migration, and high-saturation gas hydrates significantly influence slope stability.

Taken together, these systems define a spectrum of hydrate–slope interactions. The Amazonas, Krishna–Godavari, and Mississippi fans are hydrate-rich, with sand-prone channel–levees and structural pathways concentrating fluids and free gas. In these settings, hydrate formation may initially enhance sediment strength, but dissociation typically drives pore-pressure build-up and slope failure. Among them, the Amazonas Fan is the most prone to instability due to widespread hydrates, free gas, and fault-focused migration. The Krishna–Godavari Basin follows, where rift-related structures, rapid sedimentation, and hydrate-bearing sands interact to produce frequent failures. The Mississippi Fan occupies an intermediate position, with high-saturation hydrates in turbidite sands and salt-related pathways that localise overpressure and slope instability.

In contrast to these, the Danube Fan represents a hydrate-poor but gas-active endmember. Fine-grained, low-permeability levee and lobe sediments host patchy, low-saturation disseminated hydrates and shallow free-gas accumulations beneath BSR-like horizons. Here, instability is governed primarily by under-consolidation and localised gas overpressure, with hydrates exerting a secondary influence. This highlights a distinct geohazard mode in which weak mud-dominated strata and diffuse gas charging control slope failure, rather than hydrate dissociation.

In summary, the comparison underscores that while hydrate-rich systems such as the Amazonas, Krishna–Godavari, and Mississippi fans face slope instability largely through hydrate dissociation and focused fluid migration, the Danube Fan demonstrates how slope failure can arise in hydrate-poor but gas-active margins. Across all systems, the hazard mechanism, weak layers, excess pore pressure, and external triggers remain globally consistent, but the balance between hydrates, gas flux, and sediment properties dictates local failure behaviour and geohazard potential.

6 Conclusion

This study integrates high-resolution OBS-MCS seismic imaging with drilling and coring data to investigate sedimentary processes, fluid dynamics, and slope stability in the Danube Fan, northwestern Black Sea. The results demonstrate how weak sediment properties, localised fluid migration, and transient hydrate systems interact and affect the slope architecture and geohazard potential.

First, the S2 slope is dominated by rapidly deposited, under-consolidated mud-rich sequences that retain excess pore pressure due to incomplete dewatering. Elevated $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratios, low shear strength, and recurrent MTDs highlight how weak and heterogeneous layers precondition the slope to repeated failure.

Second, gas migration and hydrate behaviour exert a critical but have a secondary influence. Seismic anomalies, including chimneys, polarity reversals, and velocity pull-downs, reveal focused pathways for free gas, while hydrate occurrences are sparse, discontinuous, and low in saturation. Disequilibrium between observed BSRs and the predicted base of the GHSZ suggests a system still adjusting to post-glacial boundary conditions. In this hydrate-poor but gas-active margin, slope instability is governed more by gas overpressure within weak muds than by hydrate dissociation.

Third, the geohazard potential of the Danube Fan arises from the combined effects of under-consolidation, structural discontinuities, and localised overpressure. Within the resolution of our dataset, the lack of S-wave velocity anomalies indicates that hydrates do not act as cementing phases at grain contacts. Instead, they occur as sparse, patchy pore-filling accumulations that do not significantly alter bulk shear strength. Even moderate external triggers, sediment loading, earthquake events, or minor hydrate destabilization can initiate slope failure, underscoring the sensitivity of this margin.

Placed in a global context, the Danube Fan defines one end of a spectrum of hydrate–slope interactions. Hydrate-rich systems such as the Amazonas, Mississippi, and Krishna–Godavari fans are dominated by sand-hosted pore-filling hydrates and fracture-related accumulations that destabilise upon dissociation. By contrast, the Danube Fan is characterised by fine-grained lithologies, low-saturation disseminated hydrates, and diffuse microbial gas charging. Instability here reflects weak mud-dominated strata and gas overpressure rather than hydrate dissociation.

In conclusion, the Danube Fan exemplifies how sedimentary architecture, gas migration, and hydrate disequilibrium jointly govern slope stability in delta-fed margins. As a hydrate-poor but gas-active endmember, it advances understanding of Black Sea geohazards and provides a critical analogue for other under-consolidated, fine-grained continental margins in the early stages of hydrate system development.