Variability of a natural hydrocarbon seep and its connection to the ocean surface

Natural hydrocarbon seeps are ubiquitous along continental margins. Despite their significance, we lack a basic understanding of the long-term temporal variability of seep dynamics, including bubble size, rise velocity, composition, and upwelling and entrainment processes. The shortcoming makes it difficult to constrain the global estimates of oil and gas entering the marine environment. Here we report on a multi-method approach based on optical, acoustic, satellite remote sensing, and simulations, to connect the characteristics of a hydrocarbon seep in the Gulf of Mexico to its footprint on the sea surface. Using an in-situ camera, bubble dynamics at the source were measured every 6 h over 153 days and the integrated total hydrocarbon release volume was estimated as 53 m3. The vertical velocity was acoustically measured at 20 m above bed (mab) and found to be approximately 40% less than the dispersed-phase at the source, indicating that the measured values are reflecting the plume continuous-phase flow. Numerical simulations predict that the oily bubbles with diameters larger than 8 mm reach the surface with a small footprint, i.e. forming an oil slick origin, deflection of which with wind and surface current leads to the formation of an oil slick on the surface. Nineteen SAR images are used to estimate the oil seepage rate from GC600 for 2017 giving an average discharge of 14.4 cm3/s.

Natural hydrocarbon seeps have been reported from a broad range of oceanographic settings from the coast to the deep ocean over a wide variety of geological environments that affect the biosphere, the hydrosphere, and the atmosphere 1 . Beside the anthropogenic sources, hydrocarbon seeps are the single most important natural source of oil and methane (CH 4 ) that enter the ocean 2 . Best estimates indicate that 3 to 48 Tg/y of geo-CH 4 is released from marine seeps alone 3,4 . Recent global estimates of crude oil seepage range from 0.2 to 2 Mt with a best estimate of 0.6 Mt that accounts for 47% of the global estimate. The large uncertainties in estimates of hydrocarbons released from offshore marine sources are due to a wide range of factors including size, rise velocity, contamination by surfactants such as gas hydrates, composition, upwelling and entrainment effects, and interaction of gaseous bubbles and oil droplets with the local ocean environment that vary constantly with time.
To quantify the bubble/droplet size and rise velocity and the corresponding emission rate of hydrocarbons from the seafloor, investigators have primarily relied on snapshot acoustic [5][6][7][8][9] or optical [10][11][12][13] measurements. Based on these measurements, significant changes in hydrocarbon seepage and bubble venting have been inferred to occur over intervals of months to years 14 . Sufficiently long time series of seepage characteristics rarely exist to conclusively address the day-to-day variability of bubble size and rise velocity or the spatial changes in vent location over an extended period. Time variation in seep emissions is of particular interest. It implies variability in the local background levels against which pollution from human activities is measured, and is relevant at a global scale when seepage is scaled up to all continental margins 15,16 .
The Green Canyon 600 lease block (GC600) is located in the northern Gulf of Mexico (GoM, Fig. 1a). Satellite image records containing chronically visible surface oil slicks suggest it is one of the most prolific natural hydrocarbon seep sites in the region 17,18 . Geophysical studies affirm these findings and show abundant seep formations, including reservoirs, faults, and migration conduits found in seismic profiles 19,20 . The abundance of seeps has made this region an attractive site for studying different biological and physical aspects of natural www.nature.com/scientificreports/ methane and oil emission [20][21][22][23] . Here, we report on a comprehensive observational and modeling approach based on optics for time-lapse video measurements of bubble dynamics and size distribution, acoustics for detecting plume-induced upwelling flow, numerical modeling to predict the migration path of hydrocarbons in the water column, and satellite remote-sensing to quantify the oil discharge from its footprint on the sea surface. As a result, this paper describes an unprecedented study that explores the dynamics of a submarine seep and its interaction with the surrounding environment from the seafloor to the sea surface.

Measurement program
To identify a natural hydrocarbon seep with a visible oil slick on the surface we surveyed the seepage area (1,600 × 2,200 m 2 ) within the GC600 block using an autonomous underwater vehicle (AUV) equipped with a multibeam sonar (see "Methods"). Rising bubbles through the water column scatter acoustic energy effectively, instigating flare-like structures in hydroacoustic data. Figure 1b demonstrates some of these flares caused by gas bubbles nominally larger than 0.2 mm in diameter, which were utilized to localize their seabed source. The seep  www.nature.com/scientificreports/ flares are superimposed on a 3D view of the bathymetry obtained from the multibeam data; shading corresponds to the backscatter intensities. The bathymetry map demonstrates a salt-supported ridge trending NW-SE with two distinctive areas of high-reflectance seafloor. While the observed hydrocarbon seep locations lack distinct morphology, they are highly correlated with patches of high seafloor reflectivity that may be associated with methane-derived carbonate slabs and crusts 26,27 . Diving with the remotely operated underwater vehicle (ROV) Comanche aboard OSV Ocean Project enabled us to visit some of the hydrocarbon plumes directly. Figure 1c,d respectively display the acoustic image and a photo snapshot of the hydrocarbon seep, Mega Plume, which was selected for further measurements. The water depth at the Mega Plume, 27° 22.1915′ N, 90° 34.4580′ W, is approximately 1,180 m. It is about 10 m away from the vents described by Johansen et al. 11 and Wang et al. 10 ; these seep locations were not active during our visit, suggesting a high spatial variability of seeps in the Mega Plume region.

Hydrocarbon source characteristics
To explore the time-space variability of the seep at the seafloor, we positioned two video time-lapse cameras (VTLCs) 11 close to the seep cluster (Fig. 1d, see "Methods"). A short video made from all the last 10 frames (1/3 s) of each of 568 video bursts recorded by the VTLC-B (see Video 1 in the supplemental document) exhibits dayto-day fine-scale spatial and temporal variability of the source over the campaign period. The records indicate that the seafloor and the seepage were very dynamic on daily to monthly time scales. Initially-from Sept to late Nov 2017-gas bubbles, oil droplets, and/or a mixture of the two were released from the edge of an exposed gas hydrate deposit identified as the 'Clam Cluster' in Fig. 2a The emission from the Clam Cluster ceased completely on Dec 28, 2017, for the rest of the campaign. Visual inspection revealed that while the composition of the bubbles emanating from different vents in the Clam Cluster were changing uniformly, the multiple vents along the Fault Cluster released either purely gaseous, oily bubbles or a combination of both, simultaneously. Another interesting feature to note from the time-lapse video is that the seafloor elevation started to change over the entire period -the evident subsidence may be associated with the breakdown of gas hydrate.
The VTLC capability of autonomously recording videos over a long period comes at the cost of reduced frame rate and illumination compared to laboratory 28 and ROV-based techniques 10 . Therefore, we developed a sequence of novel semi-supervised algorithms based on standard Mathematica ® and ImageJ ® routines to resolve the rise velocity, w b , and size, d, of individual bubbles, details of which can be found in Razaz et al. 29 . Utilizing these algorithms, 1 s of each video burst was processed. Figure 3 demonstrates the VTLC-B data together with the complementary ROV-based particle tracking velocimetry (PTV) results (see "Methods") and the synoptic physical collections of hydrocarbons into an ROV-held funnel (see "Methods", Fig. 4).
At the Clam Cluster, the average bubble rise velocity (  Fig. 3b) and the overbar stands for averaging over time. Also shown are the bubble sizes determined by fitting a log-normal distribution to the bubble size PDFs calculated for each 1 s sequence and bubble volume weighted PDFs. The VTLC bubble size and rise velocity averaged over the last 24 h are also compared to the ROV-based PTV results and are respectively, 5% and 20% lower, but are within the range of the complementary measurements. In general, our findings are consistent with previous works 7,10,13,30 that reported velocities less than 30 cm/s for bubbles released within the hydrate stability zone having d < 26 mm. Figure 3c demonstrates the emission rates, Q VTLC , calculated using the bubble size, rise velocity, and frame rate of the VTLC-B. There is outstanding compliance (2% difference) between the − Q VTLC = 6.5 cm 3 /s measured 2.5 h before the synoptic physical collection of hydrocarbons with an ROV-held funnel (Fig. 4), Q FNL = 6.4 cm 3 /s. Of more importance is the significant difference between the volume of hydrocarbon estimates from long-term VTLC records and snapshot funnel measurements. That is, integrating the Q VTLC values 29 over the entire period of our campaign yields an emission volume of 51 m 3 . Had we used the constant Q FNL instead, the total emission would be overestimated by 63% (83 m 3 ) over the same period. This large divergence signals the uncertainty that might be involved in annual estimates of methane released into the hydrosphere from ROV-based synoptic measurements. The difference between the seepage characteristics of the two clusters is summarized in Table 1 which quantifies the bubble dynamics and hydrocarbon discharge. The differences may be related to a wide range of factors including sediment loading, migration pathways, hydrate formation, permeability, and porosity 31 .
Many processes can be postulated to explain short-term (days to weeks) variations in seepage rate including variations in bubble characteristics (size, shape, and composition), properties of gas-liquid systems (density, viscosity, surface tension, density difference between gas and liquid, oil coating, and formation of methanehydrate) and oceanographic conditions (temperature, pressure, and crossflow velocity) 29   www.nature.com/scientificreports/ gas hydrate deposits and subsurface oil and gas pockets potentially impinge hydrocarbon releases and could influence bubble size and gas:oil ratios.
Plume-induced upwelling velocity at 20 m above the seafloor. We measured the vertical velocity of the plume-induced flow using an acoustic scintillation flow meter (ASFM) 37-39 (see "Methods") over two weeks, Sept 3-17, 2017. In précis, the measurement concept relies upon the scattering field to advect across spatially separated acoustic paths and requires persistence of those structures to obtain useful coherence in the received signals at the receivers. As indicated by visual and acoustic imaging (Fig. 1c), the bubble plume expands with height above the bed by entraining fluid and by dispersive mechanisms; comparing the ASFM, − w ASFM = 12.9 ± 2.9 cm/s, with the VTLC results averaged over the same period suggests the increase in mass flux led to approximately 40% lower upward velocities at 20 mab compared to the source. Given the ASFM configuration at GC600 (see "Methods"), the ASFM was sensitive to turbulence length scales of approximately 0.6 m that lie within the inertial subrange of fine-scale turbulence 40 . One implication of this scale is that the ASFM treats the water (continuous-phase) and bubbles (dispersed-phase) as a mixed turbulent medium advected by a vertical velocity that is weighted towards the center of the acoustic path. Therefore, we hypothesize that the measured values are reflecting the vertical velocity of the plume continuous-phase (upwelling velocity) and do not include the dispersed-phase. It is also noted that there is no evident correlation between the upwelling flow variability and the ambient currents induced by the tide or those occurring at lower frequencies (Fig. 5).   41 to track the trajectory and the evolution of bubble size, mass, and composition as a function of depth, which provides information on the final surfacing characteristics for these bubbles. The module takes into account background stratification, advection by ambient currents (see "Methods", Fig. 6b), dissolution, and heat transfer. The bubbles were assumed to be gaseous coated with oil in our simulations with sizes ranging between 1 and 15 mm (equivalent spherical diameter with a stride of 1 mm) at the source. For further details see "Methods". The simulation spanned the period Oct 1-6, 2017, when a synthetic aperture radar (SAR) image was available for comparison. The simulations imply that depending on the bubble size and currents, bubbles take substantially different migration paths, which subsequently affect the transit time of the bubble in the water column. Figure 6c illustrates the difference between two representative bubble trajectories having sizes d = 3 mm and 15 mm and their orbital path with time in the water column as a result of diurnal tides and inertial oscillations at approximately the same frequency. The average transit time of bubbles with d > 8 mm at the source is approximately 40% of that computed for bubbles with d = 4-8 mm. An effect of the decrease in transit-time is a reduction in the size of the surfacing footprint, evident in Fig. 6d, since there is less time for dispersal by turbulence and currents 42 . Also, comparing the surface outbreak of the bubbles with different sizes with the oil slick delineated from the SAR image captured by the Sentinel-1A satellite on Oct 1, suggests that only larger bubbles contribute to the formation of the oil slick origin (OSO) 18 . The surfaced oil is then elongated as a film that drifts away with surface currents and wind (Fig. 6a) 42 . Here, it is not possible to separate the effects of the prevailing wind and near-surface current on the general shape of the oil slick as both vectors were roughly aligned to the southwest.

Surface footprint of oil seepage.
To estimate the hydrocarbon flux from sources on the seafloor using the oil slick properties on the surface, we used 19 SAR images collected over GC600 and its surrounding region between Feb 2017 and Dec 2017. In these SAR images, the wind condition was favorable for the detection www.nature.com/scientificreports/ of oil slicks. A semi-automated image processing routine, the Texture Classifying Neural Network Algorithm (TCNNA) 44,45 , was employed to delineate the oil slicks shown in Fig. 7a. Based on the TCNNA results, the average surface residence-time is computed to be t R = 8.6 h, and the average oil slick surface area is Ā= 4.45×106 m 2 over the GC600 domain, taking into account the wind and surface current conditions (see "Methods"). Choosing a conservative and uniform oil film thickness of 0.1 µm following the published standards 46 , the oil seepage on the seafloor for the GC600 domain is calculated to be Q SAR =14.4 cm 3 /s. While processing the SAR images, the deflection of the OSO from the vent is approximated by a linear relation with water depth only 47 and for GC600, the deflection distance can be as high as 2.3 km 18 . Considering the geographically dense seeps in this region, this coarse estimation degrades our ability to identify the venting source of oil slicks accurately. Figure 7b exhibits three distinct oil slicks on Oct 1, 2017, the middle of which is what we denote as coming from Mega Plume with a surface offset of approximately 1 km from the source as shown also in Fig. 6d; without extra information from the seabed, it is not possible to distinguish sources that fall within less than 2 km from each other. This outlines the value of independent measurements of the transport of oil/gas bubbles through the water column for identifying the surface footprint. If we consider the SAR image denoted in Fig. 3d as coming from the Mega Plume, given a surface residence time of 8.6 h, the discharge is estimated Q SAR−MP = 2.4 cm 3 /s, which compares favorably to our seep measurement of Q VTLC = 3.0 cm 3 /s for Oct 2, 2017.

Methods
The extensive field survey carried out in 2017 in GC600, GoM, consisted of video observations of seep source physical conditions, acoustic Doppler measurements of near-bottom horizontal flows, acoustic forward scatter measurements of the vertical flow, and acoustic backscatter measurements of the water column and seafloor. Full water column horizontal currents were obtained from the nearest NDBC station 42369. These data sets were used to model the bubble rise through the water column for comparison to sea surface oil slick observations obtained from the Sentinel-1A SAR satellite. Each of these data sets is described below.
Video-based observations. Video time-lapse camera (VTLC). Two VTLCs were positioned at approximately 60 and 90 cm from the seep by the ROV Maxx aboard MSV Ocean Intervention II on Sept 3, 2017 (Fig. 1d). VTLC-A recorded 10 s of data every 3 h during Sept 3-28, 2017, and VTLC-B logged 15 s of data every 6 h over the period of Sept 3, 2017-Feb 2, 2018; a video of the 153-day deployment is provided in supplementary data, Video 1.

ROV-based particle tracking velocimetry (PTV).
A 10-min long video of the seep at the source was recorded using the Ocean ProHD TV camera (1080i/59.4) on the ROV Millennium during the recovery cruise on Feb 2, 2018. To surmount the inherent parallax error associated with a single-camera system, a screen with size scale markings was placed behind the seep. Following the PTV methods described in Wang and Socolofsky 48 , 350 s www.nature.com/scientificreports/ of the video was processed to determine the size and rise velocity of individual bubbles for comparison to the VTLC measurements. Fig. 4 was used to measure the emission rate over a period of about 5 min-the time it took to fill the funnel volume. The funnel volume compartment was drawn with 8 visible sections, each indicating 250 cm 3 . Using the ROV Millennium arm, the funnel was kept on top of the seep about 0.15 mab so that no hydrocarbons could escape the funnel while no disturbance was made to the loosely deposited sediment. It was not possible to repeat the sampling more than once since the compartment was filled with a solid mixture of methane hydrate and oil (Fig. 4).

Synoptic physical collection of hydrocarbons. A custom-made funnel shown in
Acoustic observations. Multibeam echosounder. The primary survey sensor on the AUV aboard the Oceaneering vessel Ocean Project includes a Kongsberg EM2040 200 kHz multibeam sonar. An acoustically aided inertial navigation system consisting of an acoustic Doppler velocity log and an ultra-short baseline (USBL) positioning system on the survey vessel was used for accurate localization of the AUV. On Jun 18-19, 2017, we used the AUV to perform a systematic survey to identify natural hydrocarbon seeps in the seepage area of GC600. Throughout the survey, the AUV maintained an altitude of approximately 40 m with a speed less than 1.5 knots while cruising along 15 primary track lines spaced 100 m apart (Fig. 1b), and 3 tie lines spaced at 900 m. The water column and seafloor acoustic backscatter data were collected using 256 beams over a 140° wide swath (approximately 220 m on the seafloor in equidistance mode) at 4 Hz sampling frequency. High-resolution bottom bathymetry was also obtained from the bottom backscatter data.
Acoustic scintillation flow meter (ASFM). An abundance of reports on the scintillation method to monitor horizontal and vertical flow is available in the acoustic and oceanographic literature [37][38][39]49,50 . The scintillation method in its simplest configuration consists of two transmitters emitting acoustic pulses in succession at a fast repetition rate. Each pulse is then received at two spatially spaced receivers to resolve turbulent fluctuations focused on the Fresnel scale, √ L , where λ is the acoustic wavelength and L is the acoustic path length. A crucial requirement for the scintillation method is that there must be sufficient coherence between signals detected at the two receivers. Furthermore, the acoustic scattering must remain weak for reliable amplitude fluctuations 40 . From the time-lagged spatial coherence pattern, the path-averaged current orthogonal to the transmission path is obtained without requiring any geometric correction.
The ASFM, deployed in GC600, has been used in a variety of other deep-sea environments to monitor Mediterranean flow into the Black Sea and hydrothermal plume dynamics 38,49,50 . As depicted in Fig. 1c, using moorings carefully positioned 96.3 m apart, the transmitter and receiver arrays were deployed such that the acoustic propagation paths intersect the bubble stream at ~ 20 mab. The 100 μs transmit pulses with a central frequency of 307 kHz were emitted sequentially with 25 ms separation at 10 Hz sampling frequency. The instrument logged data in burst sampling mode, 15 min every hour. Application of toroidal transducers with a 10° beam width ensures that the sound is propagated through the bubble stream, regardless of where the stream bends.
Acoustic Doppler current profiler (ADCP). The background flow in the near-bottom region was measured using a 600 kHz TRDI Workhorse ADCP mounted on a bottom tripod. The tripod was positioned 68 m from the Mega Plume seep cluster to avoid acoustic interference with the ASFM or with the bubble stream. The range cell size was set to 2 m with the first bin at 2.64 mab with zero blanking distance. The radial beam velocities were recorded in burst mode over the period of Sept 3-Dec 7, 2017, without any ensemble averaging with 912 s (at 0.76 s sampling interval) of data every hour. The radial beam velocities were transformed into East-North-Up (ENU) coordinates using the weighted least squares approach described by Gilcoto et al. 51 . The magnetic declination at the survey site (0.24° W) obtained from the 'Magnetic Field Calculator' readily available on the NOAA website. The magnetic declination is negligible compared to the accuracy of the ADCP compass (± 2°). The hourly averaged dominant north/south flow is exhibited in Fig. 5.
The primary sources of real-time current profile data for the NOAA-NDBC program in the GoM are the ADCPs suspended from or attached to oil production platforms and drilling rigs. The current profilers may look downward (from near the surface or at mid-depth), upward (from a bottom mount or at mid-depth), or horizontally (from near the surface). The raw data and the decoded data with quality control flags are available on the NDBC website. Here we used the archived data from the NDBC station 42369 located at 27° 12.4000′ N, 90°, 16.9667′ W with 1,371.9 m depth (Fig. 1a). Almost 80 km away from the Mega Plume, this is the nearest NDBC station to our seep site with consistently high-quality data throughout the water column. At this station, three ADCPs positioned at approximate depths of 12 (down-looking), 450 (up-looking), and 450 m (down-looking) covered the range 16-132, 148-420, 470-1,050 m respectively sampling every 20 min; results are summarized in Fig. 6b. Physical properties of seawater. Ancillary physical water column data including conductivity, temperature, and depth (CTD) profiles were gathered from an array of SeaBird SM37SB MicroCATs attached to the ASFM receiver mooring. The MicroCATs logged data every 90 s throughout Sept 4-Feb 3, 2017. Temperature variability recorded at 24 mab is given in Fig. 5.
Simulating bubble dynamics with TAMOC. The composition of bubble contents in TAMOC-SBM was taken from Wang, et al. 10  The funnel experiment suggests that bubbles were coated with oil and an ice shell of methane hydrate. Thus, the dissolution rates accounting for the effect of surfactants were adopted with the same method described in Leonte et al. 52 for the gaseous components 53 . The pseudo-components required to describe the oil can be found in the ADIOS oil library (https ://noaa-orr-erd.githu b.io/ADIOS /) as "LIGHT LOUISIANNA SWEET, BP" and its incorporation method into the TAMOC is adopted from Gros, et al. 54 . Under this assumption, the volume fraction of oil at the source (compressed volumes) reaches to about 6.7%, which is smaller than the 10% limit recommended by Kvenvolden 55 for the GoM seeps. The physical, chemical, and thermodynamic properties of bubbles were updated at each time step of the TAMOC-SBM considering the background hydrographic conditions until the bubbles either reach the surface or the gas is dissolved completely.
The TAMOC model incorporates the background currents and an averaged salinity and temperature profile for the entire depth obtained during our cruises on Jun 29 and Sept 2, 2017. To investigate the effect of variability of currents in time and space on the migration path of the bubbles with different size, we chose 50 points in time (marked with dots in Fig. 6c) to preserve the oscillatory nature of the currents particularly near the surface as much as possible while keeping the computational costs reasonable. At each time-step, bubbles with sizes ranging 1-15 mm were tracked from the source to the surface.
Since wind can affect the location of the oil once on the surface, we also show the wind data from the ERA-Interim Archive provided by the ECMWF 43 at 10 m above the mean surface level on a 0.125° × 0.125° horizontal grid. This data is shown in Fig. 6a. Satellite remote sensing. At natural seeps, the oil that makes it to the sea surface spreads out into a thin layer (< 1 µm) from the oil slick origin and floats downstream with wind and surface currents 42 . SAR satellites can readily detect the oil slicks under favorable wind conditions, with speeds ranging between approximately 2-7 m/s. An oil slick will become progressively less visible and disappear over time 47 because of weathering processes. The maximum time that oil in a slick can be detected by satellite observations is the surface residence-time 42 . A flux for the oil released from the seafloor can be estimated based on the residence time of natural oil slicks 56 . Daneshgar Asl, et al. 17 estimate the daily flux of oil released from a seep using the equation Q SAR = Aζ (24/t R ) , where ζ is the oil film thickness.
In total, we obtained 41 SAR images collected by Sentinel-1A satellite from NASA's Earth Science Data and Information System (ESDIS) project which provides SAR data from a variety of satellites and aircrafts through the archives maintained by the Alaska Satellite Facility. After preliminary screening of SAR images for desirable winds, the TCNNA was applied to 19 of the images to delineate the oil slicks. An empirical linear relationship for the GoM seeps, x = 1.2346z + 796.86 as suggested by MacDonald, et al. 47 , groups OSOs together when comparing among multiple images. This relation describes the lateral displacement between an OSO and its seafloor vent location, x, as a function of depth, z. For a depth of ~ 1,200 m, OSOs are clustered over a 2.3 km radius from the MegaPlume source. To obtain the oil slick area the ArcToolbox's Conversion and Spatial Statistics tools were utilized.
To estimate the surface residence-time a set of numerical simulations was implemented considering the wind and surface current conditions. The wind history was obtained from the Cross-Calibrated Multi-Platform (CCMP) 57 project. Surface currents were obtained from the HYCOM-HYbrid Coordinate Ocean GoM 1/25° (GOM10.04) analysis experiment 32.5 Model 58 . The CCMP winds represent 10 m winds above the sea surface, with 6 h temporal resolution and are on a 0.25° × 0.25° horizontal grid. Other input parameters of the surface oil drift model are listed in Table 2. The simulations were carried out by increasing the hindcast interval in reverse time order from the image collection time in order to obtain the closest resemblance between the simulated oil pathways and the length and shape of the oil slicks observed in SAR images. The surface residence-time of the oil is estimated by using the hindcast interval that best reproduces the length and shape of the corresponding oil slick. To compare the final distribution of oil particles from the hindcast simulations against the observed oil slicks in the SAR images, we tracked the directional vectors connecting the oil slick origin location and the farthest point of each oil slick both in the model and observations.