A custom-built bubble release system (BRS) was used to release controlled amounts of CO₂ gas into the water column (Scoulding et al., 2020). A gas cylinder containing 6 kg of compressed CO₂ gas was housed in an aluminum frame (140 cm long and 120 cm tall) at a 15° upward angle to help maintain flow. The gas cylinder was connected to a dual gauge regulator and a 1/8″ screw valve to control the flow of CO₂ gas through to a PVC bubbler containing five equidistant 0.53 mm holes. Two relevant flow rates (determined by Scoulding et al., 2020) were tested, high flow (2 L min⁻¹) and low flow (0.2 L min⁻¹).

The survey was conducted on the 4th of February 2020 off of Kingston Beach in southern Tasmania. The BRS was deployed on a gently sloping sandy seafloor in a water depth of 15 m. MBES WCD was acquired from CSIRO’s 7.5 m vessel South Cape using a pole-mounted Kongsberg EM2040C 1.3° (across-track per beam) x 1.3° (along-track per beam) system. Data were acquired at 300 kHz in continuous wave, single swath, equidistant beam spacing mode (256 beams), using a “very short” pulse length of 25 μs and with 130° swath coverage (65° port and 65° starboard). WCD were recorded as separate *.wcd files with a time-varied gain (TVG) on acquisition of 30 log(r) + 20 dB offset, where r is the range in m. Lines were run at increasing distances from the BRS from both E-W and N-S directions with approximately 20 lines run per 10° distance from the BRS (e.g., 20 passes collected at angles 0–9°, 20 passes at 10–19°, up to 60°). Three speeds were tested (1 m/s, 2 m/s and 3 m/s) for each of the two bubble flow rates (0.2 L min⁻¹ and 2 L min⁻¹) (Figure 2A). A sound velocity cast was taken using a Valeport Monitor sound velocity profiler before data acquisition. The profile revealed a slight thermocline at 10 m depth (Figure 2B). Two additional casts were performed throughout the survey to confirm stable conditions throughout the water column.

Multibeam water column data are recorded natively (in the notations used in this paper) as a series of transmitted pulses (from the transmitter array) which are recorded from the receiver array in the across track direction and referred to as “pings” (Figure 1A). Each ping consists of a certain number of formed “beams” (e.g., 256 for the Kongsberg EM2040C in equidistant mode) which are further broken down into time “samples” (Figure 1A), where the signal amplitude of each sample is recorded as the return echo strength at that range. Throughout this work, “target” refers to a sample that is considered to be a real reflector within the water column (i.e., not considered noise or artifact) and “feature” refers to a coherent object (e.g., gas seep, fish school) that is made up of targets.

The sample-beam data can be geometrically processed to display each ping as a projection in depth/across track coordinates, as in Figures 1B–D. This is the typical multibeam “fan” that converts the sample/beam data into real-world XYZ coordinates. In this across-track view, the sidelobe interference manifests as a “semi-circle of radius equal to the minimum slant range (MSR) to the seafloor” (Clarke, 2006) where the sidelobe noise is due to the impact of strong targets on neighboring beams in the receiver array, primarily from the strong specular seafloor return but also due to strong reflectors within the water column. The strength of this sidelobe noise at any given range is the result of the strength of the seafloor (or other strong reflector) return combined with the beam pattern of the receiver array. For this paper, the data above the MSR which is assumed to not be impacted by receiver array sidelobe interference due to seafloor echoes is referred to as MSR_A, and the data beyond the MSR which is generally confounded by this type of sidelobe interference is referred to as MSR_B (Figure 1B). In addition to sidelobe interference, MBES WCD can be contaminated with a variety of other noise sources, such as acoustic interference (impulse noise), ship noise, weather or turbulence from the surface causing bubbles to sweep under the transducer. Examples of various types of noise and data characteristics of WCD are shown in Figure 1C and a comprehensive list of noise sources and how they appear in WCD can be found in Clarke (2006). Beam pointing angle, as reported here, is the angle of a specific beam from nadir (Figure 1D).

The sequence of a series of pings along the ship track can be combined into a 3D volume (triangular volume in Figure 1A). The insonified water volume decreases with increasing beam angle, as described by Urban et al. (2016) and Schimel et al. (2020). This presents an additional challenge when using the full extent of the MBES WCD insofar as features near the seafloor in the outer beams are truncated by the angle of the beam. In Urban et al. (2016), the insonified water volume is further reduced by excluding the volume of water beyond the MSR (Figure 3 in Urban et al., 2016) (i.e., limited to only the data within MSR_A (Figure 1B). However, this study’s goal is to extend the useable insonified water volume to include data within the full angular aperture of the swath (MSR_A + MSR_B). For targets such as gas seeps that present throughout the water column in shallow water, this limits the detection potential in the outer beams.

The raw values recorded for water column amplitudes may vary depending on manufacturer specifications and the corrections that are applied on acquisition. For the Kongsberg data used in this study, the amplitude values within the WCD datagrams (A _WC ) are recorded in 0.5 dB resolution based on Eq. 1 in Gurshin et al. (2009). The raw Kongsberg EM2040 WCD amplitudes were processed using custom Matlab scripts (MATLAB 2019b) built on the open-source CoFFee toolbox (CoFFee, 2020). In this method, water column targets were extracted automatically using sequential filtering steps, with each step further refining the extraction based on known noise patterns of the WCD to isolate real targets. “Targets” in this context, refer to individual WCD samples potentially corresponding to real features within the water column.

By running the filtering steps, the original *.wcd data files are reduced from a file size of 300 MB (approximately 500 pings collected over 60 s) to XYZ text files of 1–6 MB. The output file includes information for each sample target such as raw amplitude, S _vu , beam number, beam pointing angle, sample range, nadir bottom range, ping number and whether the sample is within the MSR_B. In addition to the XYZ point cloud, a quality control image is generated, showing the targets retained by ping after the filtering process. Figure 6 shows examples of these images for a gas seep detected at (a) nadir and (b) at 51°. In these images, the dark blue points are targets identified within the MSR_A region and orange points indicate targets identified within MSR_B. Even when fully within MSR_B (Figure 6B) the gas seep is visible as the prominent feature before classification. The reduced vertical extent of the seep at the higher angle is due to the reduced sampling volume as angle from nadir increases (i.e., the bubbles rise out of the detection range of the beams).

Individual gas seep features are automatically extracted from the raw WCD amplitudes in the filtering and classification steps. An example of one classified point cloud cluster is shown in Figure 7, where the left image illustrates the raw water column fan in 3-D and the right image shows the extracted point cloud overlaid in black. The same survey line is presented as Figure 6B, where the seep was located fully within the MSR_B region with an angle to the base of the gas seep at 51°.

Examples of gas seeps extracted as 3-D point clouds at various angles from 0° to 56° are shown in Figure 8. All of these were from the low flow rate (0.2 L min⁻¹) release with the vessel travelling at 3 m/s. The decreasing vertical detectable extent of the gas seep with increasing angle is due to the decreasing sampling volume near the edges of the water column fan. Gas seep clusters were composed of numbers of point targets between 110 and 45,292, with an average of 5,842 targets per gas seep feature. Numbers of targets within the gas seep clusters decreased with increasing speed and beam angle, due to the increased distance between pings at higher speeds and the decreased sampling volume at higher angles.

Each of the resulting gas seep point clouds was checked for accuracy against the true seep location and original water column fan images. A total of 124 survey lines with gas seeps with a verified location of the BRS were analyzed across a range of angles from 0 to 60°, along with 63 additional survey lines of data either from line turns or passes where there was no flow from the BRS to test for false positives or false negatives. Out of these 187 survey lines, 178 were correctly classified (either as true gas seep or true negative), indicating ∼95% detection success. The total number of correctly and incorrectly identified seeps binned by base angle is shown in Table 1. At base angles of greater than 50°, 19 seeps were correctly identified out of the 23 lines tested. For angles of 40–49°, 20 out of 21 were correctly identified. At angles less than 40°, all seeps were correctly identified.

For the “control” survey lines with no gas seeps a total of four false positives were detected. At least one false positive was due to anomalous loss of data in one half of the swath for a few pings which allowed vessel noise to bypass the threshold filters. For the survey lines containing the target gas seep (n = 124) the single gas seep feature was classified as two separate clusters five times. In these cases, the targets were divided between the MSR_A and MSR_B portions of the data (Figure 8D). This was likely due to the high noise at the interchange between these regions, which can cause a gap in the points greater than the 1.0 m clustering threshold. However, in these cases, the gas seep was identified correctly, with each section classified as a separate gas seep cluster. Of the five false negatives (where a gas seep should have been detected) only one occurred at angles less than 55° (occurring at 44°). Four of the five false negatives were made at the highest speed tested of 3 m/s. Only one false negative was identified for the high flow seep, which had a base angle of 58°.

The constant multiplier, k, represents the number of standard deviations from the mean from which data are identified as real targets above the background noise. Throughout the development of the method many values of k were tested. We present here two different values (one above and one below 2) to demonstrate the impact of this value on detection capabilities. These data are presented in Table 2. By decreasing k from 2 to 1.8, one additional gas seep at angles greater than 50° was correctly identified. However, the number of false positive detections increased significantly due to the lower threshold values retaining additional potential targets. Increasing k from 2 to 2.2 resulted in more true seeps (14 in total) being missed due to the more aggressive threshold values. The method presented using k = 2 was determined to provide the optimal results in terms of high confidence of detections with low chances of false positives.

Alterations to the workflow were also tested by removing certain steps of the filtering process. The filtering steps were run without first removing the samples from and below the seafloor echoes (step 2) (Alternative 3, Table 2). This resulted in a reduction in success rate from 95% to 72.0% with 32 of the 124 seeps no longer detected. Alternative 4 tested the impact of removing the neighborhood filtering step (step 5). In this case, the success rate dropped from 95% to 87.9%. Although all the seeps correctly classified by the original method were the same for this scenario, and three additional seeps were correctly detected, the removal of the neighborhood filter resulted in significantly more false positives as the noise that the filter would have removed was incorrectly classified as seep related.

Figure 9 shows a comparison of the primary and alternative methods with the results obtained by each method reported as a percentage of the correctly identified seeps depending on beam angle. This is the ratio of the positively identified seeps divided by the total number of verified gas seeps within the dataset. This does not include the control survey lines or false positives. The results were binned by angle in 5° increments between 0 and 60° and plotted at the midpoint of each bin. This figure shows that the primary method has 100% success at identifying the seeps up to 40°, with the reliability dropping to ∼50% only past 55°. The alternate methods (dashed lines) result in less consistent detections, even at angles <40°. In particular, the alternative method that did not remove the seafloor echoes prior to calculating the threshold statistics decreases significantly past the transition to MSR_B, indicating that this step is critical for successfully determining an appropriate threshold for each sample range.

Figure 10 illustrates the false positive detections for each processing scenario. It highlights the marked increase in false-positive detections for reducing k from 2 to 1.8, skipping the removal of seafloor detection step, or skipping the neighborhood filtering step. Combined with the correct percentages displayed in Figure 9, the primary method using k = 2 provides optimal results while minimizing false detections.

Changing the thresholds of the classification workflow only altered the results slightly. Distance thresholds between 0.5 and 1 m made no difference to the overall results, with the only difference being that a 0.5 m threshold resulted in more seeps being split into two separate clusters at the MSR_A/MSR_B transition. Increasing the threshold to 2 m also does not impact the overall detections, but more possibly unrelated targets are included in the classified seeps. Most gas seeps within the outer beams (>50°) contained between 100 and 600 points up to speeds of 3 m/s. Changing the point threshold up to 500 points decreased the gas seep detection rates in the outer beams (17 false negatives, overall success of 90% with no seeps detected beyond 54°), and 100 was determined to be the optimal count for retaining the most gas seeps without resulting in too many small false clusters.

The background noise results (in S _vu ) for samples within MSR_A and MSR_B were averaged for each speed using 75 pings of “empty” WCD (no gas seep present) for each file to give an overall indication of background noise levels. Within MSR_A the background levels were −66.9, −66.2, and −54.5 dB for 1, 2, and 3 m/s, respectively. Within MSR_B the background noise levels were −37.4, −36.9 and −34.2 dB for 1, 2 and 3 m/s, respectively. Background noise levels increased at higher speeds, likely due to an increase in vessel induced noise at higher speeds. This is more notable in the MSR_A region, which is expected as the sidelobe noise dominates most other noise sources in the MSR_B region.

Histograms of the S _vu and raw sample amplitudes within each region of the WCD fan are shown in Figure 11. This highlights the distinct higher background noise levels within the data MSR_B. The S _vu of the gas seep-related targets are within the upper limits of the MSR_B background levels, indicating that targets can be extracted using the filtering methods even close to the background noise levels.

The maximum, minimum and mean S _vu values were determined for each survey line for a horizontal slice of the extracted gas seep targets between 12 and 13 m depth. Figure 12A shows the mean S _vu of targets for the high flow and low flow rate gas seeps across all angles and speeds and whether the samples were detected within the MSR_A or MSR_B regions. The mean S _vu shows consistent values across all speeds and angles. There is a slight increase for the values within the MSR_B region, which is due to some lower amplitude samples being retained as targets within the MSR_A region, so the averaged values become biased at the lower angles. Figure 12B shows the maximum S _vu of the same target sets as Figure 12A. These data also indicate that beam angle does not significantly impact amplitude, although there is a slight decrease above 55°. The average of the maximum samples was −16.1 dB for the high flow (ranging between −11.8 and −21.2) and −21.5 dB for the low flow (ranging between −18.0 and −27.8). These values indicate a ∼6–7 dB detected difference between the high and low flow rates. The average of the minimum samples was −31.4 dB for the high flow (ranging between −33.1 and −28.4) and −34.0 dB for the low flow (ranging between −37.2 and −28.3). These values indicate that the detection limit for these target using this method is ∼3–5 dB above the background sidelobe noise (samples within MSR_B).

The amplitudes of the seep targets that were detected were generally well above the sidelobe noise and had consistent amplitudes across all angles, indicating that the chosen target could be reliably detected. These results show it is possible to detect seeps at a rate of 0.2 L min⁻¹ out to 55° and speeds up to 3 m/s.

Previous studies have shown that using the average signal (Schimel et al., 2020) and variance (Liu et al., 2019) of WCD samples along a common range arc can effectively suppress sidelobe noise. The method presented by Liu et al. (2019) reduced the sidelobe interference, but it did not demonstrate the full extent to which targets can be identified across the entire multibeam swath. Schimel et al. (2020) provides the most promising method to date for fully normalizing water column returns across the swath. The methods presented here build on this concept by using the variance of each sample range as a threshold to remove background noise rather than normalize it, to extract coherent features to the full extent of the multibeam swath. Other methods were also tested, such as using quantiles to determine the threshold values at each sample range rather than the mean and variance. The 95th and 97.5 quantiles were tested (substituted for step 3a), and it was found that the 95th quantile retained too much noise, while the 97.5 missed a significant number of real seeps. Therefore, this technique was not pursued further.

Some multibeam systems, such as the Kongsberg EM710, use sectors to try to minimize sidelobe interference on acquisition, which may impact how effective this method is for filtering sidelobe noise for those systems. However, it may be possible to modify the algorithm to determine thresholds based on sample range and sector and apply the thresholds separately to each sector. Zhao et al. (2017) used a similar method wherein a threshold was determined by μ+ k σ, which was calculated for the amplitude values within a whole sector rather than by sample range. This method effectively removed background noise for automated detection of gas plumes, but was only applied to data within MSR_A.

Other methods for suppressing sidelobes at the beamforming stage have been proposed (Lønmo et al., 2020), and some sonar systems specifically designed to limit sidelobe levels, such as the Simard ME70, have been developed (Trenkel et al., 2008). Such methods may lower sidelobe levels to enable additional detections within the MSR_B region. However, many MBES systems still employ a Mills Cross configuration with sidelobe suppression limited to Dolph-Chebyshev array shading on beamforming (Lurton, 2016). Thus, until suitable and affordable technological solutions can resolve the sidelobe interference issue, post-processing methods must be developed to effectively filter the sidelobe signal to take advantage of the MBES WCD from many of the MBES systems in use today. The method described here provides an easily implemented approach to optimizing data analysis within conventional Mills Cross beam forming scenarios.

The value of k in this workflow was shown to have an impact on the results. The optimal value was determined through experimentation, but comparable conclusions have been found from similar analyses (Zhao et al., 2017; Liu et al., 2019). Zhao et al. (2017) determined optimal values of k of 1.2 and 1.5 for different datasets and noted that this value may depend on depth. Liu et al. (2019) set their constant multiplier parameter (b) to 1.8, noting that increasing or decreasing this value resulted in either too many targets being removed or too many false positives, respectively. Adjustments to k may be necessary, depending on the depth of the survey area, frequency of the MBES system, or backscattering strength of the target of interest. If implemented within commercial software, the ability to adjust this parameter would improve the flexibility of optimizing the method for various datasets. As a general rule, a lower k value may be sufficient for datasets with overall low noise (e.g., calm survey conditions, low ship noise) but a higher k may be required for highly contaminated data. A lower k would result in detecting targets closer to background noise levels but may result in more false positives. Knowledge of the data quality and target of interest should help guide the choice for a suitable value of k.

In this study, the classification was tailored to shallow gas seeps, but this could be modified for any other target of interest, such as fish schools, vegetation, or wrecks, provided the amplitude of the target is at least above the MSR_B background noise levels and has a defining morphological characteristic. Gas seeps have a characteristic shape referred to in the acoustic signature as a “flare” (Judd et al., 1997), which is generally a narrow feature extending from the seafloor. Flare shape is influenced by ocean currents and the rise velocity of the bubbles, and the stability of the shape allows for these characteristics to be used as a signature to distinguish seeps from other biological targets (Judd et al., 1997; Greinert, 2008; Ostrovsky, 2009; Schneider von Deimling and Papenberg, 2012). By understanding the characteristics of different water column targets, we can implement criteria to guide classification rules. It must be noted, however, that these assumed characteristics may change due to dynamic oceanographic conditions (e.g., water currents, thermoclines) which impact the horizontal path of rising bubbles streams (McCartney and Bary, 1965; Schneider von Deimling et al., 2010; Wilson et al., 2015), and therefore may need to be modified depending on knowledge of the local oceanography. Additionally, full ray tracing of the water column samples is traditionally not applied. In locations with distinct thermo- or pycnoclines, this may impact the relative sample range of targets, the distribution of noise along a sample range, and the morphological appearance of cohesive features.

The high resolution of MBES allows repeated samples over a feature, which can then be used for the clustering algorithm to identify coherent features and classify them based on morphological characteristics. Due to the shallow water, narrow beam footprint and high ping rate of the system used here, coherent water column features were made up of many “targets” retained from the filtering steps, which allowed for a relatively high number (100) to be used for the clustering algorithms. For deep water features, sonar systems operating at lower frequencies with lower beam resolution or depending on differences in vessel speed, the number of targets that contribute to a feature may be reduced and require different classification thresholds. Reduced numbers of targets may worsen detection limits.