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 CH4 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 2 (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 cagelike 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;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 4 ]/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 4 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 chemicalphysical 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 6 t, 81% of which (52 × 10 6 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). 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.

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: 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:

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 Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 678834 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 shorttime Fourier transform and a threshold algorithm was implemented following Peregrine, 1994 andLi 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.
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; 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.
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).
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 678834 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.
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: is the ratio of gas specific heat, at constant pressure and volume, and ρ 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 4 , using Eq. (1) we can assume γ 1.32, ρ 1,030 kg/m 3 , 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 windeffect 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 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 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.

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 Fox Christopher, 2002, 2004;Tolstoy et al., 2002;Ohminato 2006;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. (2009Italiano et al. ( , 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, 2004and Italiano et al., 2009, for further details). The chemical analyses of O 2 , N 2 , CH 4 and CO 2 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%.

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.
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).

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).

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 dayslong 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.
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 2 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 seaf loor 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 andreferences 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).
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 August 2021 | Volume 9 | Article 678834 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. Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 678834 13 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][98][99][100][101][102][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 4 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.

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.

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 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 longterm 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 bottomcurrents 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.