Black Sea Methane Flares From the Seafloor: Tracking Outgassing by Using Passive Acoustics

The Black Sea bottom is well known to be earth’s largest anaerobic methane source, hosting a huge amount of cold seeps releasing significant volumes of methane of both thermogenic and biogenic origin. Taking into account the well-known effects of methane concerning global warming, including the warming up of the oceans, an effective monitoring of its output from the Black Sea is nowadays an essential target for interdisciplinary studies. We discuss the results achieved during monitoring campaigns aimed to detect and track methane flares from the seafloor of the Romanian sector of the Black Sea, in order to better constrain the possible mechanisms responsible for its injection from the marine sediments, through the water column, into the atmosphere. In the mainframe of the ENVRI-Plus project, we deployed a multidisciplinary seafloor observatory for short, mid and long time monitoring and collected samples of the water column. The multidisciplinary seafloor observatory was equipped with probes for passive acoustic signals, dissolved CH₄ and chemical-physical parameters. The collected data showed a high concentration of dissolved methane up to values of 5.8 micromol/L. Passive acoustics data in the frequencies range 40–2,500 Hz allow us to discriminate different degassing mechanisms and degassing styles. The acoustic energy associated with gas bubbling is interpreted as a consequence of the gas dynamics along the water column while the acoustic range 2–20 Hz reveals vibration mechanisms generated by gas dynamic’s along the cracks and inside the sediments.

Introduction

Methane is the second major contributor to the greenhouse effect after carbon dioxide (Shindell et al., 2009). On a 100 years time scale, it holds a global warming potential that is 20–40 times higher than CO₂ (Gentz et al., 2013). Methane is generated in sediments and sedimentary rocks by the natural destruction or microbial mediated degradation of pre-existing organic matter, (Judd et al., 2002; Judd 2003). The generation of thermogenic methane occurs when complex (long-chain organic) molecules are broken down by high-temperature and high-pressure conditions, at depths typically exceeding 1 km (Schoell 1988; Whiticar 1990; Floodgate and Judd 1992). Once generated, methane migrates towards the surface (land or seabed), although some becomes trapped to form natural gas accumulations, and some other is sequestered by the formation of gas hydrates, ice-like mixtures of water and gas. Gas hydrate molecules are trapped within a cage-like framework of hydrogen-bonded water molecules, formed under very specific temperature and pressure conditions where sediments contain both water and an adequate supply of gas. These conditions are found within seabed sediments, where the water depth exceeds 300–500 m and the water temperature is less than 5°C (Ginsburg and Soloviev 1998; Ainslie and McColm, 1998, 1999). Gas hydrates may represent the largest methane reservoir of the planet. However, the majority of this gas is safely “locked” within the hydrates, and will remain there unless there is a change in the temperature/pressure conditions critical to hydrate stability. There is a growing concern that much higher atmospheric methane concentrations may result in the future, if a substantial amount of the huge methane pools stored in the ocean (Milkov 2004) and lake/reservoir sediments will be released (Buffett 2000). Gas seeps are known to be associated with leakage from gas reservoirs, shallow gas accumulations and from gas hydrates. However, the majority of methane escaping from deep-water hydrates is likely to dissolve in the water column rather than enter the atmosphere. As bubbles rise from the seabed, a portion of methane is lost to the hydrosphere by solution and microbial oxidation. The ability of methane bubbles to survive after crossing the water column is dependent upon the bubble size, the rise velocity and the presence of surfactants on the bubble surface (Leifer and Patro 2002). Surfacing bubble plumes are most likely to be found in shallow water (Hornafius et al., 1999). Judd (2000) suggested that a total of 0.4–12.2 Tg [CH₄]/year is emitted to the atmosphere from natural seabed seeps world-wide. Methane transfer from the seafloor to the atmosphere may become an important issue affecting global climate change (Hovland et al., 1993; Etiope and Klusman, 2002; Judd, 2004; Etiope and Milkov, 2004). Moreover, methane injection into the water column may contribute to the ocean acidification due to its oxidation although in shallow shelf areas, like the Black Sea, the CH₄ release from the seafloor has a greater potential to enter the atmosphere. The Black Sea represents an ideal site for investigating the fate of marine methane as methane flares are widespread recognized from the continental shelf to the deepest part of the basin, where cold seeps (mostly biogenic as a consequence of chemosynthetic communities at the seafloor) are continuously vented into the water column. During 2019, in the mainframe of ENVRIplus, a Horizon 2020 project bringing together Environmental and Earth System Research Infrastructures, we carried out three oceanographic cruises over the Romanian sector of the Black Sea with the aim to better constrain the possible mechanisms of methane injection into the atmosphere by detecting methane flares at the seafloor and tracking them across the water column.

By the cooperation of three European Research Institutions (INGV, GeoEcoMar, IFREMER) a multidisciplinary submarine observatory (EMSO-MedIT 001) was deployed at the depth of 120 m for a short, 5 days-long, period (D1) in an area already identified by previous active acoustic surveys and characterized by the presence of a pockmark and widespread methane flares. The Observatory was then recovered, the data downloaded and then re-deployed a few hundred meters away from the first site, for a 2 months-long monitoring period (D2). Once again, the Observatory was deployed for a 4 months–long acquisition period (D3) over an area close to the D1 site, in order to obtain a clear and long frame of that exhalative site, in the absence of any anthropic noise or external disturbance. The multidisciplinary seafloor observatory, able to operate in extreme environments up to the depth of 4,000 m, has been equipped with probes for passive acoustic signals, dissolved methane, temperature, hydrostatic pressure, conductivity, pH and turbidity. All the recorded data are stored in an internal mass storage and downloaded for analysis after each recovery.

During the first cruise, the water column was also sampled by Niskin bottles allowing us to improve the knowledge of its geochemical features. This paper accounts for the information gained by the passive acoustics coupled with the geochemical features of the dissolved gases. The results show the power of multidisciplinary submarine observatories in continuous monitoring of the gases discharges as well as the chemical-physical features of submarine environments, including the temporal changes due to anthropic and natural-induced processes. In particular, this paper poses the accent in passive acoustics as a tool to track gases escaping from seafloor focusing in mechanisms, detectable by hydrophone, generating noises and seismic vibrations, related to gas dynamics inside the sediment and along the water column.

Study Area

The study area is the NW sector of the Black Sea (Figure 1) where the composite Danube and Dnieper deep sea fans constitute a prominent morphological feature. These paleo-fans extend from the shelf-break, at about 200 m water depth, down to the abyssal plain of over 2000 m and were fed by rivers that drained into this area, notably the Danube, the Dnieper, the Dniester and the Bug (Wong et al., 1997). Nowadays, those rivers discharge an average annual total sediment of 65 × 10⁶ t, 81% of which (52 × 10⁶ t) is delivered by the Danube (Popa 1993). Such a huge deposition of organic-rich sediments especially at the shelf-edge deltas and the anoxic conditions below approximately 150 m that caused microbial degradation, led to the formation of shallow gas in the sub-bottom and therefore prolific gas seepage into the water column (Pape et al., 2008; Römer et al., 2020). Anoxic conditions below 150 m come from a persistently stratified water column due to the combination of the high freshwater supply from the numerous rivers and the inflow of highly saline Mediterranean waters via the Bosporus (Özsoy and Ünlüata, 1997). The sediments below the anoxic interface are the methane source for the Black Sea water column. Marine methane derives from biogenic and thermogenic processes occurring within those sediments and enters the water column by various pathways (Judd 2003). The gas discharge occurred as a dissolved methane flare or non-dissolved gaseous amounts of methane, escaping from sediments at the sea-bottom methane is therefore transferred to the water column either dissolved in the vented fluids or, in case of over-saturation, as gas bubbles (Dziak et al., 2018). When the gas discharge is persistent and vigorous, it leads to the formation of large dissolved methane plumes. Common degassing style consists in an intermittent low-flux emanation of sub-millimetric bubble plumes, while under specific conditions the gas could be suddenly released generating almost continuous high flux emissions, characterized by even centimetric bubbles plumes (Leifer and Boles 2005; Leifer 2010; Leifer and Culling 2010). The dissolved methane is diluted by mixing with the surrounding ocean water and it is further oxidized by aerobic methanotrophs bacteria. Only in cases where dissolved methane reaches the surface-mixed layer in concentrations above saturation, it can be transferred to the atmosphere via sea-air gas exchange (Cynar and Yayanos, 1993; Mau et al., 2017).

FIGURE 1

FIGURE 1. Location map of the study area extracted from Google Earth. In the insert figure: D1- 1st deployment; D2 - 2nd deployment; D3 - 3rd deployment; X and B - vertical cast; GC1 - 1st gravity core; GC2 - 2nd gravity core.

The selected area shown in Figure 1 (the Romanian sector of the Black Sea), hosts a large number of sites of methane-rich gases, from the continental shelf to the deepest part of the basin. It is consequently a good candidate for investigating the fate of marine methane, from the sedimentary layers throughout the water column to the atmosphere.

Methods

Multidisciplinary Seafloor Observatory

The multidisciplinary seafloor observatory (built in the mainframe of the EMSO-Medit project) is a versatile underwater system, able to operate in extreme environments up to the depth of 4,000 m. It is able to automatically record and store a large spectrum of long-term data, coming from acoustic and chemical-physical sensors, using an open hardware architecture that allows to integrate any commercial underwater instrument. The observatory is composed of three hard anodized aluminum vessels, hosting electronic boards and batteries, connected together by underwater high-pressure type cables and connectors (Figure 2). Vessels and sensors are hosted on the nylon-made, main frame structure (1 m × 1 m x 1.6 m). The set of sensors used during these cruises is composed by:

• Hydrophone: HTI 94-SSQ

• Temperature: SeaBird SBE 3F

• Conductivity: SeaBird SBE 4C

• Depth (Pressure): SeaBird SBE 29

• pH: SeaBird SBE 18

• Turbidity: STM11

• Methane: Franatech K-METS, S/N:K4

FIGURE 2

FIGURE 2. Multidisciplinary seafloor observatory EMSO MedIT-001.

Many efforts were made in order to have a modular and fully customizable hardware and software system architecture, allowing to easily set up the data collection configuration of the multidisciplinary observatory, e.g. to adapt the acquisition rate and the power consumption to the battery capacity and sensors needs.

The electronic system was designed to be able to acquire from both analog and digital signals. It consists of eight analog synchronous channels with configurable front-end (current, voltage or high impedance input), each channel can accept wide band signals (0.1 Hz–52 kHz) with 24-bit resolution. Digital inputs can communicate by serial protocols or Ethernet. The modular system also includes a dedicated power management board, allowing to selectively power each sensor or all sensors together. The special HW/SW design allows extremely low power consumption, thus obtaining extended autonomy, even with constrained battery volume resources. It is therefore possible to properly configure period and acquisition time to schedule sequences lasting even more than 1 year. Methane gas concentration was measured using Franatech K-METS, S/N:K4 sensor, calibrated to operate in the range from 100 nmol/L to 10 μmol/L. Following datasheet specifications in order to collect good quality data, prior to each acquisition the observatory electronics management software was configured to execute a preliminary 20 min warm-up, to stabilize its electronics response at the deployment conditions. It is reasonable to presume that there was no need for such a long warm-up period, since the sensor was submerged constantly during the three deployments and the relative humidity of the detector is supposed to be in a full stable condition.

The main data set used for this study was acquired during three oceanographic campaigns in 2019.

Deployment has been done by free fall from the research vessel, using a 1 ton concrete dead weight, while recovery has been guaranteed by acoustic releaser using the buoyancy effects of seven benthospheres.

The observatory was deployed and recovered three times, each time varying the distance from the flares area and the length of the observation period. The first and third deployment sites were indeed a few meters away each other, close to flares and surrounded by several gas emissions, while the second deployment site was at a distance of 275 m from the main emission area:

D1: Deployed on 3^rd April 2019 at 120 m depth - Recovered on 07^th April 2019;

D2: Deployed on 7^th April 2019 at 117 m depth - Recovered on 04^th June 2019;

D3: Deployed on 5^th June 2019 at 121 m depth - Recovered on 26^th September 2019.

Acoustic Records

Spectral Data Extrapolation

Hydrophone audio records along the three deployments revealed a variety of sounds, produced by different mechanisms that characterize distinct bands in the frequency spectrum analysis. Acoustic signals were sampled at 105 kHz with a 24-bit resolution A/D converter. The used hydrophone, based on a piezo-ceramic sensor element, was HTI-94-SSQ with a -198 dB re 1 V/µPa omnidirectional response sensitivity within 2 Hz–30 kHz frequency band, connected to a 40 dB fixed gain input channel. In order to perform long-term acquisition with a limited amount of storage capacity, audio recordings were bounded to segments of short duration. Segment length and data acquisition frequency were varied at each deployment (a total amount of 72 records of 10 s each were gathered per day during D2, while 12 daily records of 40 s each were acquired for D3). After a pre-processing down-sampling technique, that included anti-aliasing filtering for each record, two main approaches were used in order to inspect audio frames: at first, daily spectrograms of the resulting merged audio records were visually inspected in order to detect spectral signatures of the degassing process; then, each file was transformed to pressure signals, evaluating both manufacturer specified hydrophone sensitivity and the overall system gain; each acoustic record was then processed using spectral analysis, to extract power magnitude time-series for different acoustic frequencies bands of interest. Every acoustic record was down-sampled, in order to exclude high frequency induced noise, then low-pass filters were applied to remove inducted aliasing effects. Power Spectral Density (PSD) was then estimated for each band applying Welch’s method (Gupta and Mehra 2013), on a 10 s long segment window with 50% overlapping factor, multiplied by Hanning window, together with a 6-pole band-pass Butterworth filter to selectively analyze each band. Acoustic average pressure (intensity magnitude) modulation, and frequency bin peaks values were then extracted from every selected frequency range. Chronological metadata series were used to analyze the magnitude modulation along time and thus provide a relative estimation of the flares flow rate changes. For weaker and intermittent or impulsive energy events, once visually recognized, the relative segments were analyzed using a short-time Fourier transform and a threshold algorithm was implemented following Peregrine, 1994 and Li et al., 2021 in order to identify and count bubbles events, obtaining a bubble size distribution for each deployment.

The power spectrum analysis, performed on the acoustic records, revealed that the areas surrounding the observatory deployment sites radiated significant acoustic energy at a wide frequency band, apart from the ones of the typical ambient noise, as depicted by Figure 3, where PSD representative of the entire month of June is shown. Nevertheless, ambient noise energy levels often overcome weaker seep-related signals, representing a challenge in using this methodology to fully track the release process behaviour over time.

FIGURE 3

FIGURE 3. Monthly averaged Power Spectrum Density (PSD) of June 2019, depicting the typical envelope of the ambient noise background of the area D3. Most energetic spikes in the range 30–1,000 Hz are related to intermittent ship traffic and wind noise. Infrasonic frequencies highlight the seismic events related to methane seep release. Single, less energetic spots above 1,000 Hz are likely related to bubble streams.

Indeed, in a typical shallow underwater environment, ambient noise intensity regularly decreases from low frequencies with a slope varying from −8 to −10 dB/dec, notwithstanding that both the intensity and signature of spectral curves undergo variability in the presence of meteorological, biological and anthropogenic intermittent sources, especially in shallow sites (Wenz 1962; Wilcock et al., 2014). This represents the major difficulty of passive acoustic sources detection. A multiplicity of studies has been carried out to understand sources of underwater ambient noise that can be grouped in four representative frequency bands, each one related to different source mechanisms (Wenz 1962).

The low frequency band from 0.01 to 5 Hz is dominated by microseisms (seismic noise created by the nonlinear interaction of ocean waves) and represents the most energetic band. Spectral noise ranging from 5 to 200 Hz has a predominant anthropogenic source, produced mainly by shipping. Ships radiate a distinctive signature containing multiple harmonics, relative to the rotary machinery involved in the propulsion system (Wilcock et al., 2014 and references therein). In this band also wind induced noise can play a significant role in shallow sites. The central frequency, from 100 to 10,000 Hz, is therefore influenced by sea state and wind condition, where noise can be induced by rain droplets, shoaling gravity and breaking waves. Finally, spectral power above 10 kHz is produced by thermal agitation of water molecules. Also biological sources can be part of ambient noise; dolphins for example, can generate intense short-duration sounds from 100 Hz to 50 kHz.

Most of the spectral energy received by the observatory during the three deployments lies within the 40–2,500 Hz interval (Figure 4). The power spectral signature, composed of narrow and broadband components, testify the existence of different mechanisms responsible for sound production in this submarine system. In the first site, due to the proximity of the module to the methane flares, we assumed that the contribution of the bubbles vibration could represent part of the spectral signature, apart from the high anthropic noise due to the oceanographic operations that fully blanks the signals in some cases. In this study, we focused our analysis on frequencies up to 2,500 Hz and separately analyzed the audible frequency range (400–2,500 Hz) and the infrasonic frequency range (2–20 Hz), as the corresponding sound waves are reasonably produced by different physical mechanisms.

FIGURE 4

FIGURE 4. Typical spectra and frequency distribution of underwater sound sources (Erbe, 2008, adapted from; Wenz (1962)). Black arrows indicate the frequency range covered by the sounds emitted from various sources. Thick black lines indicate limits of prevailing ambient noise. Sea states unit describe different wave heights (e.g. wind at Beaufort 4 = approx. sea state 3).

Audible “Bubbles” Frequency Band

As mentioned above, in a submarine environment most of the energetic radiated noise could be produced by gas bubbles, due to the volume oscillation motion of the ‘bubbles’ wall after their formation at the vents. As demonstrated by Minnaert (1933) and Strasberg (1956), the fundamental peak frequency (the zeroth oscillatory mode) f depends on its equivalent radius r and ambient pressure P and can be analytically expressed by the equation:

$f=\frac{1}{2\pi r}\sqrt{\left(\frac{3P\gamma }{\rho }\right)}(1)$

where $[\gamma ]$ is the ratio of gas specific heat, at constant pressure and volume, and $\rho$ is the water density. Comparisons of Minnaert’s simple expression with experimental values and with theoretical expressions, which take into account the effects of viscosity, surface tension and thermal conductivity (Devin 1959), demonstrate that for a bubble radius greater than 0.001 cm this simple expression is valid within 5%, for the bubble fundamental frequency (also known as breathing mode).

Considering that the bubbling gas is mainly composed of CH₄, using Eq. (1) we can assume γ = 1.32, ρ = 1,030 kg/m³, attaining a product frequency by radius of fr ≈ 11.33 Hz m. Thus, for radius spanning from 0.47 to 5.6 cm, as reported in several studies based on active acoustics surveys, conducted in the same area even if at shallower depth (McGinnis et al., 2006), the ‘breathing-mode’ of bubbles should be in the range 2,400–200 Hz ± 5%. Coherently, our data revealed that the most energetic peaks emerging from the spectral analysis covered this range, rising even higher frequencies. The presence of peaks at lower frequencies, bounded between 60 and 400 Hz, may also be ascribed to bubble clouds, as bubbles together act as coupled oscillators and the resulting gas-water mixture region behaves with its own ‘breathing-mode’ frequency (Lu et al., 1990) mimicking bigger bubbles oscillation mode. As reported in several optical and passive acoustics studies, there could be three different types of bubbles plumes depending on flux: minor plumes characterized by low flow rate and bubble size strictly dependent on flux, major plume characterized by high flow rate showing broad bubble size distribution with small and large bubbles as well and intermediate plumes that show characteristics of both minor and major plumes. Depth and type of sediments (mud covered by centimetric sandy layer, as retrieved by gravity core performed during the cruise and also reported in letterature Commission on the Protection of the Black Sea Against Pollution/State of the Environment of the Black Sea (2001)) play also a fundamental role in bubbles formation mechanisms, release style and bubbles dimension (Leifer and Boles, 2005; Leifer 2010; Leifer and Culling, 2010).

Turbulent jet stream can generate a shift of the bubble frequencies due to overlapping of bubbles of various sizes (Vazquez et al., 2005; Leifer and Tang 2007; Wiggins et al., 2015).

Visual inspection of the spectrograms related to the collected data within 40–2,500 Hz confirmed the existence of pulsating sounds, having bandwidth in the same range of the one analytically assumed, likely related to bubbling oscillatory mechanism, mimicking the gas outflow behaviour over time. As the bubble stream bandwidth overlap the ambient noise one, this represented a concrete obstacle to extract clear bubble signals (Li et al., 2020). Since sea roughness could represent a relevant noise source, in order to evaluate this effect upon acoustic intensity level, we looked for any correlation between wind speed, also intended as a proxy of sea roughness variation, and bubbling related signatures, with the aim to exclude any mismatches. A weak or no correlation was found in the first site, while a clear wind-effect was found in the second and third deployment.

Acoustical Bubble Identification and Thresholding Algorithm

At low flow rates, the signal of bubble emission is not stationary and distinct transient events can be recognised. A time-frequency approach to discriminate between bubble sound emission and noise can be implemented. We took into account what was already proposed and successfully used by Li et al. (2021) and implemented a single bubble identification method based upon bubble characteristics recognition, that includes pulsation time interval, frequency bandwidth and energy strength. As the signature of a bubble is a transient event with an exponentially decaying coda, each discrete time series of the acquired acoustic record was firstly divided into 75% overlapping segments, where the length of each one was chosen in order to approximately match the average duration of the bubble signature and simultaneously guarantee a suitable frequency resolution. On each segment, a windowing Hanning function is applied before the use of Short Time Fourier Transform, in order to get spectral characteristics of the investigated segment. Once achieved the spectral computation, the algorithm was instructed to consider both the energy threshold and the frequency bandwidth threshold for identifying the bubble sounds. The damping sinusoid duration of a bubble, its dependence with frequency, gas composition within the bubble’s wall and depth, has been firstly studied by Devin (Devin, 1959) and further refined by Ainslie & Leigthon (Ainslie and Leighton, 2011). Experimental test taken by Walton et al., 2005 demonstrates that the calculated exponential decay of methane bubbles was characterised by a quality factor Q = 24 at atmospheric pressure. As a rough generalization, the value holds up to within about 50% of its surface value down to a depth of about 100 m. Having Q expressed, it is possible to roughly estimate the typical duration in ∼5ms for methane bubbles having resonance frequency centered at 1 kHz at a depth of 120 m. Thus, once Fourier parameters were determined, the thresholding technique was applied in order to identify spectral peaks above a certain energy threshold Th_E. After the peaks have been identified, all adjacent peaks were collected in a vector, and the difference between the minimum and the maximum frequency constituting the vector were compared with a bandwidth threshold Th_band (detailed algorithm conceptualization can be found in Li et al., 2021). If the resulting difference exceeds the bandwidth threshold, the sound event is recognised as a bubble signal, excluded otherwise. Such thresholding technique improves bubble signature recognition, as it eliminates isolated noise peaks. From Leighton (Peregrine, 1994. The Acoustic Bubble), the predicted bandwidth of a bubble signature, calculated upon the quality factor Q = 12 is 66 Hz for a bubble peaked at 800 Hz, rising up to 200 Hz for bubbles peaked at 2,400 Hz. However, as evidenced by Walton (Walton et al., 2005), a certain percentage of difference between calculated quality factor and measured quality factor of 3.5–40% place this value as a lower limit, It is worth mentioning that setting low thresholds, notwithstanding the increase of the probability to detect the majority of bubble events, will mostly affect the detection of false events, as many noise pulses may be identified as bubbles. Setting high thresholds prevents false alarms, but may filter too much, producing a dataset with reduced detection.

In this section we discuss the steps carried out to refine the algorithm thresholding techniques and the obtained results. The implemented algorithm has to be refined and “tuned” using experimental criteria in order to suppress the noise background from bubble events, paying attention to not exclude bubbles together with the noise. Therefore several threshold bandwidths (Th_band) have been investigated, moving by steps of 50 Hz, starting to 50 Hz (chosen as lower limit, due to the concepts discussed before) up to 200 Hz. In the same way, several energy Thresholds (Th_E) were examined by steps of 3 dB (corresponding to a doubling in energy). choosing a specific energy depending upon the deployment site and the degassing style. (see sections below). The Figure 5 describes the distribution matrix of the events detected along the whole period of acoustic records acquired in April during D2, and how it evolves as the bandwidth threshold increases. Increasing both bandwidth and energy thresholds produce a distribution of bubbles size that move from a power law distribution and become even more assimilable to a gaussian distribution that represents the best agreement between detection of bubbles events and false positive counting due to background noise. On the other hand, choosing higher energy and bandwidth thresholds could result in a misreading of detection at higher frequency, thus underestimating bubbles count and, as a consequence, the amount of released gas over time. Looking at the distribution, it becomes clear that the shape is insensitive with respect to the choice of the energy threshold, while it changes as the bandwidth threshold increases. The figure shows also how a higher energy threshold represents a good filter, but at the same threshold, setting a good bandwidth threshold improves discrimanationn of bubble events. Some bubble events could however be excluded, thus, in order to choose the right thresholds, it is necessary to investigate how the distribution for all pulses is affected by noise other than bubbles. For each deployment period, we generated a distribution matrix where each row represents the investigation of bandwidth threshold, considering the same energy threshold. Figure 6 describes distribution matrix of the events detected along the whole period of acoustic records acquired during D3

FIGURE 5

FIGURE 5. Bubble size distribution matrix at increasing Th_E [97–103 dB] and TH_Band [50–200 Hz] relative to April 2019. Each barchart plot shows, once Th_E and Th_Band fixed, the evolution of %vol related to detected “bubble” events vs radius size as a function of frequency. Barchart plot in the black box represents the chosen thresholds where the gaussian function describes the distribution at the best. First Barchart plot on the top right indicates a distribution described by a power law function giving an idea of not optimized threshold choices. Line plot on the right shows Number of detected bubbles and relative computed l/m vs increasing thresholds. Black square inside the line plot hilights the best agreement between false positive and bubble detection

FIGURE 6

FIGURE 6. Bubble size distribution matrix at increasing ThE 85–90 dB] and THBand [50–200 Hz] relative to the entire D3. Each barchart plot shows, once ThE and ThBand fixed, the evolution of %vol related to detected “bubble” events vs radius size as a function of frequency. Barchart plot in the black box represents the chosen thresholds where the gaussian function describes the distribution at the best. First Barchart plot on the top right indicates a distribution described by a power law function giving an idea of not optimized threshold choices. Line plot on the right shows Number of detected bubbles and relative computed l/m vs increasing thresholds. Black square inside the line plot hilights the best agreement between false positive and bubble detection

The Figure 7 highlights how the filtering technique helps step by step to improve the detection of the “bubbles” events with respect to noise. Fixing the optimal energy threshold, the resulting spectrograms depict the improvement of the performance of the detection as the bandwidth threshold increases.

FIGURE 7

FIGURE 7. Spectrogram matrix using short-time Fourier transform, where different thresholds are applied; (A) Spectrogram of raw data recorded; (B) data filtered applying Th_E = 88 dB, without bandwidth thresholding (C) data filtered applying Th_E = 88 dB, Th_B = 50 Hz, from which is possible to note the coexistence of bubbles and noise events; (D) data filtered with Th_E = 88 dB, Th_B = 100 Hz, the resulting data present still false positives; (E) Th_E = 88 dB, Th_B = 150 Hz, identified “bubbles” are clearly shown as orange dots dominantly at frequency interval 800–2,200 Hz; (F) defined Th_eng = 88 dB, Th_Band = 200 Hz as best couple of values from the distribution matrix, identified single “bubbles” are easily identified from the spectrogram.

Infrasonic Frequency Band: 2–20 Hz

Over the three deployment periods analyzed in this paper, low acoustic frequencies within the range 2–20 Hz were investigated as well in order to detect and follow the vibrational signals evolution. Such “vibrational” signals depend on a totally different acoustics mechanism with respect to the one generating bubbling noise, being the conversion into acoustic waves at the solid-liquid interfaces of seismic energy (Wilcock et al., 2014 and references therein) generated below the seafloor by dynamics of ascending fluids along cracks (Dziak and Fox Christopher, 2002, 2004; Tolstoy et al., 2002; Ohminato 2006; Dziak and Haxel, 2008, 2012; Alparone et al., 2010; Milluzzo et al., 2010; Longo et al., 2021). Searching for acoustics sources from these kinds of mechanisms results from the assumption that if there is gas in the water column producing bubbles noise (named bubbling, whose oscillatory mode is acoustically detected at higher frequencies within 200 and 2,500 Hz), there necessarily should be a gas flowing inside the sediments along a network of pathways that links the gas source area with the sea bottom. Gas flow dynamic itself makes the pathways vibrating and produces a seismic signal converted in infrasonic frequency acoustic energy at the solid-liquid interface just at the seafloor.

Vertical Casts

During the first cruise we performed vertical casts along the water column using a Rosette with 16 Niskin sampling bottles. We collected water samples at various depths in two localities (named B and X) nearby the deployment site (Figure 1 and Table 1). The results allowed us to define the geochemical features of the water column and its dissolved gases up to the sea bottom. The samples, specifically collected for the extraction of the whole gas phase for chemical and isotopic analyses, were stored into 240 ml Pyrex bottles, sealed onboard using rubber or Teflon septa and purpose-built pliers, and analyzed within 2 weeks. Details on the sampling methodology are reported in Italiano et al. (2009, 2014). During the cruise, the sampled bottles were stored upside down, keeping the necks with the rubber septa submerged in seawater until they were transferred to the laboratory. The collected seawater samples underwent laboratory procedures for chemical (concentration of dissolved gas species) determinations. In the laboratory, the dissolved gases were extracted after equilibrium was reached at constant temperature with a host gas (high-purity argon) injected in the sample bottle (see Capasso and Inguaggiato, 1998, Inguaggiato and Rizzo, 2004 and Italiano et al., 2009, 2014, for further details). The chemical analyses of O₂, N₂, CH₄ and CO₂ were carried out by gas chromatography (Agilent 7800B equipped with a double thermal conductivity detector–flame ionization detector; TCD–FID) using argon as a carrier gas. Typical analytical uncertainties were within 5%.

TABLE 1

TABLE 1. List of sites for: deployment multidisciplinary seafloor observatory (D1 - D2 - D3); sampling gravity cores (GC-01 - GC-02) and sampling seawater columns (B - X) during the first surveys.

Results and Discussion

Geochemical Features of Seawater Columns

The chemical-physical features and the analytical results in terms of dissolved gases composition of the water samples collected at various depths are listed in Table 2.

TABLE 2

TABLE 2. Chemical-physical features and dissolved gas composition of the seawater columns from the studied localities in the Black Sea (B, X). T. BDL: below detection limit, na: not analyzed.

The seawater samples collected along the water column show a chemical composition similar to the ASSW (Air Saturated Sea Water) up to a depth of 60 m. From 60 m to the bottom, the dissolved methane concentration started to increase, reaching the highest concentration in the range of 0.2 micromol/L close to the bottom at a depth between 100 and 110 m. The pH values decreased up to 7.8 value, lower than normal seawater. The dissolved carbon dioxide shows an increased concentration as well, reaching values around 0.08 millimol/L. The rise of carbon dioxide content was coherent with the pH decrease and could be related to methane oxidation processes. The Eh values exhibited values around 0 or negative close to the sea bottom accounting for enhanced anoxic conditions (Figure 8).

FIGURE 8

FIGURE 8. Cast profiles obtained by discrete sampling from niskin bottles collected at different depths. Data highlight increasing dissolved methane concentrations from the depth of 60 m to the bottom with the highest values ranging about 0.2 micromol per liter, negative values of Eh and pH about 7.78 at 110 m water depth. Yellow marker represents cast profile data from site B, blu marker represents cast profile data from site X (see Figure 1 for locations of both B and X sites).

Chemical-Physical Parameters

The continuous recording of chemical-physical parameters (temperature, conductivity, pH, hydrostatic pressure), performed by the observatory over the three deployment periods (Figure 8), showed how the main features of the bottom water did not suffer significant variations over the time. Conductivity being not T-compensated shows the same modulation over time with respect to the temperature that changes in the narrow range of 8.5 ± 0.1°C. Regarding the pH value, it ranges from 7.6 to 8, mostly persisting on lower values (7.8) with respect to the normal seawater, in agreement with the values recorded near the sea-bottom during cast profiles (Figure 9).

FIGURE 9

FIGURE 9. Chemical-physical parameters recorded by seafloor observatory from April 2019 to September 2019. The diagram shows, from the top to the bottom, variation over time of temperature (black line), pH (green line), electric conductivity (red line), hydrostatic pressure (dark blue line). D1, D2, D3 vertical black bars separate the three different deployments in the three different sites.

First Deployment (D1)

Investigating the acoustic energy either as a gas detecting tool or as a proxy of gas flow is the main goal of the passive acoustic measurements. Methane bubbles coming from flares radiate significant broadband and narrowband acoustic pressure waves. In addition, dissolved methane results in chemical anomalies detectable by the methane sensor. The first 4 days-long deployment (D1) allowed us to verify the acoustic response of an abrupt methane injection into the water column. For D1, the observatory was set up to promptly start acquisition, recording all parameters every 20 min, once reached the sea-bottom. The very first data of dissolved methane, turbidity and acoustic energy in the typical bubbling frequency range as predicted by the Minneart equation (400–2500 Hz) have shown low values that we could consider as the local background.

All those three parameters kept low values until the first planned gravity core (GC-01) was performed less than 200 m away from the site of the observatory deployment (Figure 1). After GC-01, all the measured parameters suddenly increased. The gravity core activity likely induced the stripping of gas trapped in the shallower sediments. As a secondary effect, turbidity, produced by sediments moved by gas jets and/or bubbles clouds, grew up. Dissolved methane within 24 h reached higher values around 3.5 micromol/L. The increasing trend of dissolved methane values was exactly matched with the trend of turbidity ones. The acoustic energy levels followed the trends of both turbidity and methane as well, as the enhanced high gas flow likely produced turbulent gas-jets and plumes of bubbles characterized by different sizes (from millimetric up to even centimetric) which oscillation, increased the acoustic intensity (Figure10). Spectral short-time Fourier transform analysis was performed over 512 points, corresponding to 9.7 ms observing window size, followed by a threshold denoising algorithm. This allowed us to discriminate, during loud less periods (e.s. records acquired overnight, when the survey activities were stopped), bubbles streams, as shown in Figure 11. Moreover, acoustic energy levels in infrasonic frequency ranges (2–20 Hz) showed variations due to ascending methane through the pathways inside the sediments below the seafloor. From the beginning of April 5th all parameters showed a synchronous decreasing trend. On the same day at 05:15 (UTC) the second planned gravity core (GC-02) was performed more than 400 m away from the observatory. After the GC-02, only acoustics data recorded a weak simultaneous variation related to the gravity core activity. In fact, due to both effective double distances of GC-02 from the observatory with respect the previous (GC-01) and the likely weak presence of methane still trapped in the sediments, only a little and delayed increase was recorded by acoustics and turbidity, while dissolved methane concentration showed a slight increase up to 1.9 micromol/L. All parameters went back to their initial values at the end of April 6th.

FIGURE 10

FIGURE 10. Evolution over 5 days of the parameter recorded by the seafloor observatory during D1 from the bottom variation over time of dissolved methane (dark green line), turbidity (brown line), low acoustic energy levels 2–20 Hz (dark red line), acoustic energy level 500–2,500 Hz (orange line), wind speed (red line). Dashed black vertical lines mark the time of gravity core activity (GC-01 and GC02). On the bottom, from left to right, correlation plots of Turbidity vs Dissolved methane, 500–2,500 Hz vs Dissolved methane, Turbidity vs 500–2,500 Hz and 2–20 Hz vs Dissolved methane

FIGURE 11

FIGURE 11. Detailed short-time Fourier transform (STFT) and waveform relative to a single 10-second-long audio segment acquired on 5th April, where diffuse bubble streams are recorded. STFT was calculated over 120 points, 75% overlapping factor, corresponding to 20 ms window size.

Due to coring activity D1 is therefore characterized by high flow and plume of bubbles that produce a stationary acoustic signal. Identification bubbles could be performed using the acoustic inversion method proposed by Leighton and White, 2011 but, considering the concurring high level of anthropic noise that often superimposes over the spectral bandwidth related to bubble streams, it is hard to apply. Despite all, it was however possible to track and extract the parameters described above, that showed high correlation, round-crossing all parameters one with each other, with R² factor spanning from 0.7 up to 0.86 as shown in the correlations plots of Figure 10, therefore revealing a strict dependence among the physical phenomena observed.

Second Deployment (D2)

On April 7th the seafloor observatory was recovered, and after data downloading and battery recharging operations, was again deployed at the second site (D2), 275 m to the SW of the first site. In order to collect data for a longer period the observatory was set to acquire every 40 min and record 10 seconds-long audio segments each time, obtaining a total of 72 audio records per day. Data acquisition of this deployment lasted for 2 months and showed pulsating variation of dissolved methane at the sea bottom that never reached values higher than 1.5 micromol/L. Pulsating behavior suggested a cyclic, although not-regular, disappearance of the methane emissions. This is coherent with the observations reported by other authors (Judd, 2003 and references therein, Dziak et al., 2018) that described, for similar areas, an exhalation style characterized by cyclic accumulation inside the sediment followed by low-flux methane release. With respect to D1, in which turbulent high flow gas release has been induced by anthropic activity (Gravity Core activities), in D2 gas release likely occurred every time the sediments were not able to retain pressure induced by gas accumulation (Figure 12).

FIGURE 12

FIGURE 12. Evolution over 2 months of the parameter recorded by the seafloor observatory during D2 from the bottom variation over time of low acoustic energy levels 2–20 Hz (dark red line),, %vol released during April (dark-blue line) and %vol released during May (light-blue line), turbidity (brown line), dissolved methane (dark green line),and wind speed (red line).

In such conditions, during a month long period starting from the second week of April and the beginning of May, we observed an increasing stage of energetic burst-like and intermittent broadband events, that slowly decreased until totally disappearing from May 9th on. From both visual and acoustical inspection, we are confident that such events are likely related to methane bubbles streams, as shown in Figure 13A, accounting for a sudden release of variously sized, even centimetric, bubbles. In order to investigate the bubble size distribution of the events, the thresholding analysis was performed for all the records in the mentioned period, using the thresholding algorithm described in section 3.2.3 above. As the visual inspection of the spectrogram returned energetic broadband signals, we investigated frequency windows within 400 and 1,200 Hz. The algorithm was instructed to better analyze frequency content and retrieve several distribution histograms, by the inspection of different energy thresholds Th_E, in the range [97–103] dB, by step of 3 dB, corresponding to a doubled energy. For each fixed Th_E, different frequency bandwidths were inspected, from 50 to 200 Hz. Figure 5 shows the resulting bubble size distribution matrix, highlighting how the peaks centered on frequency 500 Hz, corresponding to radius 2.2 cm, persist even though the increase of threshold values. Looking at the distribution matrix, it is possible to observe how the application of high threshold values may produce an overfiltering, excluding bubbles from the distribution dataset, especially the ones less energetic with higher frequency, while the distribution shape doesn’t change significantly. On the other hand, setting a low Th_Band, such as 50 Hz, retrieves false positive detections. Infact, as can be seen from the insert graph in Figure 12, setting low Th_Band threshold produces a local peak of detection. The investigation of energy thresholds suggests that, from the 100 dB, the slope of the curve smoothly decreases, suggesting that the detected events preserve a consistent bandwidth and are thus assimilable to bubble events. Moving forward from this point on, may result in an excessive filtering of bubbles events. Therefore, the best agreement between bubble detection and false positives is reached choosing Th_E = 100 dB, Th_B = 200 Hz, corresponding to the distribution highlighted by the solid-black square in the matrix. Once selected the best suited histogram to express bubble size distribution, the related evolution of volume percentage were compared to wind speed (as proxy of sea roughness), low frequency averaged spectral energy and chemical-physical time series of Figure 11, revealing a high coherence over time. t The correct tuning of the algorithm made therefore possible to acoustically track the energy signature related to the bubbles’ walls motion, just nearby the time corresponding with the higher values reached by dissolved methane (Figure 13B). The waveform of the reconstructed signal further highlights the chronological succession of the signals, suggesting a foregoing release of millimetric and sub-centrimetic bubbles followed by larger ones. Lower frequency content is likely produced by the oscillatory mode of the bubble cloud that surrogate a larger bubble oscillatory motion, as previously mentioned (Lu et al., 1990) or even the frequency shift produced by coalescing effect of small chained bubbles (Figure 10B) (Vazquez et al., 2005; Leifer and Tang 2007; Wiggins et al., 2015). At the same time, infrasonic signal levels, acting in pulsating mode as well, synchronously reached high intensity levels, revealing the best agreement over time with dissolved CH₄ in the water column. Also turbidity showed variation over time correlated to methane, somehow testifying the suspension of sediments due to bubbles rise dynamic (Figure 12). Comparisons of infrasonic frequency with wind time series (as a proxy of sea state) exclude any correlation. It is therefore plausible to exclude the contribution of microseismicity induced by sea roughness. It is a matter of fact that the infrasonic level rise (which energetic intensities gained the highest values of the entire period) is exactly in accordance with the bursts increase, testifying that the two parameters describe consequential mechanisms. Such anomalous events of gas release slowly exhausted until totally disappearing on May 10th, following a rather quiescent period. Since the beginning of May, a new gas emission style was recorded by the hydrophone, following seismic and methane concentration changes as well. Bubble size distribution performed on the entire D2 period (Figure 14) confirmed the change in gas release style. After a visual inspection of acoustic records, the threshold algorithm has been tuned to analyze a frequency window within 800 and 2,400 Hz, looking for, by Equation (1), oscillatory bubbles’ mode likely sourced by bubbles which radii laid in the range 0.4–1.5 cm. Therefore, applying Th_E by steps of 3 dB from 85 to 90 dB and Th_Band from 50 to 200Hz, the resulting matrix highlights how the most confident distribution is obtained choosing Th_E 88 dB and Th_Band 200 Hz (Figure 14). Once again plotting the relative Volume percentage over time we obtained the evolution of bubbles released with a high degree of coherence with respect to dissolved methane and infrasonic energy levels. The recorded bubble release process was definitely different from the previous one, being characterized by a single smaller bubble emission, whose rate slowly increased until rising its maximum during the period May 30th–June 2nd.

FIGURE 13

FIGURE 13. (A) Daily resulting spectrogram, calculated on 2048 point-sized window (38.6 ms), of merged 10 seconds-long audio records, recorded on 26th of April; the bandwidth shown is bounded to the inspected 20–2,500 Hz band, where energetic bubble bursts signals emerge, together with ship related noise. (B) Detailed short-time Fourier transform (STFT) and waveform of a single burst event recorded. STFT was calculated over 120 points, 75% overlapping factor, corresponding to 20 ms window size.

FIGURE 14

FIGURE 14. Graphs describing the D2 features resulting from the investigation. (A) Line plot describing “Number of events” and related l/min vs. threshold. (B) Bar chart plot showing the distribution chosen after the thresholds investigation. (C) Detailed short-time Fourier transform (STFT) spectrogram and waveform (Th_E = 88 dB, Th_Band = 200 Hz) relative to a single record acquired on May 15th.

Third Deployment (D3)

The third deployment site was 20 m away from the first site, already described as an extensive exhalative area. This was the longest monitoring period thus, in order to further reduce power consumption, the observatory was configured to acquire every 2 hours, obtaining a total of 12 samples per day. Audio records were set to acquire 30 s-long segments. Data recording lasted 4 months and revealed, as expected, a significant and almost continuous methane release with dissolved concentrations that in certain periods reached values up to 5.8 micromol/L (Figure 13). Higher dissolved methane values were somehow preceded by an increase of the acoustic energy in infrasonic frequencies levels with a synchronous increase of the bubbles emission rate. Bubble size distribution (Figure 6) obtained through the application of the threshold algorithm in the frequency range [ 800–2,400] Hz returned similar characteristics with respect to D2, with the best distribution obtained using Th_Band 200 Hz and Th_E 88 dB. The investigation process revealed the existence of two gaussian distributions, the main one peaked at 1.2 cm, while the second one characterized by a moda at 0.9 cm. Such evidence likely suggests the detection of almost two different degassing features, characterized by small bubbles gently released. Figure 15 shows the evolution over 4 months of volume percentage, obtained plotting the time series of the number of detected events, using the threshold values that gave the best distribution. The series is compared with dissolved methane, wind speed and infrasonic energy levels. Looking at the figure it is possible to note a temporal shift among acoustics data and methane concentration, suggesting that the deployment site is quite distant from the seeps, differently with respect to D1 and D2, and resulting in a sort of hysteresis in the chemical response on the water column. In light of this scenario, choosing a lower energy threshold might retrieve a more realistic quantification of volume of gas released, but might also include too much noise, resulting in a higher percentage of false detection as the longer distance attenuates less intense seeps. Our hypothesis could be further confirmed if we take into account the energetic levels received: the burst-like events generated in D2 april raised 120 dB, while the sound of bubble detachment never overcame 107 dB. Acoustic record signatures describe a natural and gentle exhalation style, characterized by an almost continuous release of small (millimetric) bubbles, centered in the 900–1,000 Hz range, as shown in Figure 16. During the increasing stages, the presence of larger, but still subcentimetric bubbles (peaked in the frequency range 700–900 Hz) was detected as well. The overall acoustic results describe a general low emission rate.

FIGURE 15

FIGURE 15. Evolution over 4 months of the parameter recorded by the seafloor observatory during D3 from the bottom variation over time of low acoustic energy levels 2–20 Hz (dark red line), %vol released over the entire D3 (blu line), dissolved methane (dark green line). and wind speed (red line)

FIGURE 16

FIGURE 16. (A) Example of raw daily resulting spectrogram, calculated on 2048 point-sized window (38.6 ms), of merged 30-second-long audio records, recorded on 24th of June; the bandwidth shown is bounded to 700–2,500 Hz band; white rectangular highlight single bubbles’ oscillatory signals; red rectangle details a specific bubble stream. (B) Line plot (relative to D3 investigation) describing “Number of events” and related l/min vs. threshold. (C) Bar chart plot showing the distribution chosen after the investigation made on D3 thresholds investigation. (D) Detailed short-time Fourier transform (STFT) spectrogram and waveform (Th_EE = 88 dB, Th_Band = 200 Hz) relative to a single record acquired on August 26^thc during D3.

Bubbles Volume and Rate Estimation

Bubble size distribution method proposed by Leighton allows to obtain a quantification of bubbles detection versus frequency bands. Using the Minneart equation, each frequency band corresponds to a bubble radius. Under the assumption of spherical bubble shape, it is therefore possible to convert events detected to volume of gas if we simply consider the volume corresponding to a spherical bubble with the detected radius. Hence, using each distribution, it is possible to estimate the degassing volume emitted per day. Looking at the distribution matrix of D3 Figure 6), obtained by the cumulative detected bubbles over 4 month of observation, it returns major events peaking at radius equals to 1.2 cm, corresponding to bubbles of about 7cc. Hence, if we consider the entire distribution of radii and their corresponding volume per number of detected bubbles, under the assumption of steady degassing rate we obtain the total averaged gas emitted during the D3 period. Thus, D3 bubble size distribution gives back about 0.7 L/min, which results to a total degassing volume of 1,368 L/day corresponding to 771 mol/day. In the same way we can compute the outgassing values relative to D2 (Figure 14). As already described in section 4.4, D2 is characterized by a general degassing stage, similar to the one observed in D3 during 4-month long deployment. During D2, within some limited period, particularly during April, burst style degassing events were observed. Analyzing the bubble size distribution of the entire D2 it is worth to note that, despite the shorter length of the observation temporal window (2 months long), a greater overall averaged volume of gas per day was observed, estimated at 3,054 L/day corresponding to 1721 mol/day. Moreover, analyzing the size distribution relative to burst events of April (Figures 5, 13), characterized by bubbles radii of about 2.3 cm, we obtain a further outgas value per day equals to 547 L/day corresponding to 308 mol/day. Adding this estimation to the previous value, a total value of 3,601 L/day equals to a 2030 mol/day is obtained, about three times greater than the estimation retrieved from D3. The comparison between the obtained volume estimations and the different mean concentration of the dissolved methane in both the areas, allows us to assume that, even though D3 shows a gentle emission style and released small size bubbles, albeit continuously, it shows the highest concentration values (5.8 micromol/L). Such evidence is likely due to the position of the observatory with respect to the flares emitting methane. Specifically, it may be located downstream with respect to the direction of the currents that possibly transport the gas just in the direction of the chemical sensor. Differently, despite the two combinated degassing style, continuously diffused and punctual burst-like release, D2 could be likely sited upstream with respect to the dominant currents. Hence, while acoustic records could detect bubble signatures, chemical sensor acquisitions may gain less concentration of dissolved methane (1.5 micromol/L). The lack of a current meter in our equipment represents a limit in the robustness of the above assumptions.

Conclusion

Coupling the geochemical features of the water column with the temporal changes recorded by a sea-floor multidisciplinary observatory and the spectra analysis of acoustic records we gained further information useful to better constrain the processes and the interactions governing the injection of methane from the sea bottom sediments to the water column in the Black Sea.

The Romanian sector of the Black sea, hosting a huge amount of methane of various origin, has been a perfect study area where we deployed a multidisciplinary submarine observatory (EMSO-MEDIT_001) in the main frame of the ENVRI-Plus activities. The geochemical features of the dissolved gases along the water column show a methane content which increases in the deep water layers besides lowering of pH down to 7.8 or negative Eh values.

The data of the continuous monitoring of temperature, pH, EC, turbidity, dissolved methane collected by the sea-floor observatory during short, mid and long-term deployments, depicted temporal changes induced by anthropic disturbances (coring) as well as by natural forces (e.g. wind, sea waves) and natural events (e.g. earthquakes). The information from passive acoustics allowed us to detect and describe several mechanisms related to the gas release inside the sediments and along the water column. Albeit the short duration of every single acoustic record, it was yet possible to identify different frequency ranges (within the 200–2400 Hz range) related to the energy produced by the motion of methane bubble’s walls of various sizes.

As this feature is a consequence of the emission style, we discriminated sudden and impulsive gas bursts from continuous gentle gas emission. Moreover we d roughly estimated volume of the released gas, based on bubble size distribution and the bubbles acoustic signature in the different areas. We estimated a methane flow rate in the range of 771-2030 mol/day. The acoustic technique is thus capable of identifying single bubbles in noisy environments and determining bubble size distribution and gas flux. The implemented single bubble identification technique can be performed automatically for any length of acoustic signals. This allows long-term underwater acoustic detection and quantification of different gas releases over a variety of marine environments in the world-wide ocean seafloor. The technique is applicable when bubbles are generated as a single stream or multiple streams. However, when the bubbles are generated as plumes or clouds rather than single stream, it is difficult for the technique to identify single bubbles, which may make the estimate results inaccurate and an alternative way would be the acoustic inversion method (Leighton and White, 2011). The use of multiple hydrophones in an array configuration may however improve the detection performance, especially in noisy environments, and even extend the range of observation. Furthermore, equipping the submarine observatories with bottom-currents measurement devices could result in a more complete and exhaustive information concerning the water column features. The hydro-acoustic data gathered allowed also to detect the “vibrations” of the pathways network induced by the escaping gas dynamics through the sediments below the sea bottom (in the infrasonic band 1–20 Hz), providing an additional perspective to better bind the mechanisms responsible for the generation of methane degassing.

As the methane release in the submarine environment is a cause of growing concern worldwide, due to the effects related to global warming and ocean acidification, the results we gained provide new tools to better constrain and quantify the acceleration of production and release in marine environments.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

ML and GL contributed to conception and design of the study. CC organized the database. ML and GL wrote the first draft of the manuscript. CC, SS, VR, RR, FI and DB wrote sections of the manuscript. All authors contributed to article revision, read, and approved the submitted version.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors are grateful to all crew members of R/V MareNigrum, for their precious help and professionality. All the valuable results couldn’t have been possible without the financial support of the “European Multidisciplinary Seafloor and water column Observatory” consortium (EMSO), which under EMSO-Medit and EMSO-Miur projects enabled infrastructure maintenance and data analysis. The authors would also like to tank the reviewers Cuiling Xu, Jianghui Li and the third anonimous reviewer for their precious suggestions and comments that made possible a significant improvement of the manuscript.

References

Ainslie, M. A., and Leighton, T. G. (2011). Review of Scattering and Extinction Cross-Sections, Damping Factors, and Resonance Frequencies of a Spherical Gas Bubble. The J. Acoust. Soc. America. 130, 3184–3208. doi:10.1121/1.3628321

CrossRef Full Text | Google Scholar

Ainslie, M. A., and McColm, J. G. (1998). A Simplified Formula for Viscous and Chemical Absorption in Sea Water. J. Acoust. Soc. America. 103, 1671–1672. doi:10.1121/1.421258

CrossRef Full Text | Google Scholar

Alparone, S., Cannata, A., Gambino, S., Gresta, S., Milluzzo, V., and Montalto, P. (2010). Time-Space Variation of Volcano-Seismic Events at La Fossa (Vulcano, Aeolian Islands, Italy): New Insights into Seismic Sources in a Hydrothermal System. Bull. Volcanol. 72, 803–816. doi:10.1007/s00445-010-0367-6

CrossRef Full Text | Google Scholar

Buffett, B. A. (2000). Clathrate Hydrates. Annu. Rev. Earth Planet. Sci. 28, 477–507. doi:10.1146/annurev.earth.28.1.477

CrossRef Full Text | Google Scholar

Capasso, G., and Inguaggiato, S. (1998). A Simple Method for the Determination of Dissolved Gases in Natural Waters. An Application to thermal Waters from Vulcano Island. Appl. Geochem. 13, 631–642. doi:10.1016/s0883-2927(97)00109-1

CrossRef Full Text | Google Scholar

Commission on the Protection of the Black Sea Against Pollution/State of the Environment of the Black Sea (2001). Available at: http://www.blacksea-commission.org/_publ-SOE2009.asp#_Toc225838294.

Cynar, F. J., and Yayanos, A. A. (1993). Biogeochemistry of Global Change: Radiatively Active Trace Gases. Editor R. S. Oremland (Menlo Park, CA: Chapman & Hall), 551–573.

Devin, C. (1959). Survey of thermal, Radiation, and Viscous Damping of Pulsating Air Bubbles in Water. J. Acoust. Soc. America. 31, 1654–1667. doi:10.1121/1.1907675

CrossRef Full Text | Google Scholar

Dziak, J. J., Nahum-Shani, I., and Collins, L. M. (2012). Multilevel Factorial Experiments for Developing Behavioral Interventions: Power, Sample Size, and Resource Considerations. Psychol. Methods. 17 (2), 153–175. doi:10.1037/a0026972

PubMed Abstract | CrossRef Full Text | Google Scholar

Dziak, R. P., and Fox Christopher, G. (2002). Evidence of Harmonic Tremor From a Submarine Volcano Detected across the Pacific Ocean basin. JGR Solid Earth. 107, ESE 1-1–ESE 1-11. doi:10.1029/2001jb000177

CrossRef Full Text | Google Scholar

Dziak, R. P., and Haxel, J. H. (2008). Observations of Regional Seismicity and Local Harmonic Tremor at Brothers Volcano, South Kermadec Arc, Using an Ocean Bottom Hydrophone Array. J. Geophys. Res. Solid Earth. 113, B08S04. doi:10.1029/2007jb005533

CrossRef Full Text | Google Scholar

Dziak, R. P., Matsumoto, H., Embley, R. W., Merle, S. G., Lau, T.-K., Baumberger, T., et al. (2018). Passive Acoustic Records of Seafloor Methane Bubble Streams on the Oregon continental Margin. Deep Sea Res. Part Topical Stud. Oceanography. 150, 210–217. doi:10.1016/j.dsr2.2018.04.001

CrossRef Full Text | Google Scholar

Dziak, R. P., Smith, D. K., Bohnenstiehl, D. R., Fox, C. G., Desbruyeres, D., Matsumoto, H., et al. (2004). Evidence of a Recent Magma Dike Intrusion at the Slow Spreading Lucky Strike Segment, Mid-Atlantic Ridge. J. Geophys. Res. 109, B12102. doi:10.1029/2004JB003141

CrossRef Full Text | Google Scholar

Erbe, C. (2008). Critical Ratios of Beluga Whales (Delphinapterus leucas) and Masked Signal Duration. J. Acoust. Soc. America. 124, 2216–2223. doi:10.1121/1.2970094

PubMed Abstract | CrossRef Full Text | Google Scholar

Etiope, G., and Klusman, R. W. (2002). Geologic Emissions of Methane to the Atmosphere. Chemosphere. 49 (8), 777–789. doi:10.1016/s0045-6535(02)00380-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Etiope, G., and Milkov, A. V. (2004). A New Estimate of Global Methane Flux From Onshore and Shallow Submarine Mud Volcanoes to the Atmosphere. Env Geol. 46 (8), 997–1002. doi:10.1007/s00254-004-1085-1

CrossRef Full Text | Google Scholar

Floodgate, G. D., and Judd, A. G. (1992). The Origins of Shallow Gas. Continental Shelf Res. 12 (Issue 10), 1145–1156. doi:10.1016/0278-4343(92)90075-u

CrossRef Full Text | Google Scholar

Gentz, T., Damm, E., Schneider von Deimling, J., Mau, S., McGinnis, D. F., and Schlüter, M. (2014). A Water Column Study of Methane Around Gas Flares Located at the West Spitsbergen continental Margin. Continental Shelf Res. 72 (1 January 2014), 107–118. doi:10.1016/j.csr.2013.07.013

CrossRef Full Text | Google Scholar

Ginsburg, G. D., and Soloviev, V. A. (1998). Submarine Gas Hydrates. Russia/Statoil, Stavanger, Norway: VNII Okeanologia, St. Petersburg, 216

Gupta, H. R., and Mehra, R. (2013). Power Spectrum Estimation Using Welch Method for Various Window Techniques. Int. J. Scientific Res. Eng. Technology (Ijsret). 2 (6), 389–392. ISSN 2278 - 0882

Google Scholar

Hornafius, J. S., Quigley, D., and Luyendyk, B. P. (1999). The World's Most Spectacular Marine Hydrocarbon Seeps (Coal Oil Point, Santa Barbara Channel, California): Quantification of Emissions. J. Geophys. Res. 104, 20703–20711. doi:10.1029/1999jc900148

CrossRef Full Text | Google Scholar

Hovland, M., Judd, A. G., and Burke, R. A. (1993). The Global Flux of Methane From Shallow Submarine Sediments. Chemospher. 26 (1-4), 559–578. doi:10.1016/0045-6535(93)90442-8

CrossRef Full Text | Google Scholar

Inguaggiato, S., and Rizzo, A. (2004). Dissolved Helium Isotope Ratios in Ground-Waters: a New Technique Based on Gas-Water Re-equilibration and its Application to Stromboli Volcanic System. Appl. Geochem. 19, 665–673. doi:10.1016/j.apgeochem.2003.10.009

CrossRef Full Text | Google Scholar

Italiano, F., Bonfanti, P., Ditta, M., Petrini, R., and Slejko, F. (2009). Helium and Carbon Isotopes in the Dissolved Gases of Friuli Region (NE Italy): Geochemical Evidence of CO2 Production and Degassing Over a Seismically Active Area. Chem. Geology. 266, 76–85. doi:10.1016/j.chemgeo.2009.05.022

CrossRef Full Text | Google Scholar

Italiano, F., Yuce, G., Uysal, I. T., Gasparon, M., and Morelli, G. (2014). Insights Into Mantle-type Volatiles Contribution from Dissolved Gases in Artesian Waters of the Great Artesian Basin, Australia. Chem. Geology. 2014, 75–88. doi:10.1016/j.chemgeo.2014.04.013

CrossRef Full Text | Google Scholar

Judd, A. G. (2000). “Geological Sources of Methane,” in Atmos. Methane.. Editors M. A. K. Khalil (Berlin, Heidelberg: Springer, 280–303. doi:10.1007/978-3-662-04145-1_15

CrossRef Full Text | Google Scholar

Judd, A. G. (2003). The Global Importance and Context of Methane Escape from the Seabed. Geo-mar. Lett. 23, 147. doi:10.1007/s00367-003-0136-z

CrossRef Full Text | Google Scholar

Judd, A. G. (2004). Natural Seabed Gas Seeps as Sources of Atmospheric Methane. Env Geol. 46, 988–996. doi:10.1007/s00254-004-1083-3

CrossRef Full Text | Google Scholar

Judd, A. G., Hovland, M., Dimitrov, L. I., García Gil, S., and Jukes, V. (2002). The Geological Methane Budget at Continental Margins and its Influence on Climate Change. Geofluids. 2 (2), 109–126. doi:10.1046/j.1468-8123.2002.00027.x

CrossRef Full Text | Google Scholar

Konstantinou, K. I. (2015). Tornillos Modeled as Self-Oscillations of Fluid Filling a Cavity: Application to the 1992-1993 Activity at Galeras Volcano, Colombia. Phys. Earth Planet. Interiors. 238, 23–33. doi:10.1016/j.pepi.2014.10.014

CrossRef Full Text | Google Scholar

Kvenvolden, K. A. (1999). Potential Effects of Gas Hydrate on Human Welfare. Proc. Natl. Acad. Sci. 96 (7), 3420–3426. doi:10.1073/pnas.96.7.3420

PubMed Abstract | CrossRef Full Text | Google Scholar

Leifer, I., and Boles, J. (2005). Measurement of marine Hydrocarbon Seep Flow through Rock and Unconsolidated Sediment. J. Mar. Petrol. Geol. 22 (4), 551–568. doi:10.1016/j.marpetgeo.2004.10.026

CrossRef Full Text | Google Scholar

Leifer, I. (2010). Characteristics and Scaling of Bubble Plumes from marine Hydrocarbon Seepage in the Coal Oil Point Seep Field. J. Geophys. Res. 115, C11014. doi:10.1029/2009JC005844

CrossRef Full Text | Google Scholar

Leifer, I., and Culling, D. (2010). Formation of Seep Bubble Plumes in the Coal Oil point Field. Geo-mar. Lett. 30 (3–4), 339–353. doi:10.1007/s00367-010-0187-x

CrossRef Full Text | Google Scholar

Leifer, I., and Patro, R. K. (2002). The Bubble Mechanism for Methane Transport From the Shallow Sea Bed to the Surface: A Review and Sensitivity Study. Continental Shelf Res. 22 (16), 2409–2428. doi:10.1016/S0278-4343(02)00065-1

CrossRef Full Text | Google Scholar

Leifer, I., and Tang, D. (2007). The Acoustic Signature of Marine Seep Bubbles. J. Acous. Soc. doi:10.1121/1.2401227

CrossRef Full Text | Google Scholar

Leighton, T. G., and White, P. R. (2011). Quantification of Undersea Gas Leaks From Carbon Capture and Storage Facilities, from Pipelines and From Methane Seeps, by Their Acoustic Emissions. Proc. R. Soc. A. 468, 485–510. doi:10.1098/rspa.2011.0221

CrossRef Full Text | Google Scholar

Li, J., Roche, B., Bull, J. M., White, P. R., Davis, J. W., Deponte, M., et al. (2020). Passive Acoustic Monitoring of a Natural CO2 Seep Site - Implications for Carbon Capture and Storage. Int. J. Greenhouse Gas Control. 93, 102899. doi:10.1016/j.ijggc.2019.102899

CrossRef Full Text | Google Scholar

Li, J., White, P. R., Roche, B., Bull, J. M., Leighton, T. G., Davis, J. W., et al. (2021). Acoustic and Optical Determination of Bubble Size Distributions - Quantification of Seabed Gas Emissions. Int. J. Greenhouse Gas Control. 108, 103313. doi:10.1016/j.ijggc.2021.103313

CrossRef Full Text | Google Scholar

Longo, M., Lazzaro, G., Caruso, C. G., Corbo, A., Sciré Scappuzzo, S., Italiano, F., et al. (2021). “Hydro-acoustic Signals From the Panarea Shallow Hydrothermal Field: New Inferences of a Direct Link With Stromboli, Geological Society, London, Special Publications. From,” in Volcanic Island: From Hazard Assessment to Risk Mitigation. Editors E. Marotta, L. D’Auria, F. Zaniboni, and R. Nave (London, UK: Geological Society, London, Special Publications), 519, SP519–2020. doi:10.1144/SP519-2020-184

CrossRef Full Text | Google Scholar

Lu, N. Q., Prosperetti, A., and Yoon, S. W. (1990). Underwater Noise Emissions from Bubble Clouds. IEEE J. Oceanic Eng. 15, 275–281. doi:10.1109/48.103521

CrossRef Full Text | Google Scholar

Mau, S., Römer, M., Torres, M. E., Bussmann, I., Pape, T., Damm, E., et al. (2017). Widespread Methane Seepage Along the Continental Margin off Svalbard - from Bjørnøya to Kongsfjorden. Sci. Rep. 7, 42997. doi:10.1038/srep42997

PubMed Abstract | CrossRef Full Text | Google Scholar

McGinnis, D. F., Greinert, J., Artemov, Y., Beaubien, S. E., and Wüest, A. (2006). Fate of Rising Methane Bubbles in Stratified Waters: How Much Methane Reaches the Atmosphere?. J. Geophys. Res. 111, C09007. doi:10.1029/2005jc003183

CrossRef Full Text | Google Scholar

Milkov, A. V. (2004). Global Estimates of Hydrate-Bound Gas in Marine Sediments: How Much Is Really Out There?. Earth-Science Rev. 66, 183–197. doi:10.1016/j.earscirev.2003.11.002

CrossRef Full Text | Google Scholar

Milluzzo, V., Cannata, A., Alparone, S., Gambino, S., Hellweg, M., Montalto, P., et al. (2010). Tornillos at Vulcano: Clues to the Dynamics of the Hydrothermal System. J. Volcanology Geothermal Res. 198, 377–393. doi:10.1016/j.jvolgeores.2010.09.022

CrossRef Full Text | Google Scholar

Minnaert, M. (1933). XVI.On Musical Air-Bubbles and the Sounds of Running Water. Lond. Edinb. Dublin Philosophical Mag. J. Sci. 16, 235–248. doi:10.1080/14786443309462277

CrossRef Full Text | Google Scholar

Ohminato, T. (2006). Characteristics and Source Modeling of Broadband Seismic Signals Associated With the Hydrothermal System at Satsuma-Iwojima Volcano, Japan. J. Volcanology Geothermal Res. 158, 467–490. doi:10.1016/j.jvolgeores.2006.08.004

CrossRef Full Text | Google Scholar

Özsoy, E., and Ünlüata, Ü. (1997). Oceanography of the Black Sea: A Review of Some Recent Results. Earth-Science Rev. 42 (4), 231–272. doi:10.1016/s0012-8252(97)81859-4

CrossRef Full Text | Google Scholar

Pape, T., Blumenberg, M., Seifert, R., Bohrmann, G., and Michaelis, W. (2008). Marine Methane Biogeochemistry of the Black Sea: A Review, Links between Geological Processes, Microbial Activities & Evolution of Life: Microbes and Geology. Dordrecht: Springer Netherlands, 281–311.

Peregrine, D. H. (1994). The Acoustic Bubble. By T. G. Leighton. Academic Press, 1994. 613 Pp. £95. ISBN 0-12-441920-8. The Acoustic Bubble. From Peregrine, D. H. Academic Press Pp 613. ISBN 0-12-441920-8. J. Fluid Mech. 272, 407–408. doi:10.1017/s0022112094214519

CrossRef Full Text | Google Scholar

Popa, A. (1993). Liquid and Sediment Inputs of the Danube River Into the Northwestern Black Sea. Mitt. Geol.-Palaont. Inst. Univ. Hamburg. 74, 137–149.

Google Scholar

Römer, M., Sahling, H., dos Santos Ferreira, C., and Bohrmann, G. (2020). Methane Gas Emissions of the Black Sea-Mapping From the Crimean continental Margin to the Kerch Peninsula Slope. Geo-mar Lett. 40, 467–480. doi:10.1007/s00367-019-00611-0

CrossRef Full Text | Google Scholar

Sano, Y., and Wakita, H. (1988). Precise Measurement of Helium Isotopes in Terrestrial Gases. Bcsj. 61, 1153–1157. doi:10.1246/bcsj.61.1153

CrossRef Full Text | Google Scholar

Schmale, O., Greinert, J., and Rehder, G. (2005). Methane Emission From High-Intensity Marine Gas Seeps in the Black Sea into the Atmosphere. Geophys. Res. Lett. 32, a–n. doi:10.1029/2004GL021138

CrossRef Full Text | Google Scholar

Schoell, M. (1988). Multiple Origins of Methane in the Earth. Chem. Geology. 71, 1–10. doi:10.1016/0009-2541(88)90101-5

CrossRef Full Text | Google Scholar

Shindell, D. T., Faluvegi, G., Koch, D. M., Schmidt, G. A., Unger, N., and Bauer, S. E. (2009). Improved Attribution of Climate Forcing to Emissions. Science. 326 (5953), 716–718. doi:10.1126/science.1174760

PubMed Abstract | CrossRef Full Text | Google Scholar

Stadnitskaia, A., Ivanov, M. K., Poludetkina, E. N., Kreulen, R., and van Weering, T. C. E. (2008). Sources of Hydrocarbon Gases in Mud Volcanoes From the Sorokin Trough, NE Black Sea, Based on Molecular and Carbon Isotopic Compositions. Mar. Pet. Geology. 25 (10), 1040–1057. doi:10.1016/j.marpetgeo.2007.08.001

CrossRef Full Text | Google Scholar

Strasberg, M. (1956). Gas Bubbles as Sources of Sound in Liquids. J. Acoust. Soc. America. 28, 20–26. doi:10.1121/1.1908212

CrossRef Full Text | Google Scholar

Tolstoy, M., Vernon, F. L., Orcutt, J. A., and Wyatt, F. K. (2002). Breathing of the Seafloor: Tidal Correlations of Seismicity at Axial Volcano. Geol. 30, 503–506. doi:10.1130/0091-7613(2002)030<0503:botstc>2.0.co;2

CrossRef Full Text | Google Scholar

Vazquez, A., Sanchez, R. M., Salinas-Rodríguez, E., Soria, A., and Manasseh, R. (2005). A Look at Three Measurement Techniques for Bubble Size Determination. Exp. Therm. Fluid Sci. 30, 49–57. doi:10.1016/j.expthermflusci.2005.03.018

CrossRef Full Text | Google Scholar

Walton, A. J., Gunn, M. G., and Reynolds, G. T. (2005). The Quality Factor of Oscillating Bubbles as an Indication of Gas Content With Particular Reference to Methane. IEEE J. Oceanic Eng. 30, 924–926. doi:10.1109/joe.2004.834171

CrossRef Full Text | Google Scholar

Wenz, G. M. (1962). Acoustic Ambient Noise in the Ocean: Spectra and Sources. J. Acoust. Soc. America. 34, 1936–1956. doi:10.1121/1.1909155

CrossRef Full Text | Google Scholar

Whiticar, M. J. (1990). A Geochemical Perspective of Natural Gas and Atmospheric Methane. Org. Geochem. 16 (Issues 1–3), 531–547. doi:10.1016/0146-6380(90)90068-B

CrossRef Full Text | Google Scholar

Wiggins, S. M., Leifer, I., Linke, P., and Hildebrand, J. A. (2015). Long-Term Acoustic Monitoring at North Sea Well Site 22/4b. Mar. Pet. Geology. 68, 776–788. doi:10.1016/j.marpetgeo.2015.02.011

CrossRef Full Text | Google Scholar

Wilcock, W. S. D., Stafford, K. M., Andrew, R. K., and Odom, R. I. (2014). Sounds in the Ocean at 1-100 Hz. Annu. Rev. Mar. Sci. 6, 117–140. doi:10.1146/annurev-marine-121211-172423

PubMed Abstract | CrossRef Full Text | Google Scholar

Wong, H. K., Winguth, C., Panin, N., Dinu, C., Wollschläger, M., Ungureanu, G., et al. (1997). The Danube and Dniepr Fans: Morphostructure and Evolution”. GeoEcoMarina. 2, 77–102. doi:10.5724/gcs.12.32.0701

CrossRef Full Text | Google Scholar