Acoustic Predictors of Active Fluid Expulsion From a Hydrothermal Vent Field, Offshore Taupō Volcanic Zone, New Zealand

Understanding fluid expulsion is key to estimating gas exchange volumes between the seafloor, ocean, and atmosphere; for locating key ecosystems; and geohazard modelling. Locating active seafloor fluid expulsion typically requires acoustic backscatter data. Areas of very-high seafloor backscatter, or “hardgrounds,” are often used as first-pass indicators of potential fluid expulsion. However, varying and inconsistent spatial relationships between active fluid expulsion and hardgrounds means a direct link remains unclear. Here, we investigate the links between water-column acoustic flares to seafloor backscatter and bathymetric metrics generated from two calibrated multibeam echosounders. Our site, the Calypso hydrothermal vent field (HVF) in the Bay of Plenty, Aotearoa/New Zealand, has an extensive catalogue of vents and seeps in <250 m water depth. We demonstrate a method to quantitatively link active fluid expulsion (flares) with seafloor characteristics. This allows us to develop predictive spatial models of active fluid expulsion. We explore whether data from a low (30 kHz), high (200 kHz), or combined frequency model increases predictive accuracy of expulsion locations. This research investigates the role of hardgrounds or surrounding sediment cover on the accuracy of predictive models. Our models link active fluid expulsion to specific seafloor characteristics. A combined model using both the 30 and 200 kHz mosaics produced the best results (predictive accuracy: 0.75; Kappa: 0.65). This model performed better than the same model using individual frequency mosaics as input. Model results reveal active fluid expulsion is not typically associated with the extensive, embedded hardgrounds of the Calypso HVF, with minimal fluid expulsion. Unconsolidated sediment around the perimeter of and between hardgrounds were more active fluid expulsion sites. Fluids exploit permeable pathways up to the seafloor, modifying and refashioning the seafloor. Once a conduit self-seals, fluid will migrate to a more permeable pathway, thus reducing a one-to-one link between activity and hardgrounds. Being able to remotely predict active and inactive regions of fluid expulsion will prove a useful tool in rapidly identifying seeps in legacy datasets, as well as textural metrics that will aid in locating nascent, senescent, or extinct seeps when a survey is underway.

A direct link between active fluid expulsion and acoustic hardgrounds however remains unclear. Varying and inconsistent spatial relationships between observable fluid expulsion, hardgrounds, and the surrounding seafloor have been noted previously (Naudts et al., 2008;Naudts et al., 2010;Thorsnes et al., 2019). This unclear relationship reflects the complexity of hardground genesis, gas migration, fluid expulsion processes, and spatio-temporal differences in expulsion and embedment. Active fluid expulsion can be both short-(decades) or long-lived (thousands of years) (Judd and Hovland, 2009;Levin et al., 2016). Fluids exploit permeable pathways upward to the seafloor, with preliminary seafloor expression (either liquid or gas bubble expulsion) often occurring in areas of soft unconsolidated sediment due to an unrestricted pathway (Levin, 2005;Judd and Hovland, 2009). Hardground formation begins within the subsurface and, given enough time, can become embedded on the seafloor. Once the conduit closes or ceases activity, hardgrounds will remain at the surface or in the subsurface, with surface hardgrounds thenceforth observable in seafloor acoustic backscatter (Judd and Hovland, 2009;Levin et al., 2016). Formalizing the link between gas flares and hardgrounds has thus remained elusive due to spatiotemporal differences between active expulsion and hardground formation, which can post-date expulsion activity.
Here, we present a study from the active Calypso hydrothermal vent field (HVF), in the offshore Taupō Volcanic Zone, Aotearoa/New Zealand to address whether there are measurable links between gas flares, as a proxy for active fluid expulsion, and seafloor characteristics, in particular hardgrounds. Calypso HVF represents a rich laboratory of active fluid expulsion, in gas and liquid form, with a mix of vents, seeps, and hardground features, with a formal catalogue of active fluid expulsion spanning over five decades (Duncan and Pantin, 1969;Glasby, 1971;Sarano et al., 1989;Stoffers et al., 1999;Botz et al., 2002;Mitchell et al., 2004;Lamarche and Barnes, 2005;Schwarz-Schampera et al., 2007;Hocking et al., 2010;Wysoczanski et al., 2015).
Our dataset encompasses bathymetry and backscatter from two calibrated multibeam echosounders (MBES) operating simultaneously at 200 and 30 kHz (nominal) frequencies.
Seafloor backscatter reflectivity and penetration depends, in part, on the signal frequency or aperture of the chosen MBES, as well as the acoustic impedance characteristics of the seafloor (a function of the medium's density and consequent sound speed) (Jackson et al., 1986;Guillon and Lurton, 2001;Huff, 2008). The acoustic response, or backscatter intensity returned, is dependent on interactions between wavelength, seafloor roughness, seafloor scattering regime (a function of substrate particle size), and volume scattering effects (a function of signal penetration), once radiometric and geometric compensations are applied (Urick, 1956;Medwin and Clay, 1997;Hughes Clarke, 2015). The seafloor backscatter mosaic produced may therefore not directly correlate to the feature of interest if the MBES penetration potential is too high or low. Specific to hardgrounds, the MBES frequency may be unable to reveal hardgrounds embedded beneath the surface unless it can penetrate to the appropriate depths or may not accurately capture surface expressions if the frequency is too low.
Here, we use a classification approach, involving image segmentation and random forest modelling, to: 1 Quantitatively link water-column acoustic flares (here considered as a proxy for active fluid expulsion) with specific seafloor characteristics to determine potential links to hardgrounds in this region; 2 Assess whether these specific seafloor characteristics can be used to develop predictive spatial models of active fluid expulsion; 3 Evaluate if a low (nominally 30 kHz), high (200 kHz), or a combined frequency model increases predictive accuracy of seafloor fluid expulsion; and 4 Investigate the role hardgrounds or surrounding sediment cover may have on the accuracy of a predictive model.

The Calypso Hydrothermal Vent Field
The Calypso HVF is a well-known gas-hydrothermal region of the offshore Taupō Volcanic Zone (TVZ), in the center of the Bay of Plenty, Aotearoa/New Zealand. It lies ∼40 km offshore Te Ikaa-Māui/North Island and 7 km SW of the Whakaari/White Island active volcano ( Figure 1A) (Sarano et al., 1989;Stoffers et al., 1999;Botz et al., 2002). The Calypso HVF is part of the Central Volcanic Region, that represents the volcanically active back-arc region of the Hikurangi subduction system located ca. 300 km to the east. The extensional tectonic setting associated with the back-arc has resulted in a dense network of NE-SW trending normal faults (Wright, 1992), expressed at the seafloor by a pervasive system of scarps (Lamarche and Barnes, 2005). The TVZ extends offshore from the Bay of Plenty coast into the NE-trending depression of the Whakatāne Graben. Active fluid expulsion was formally catalogued at Calypso HVF from water-column acoustic scattering and visible sea surface bubbles (Duncan and Pantin, 1969;Glasby, 1971) and then from video observations (Sarano et al., 1989;Hocking et al., 2010). These fluid expressions are evidence of both diffuse and discrete gas emission, predominantly CH 4 and CO 2 , from submarine geothermal activity (Stoffers et al., 1999;Botz et al., 2002). Calypso HVF fluid expulsion is suggested to be structurally controlled by the pervasive NE-SW trending fault system, with observed fluid expulsion ephemerality inferred to be related to seismic activity (Pantin and Wright, 1994). Stoffers et al. (1999) observed numerous active fluid expulsion areas aligned with fault scarps. In addition to the North Calypso HVF (including the original Calypso site from Sarano et al. (1989)) a second site, now known as the South Calypso HVF, was also detected with two discrete zones of fluid expulsion (Stoffers et al., 1999). Subsequent voyages produced bathymetric maps of the area (25 m grid resolution) using a Kongsberg EM302 (30 kHz) MBES and generated a good knowledge of the tectonic fabric along with surface, water-column, and seafloor evidence for fluid expulsion and mineral precipitation (Mitchell et al., 2004;Schwarz-Schampera et al., 2007;Wysoczanski et al., 2015).

Multibeam Data Acquisition and Processing
Data were acquired during the 2018 RV Tangaroa voyage TAN 1806, a component of the international project QUOI: "Quantitative Ocean-column Imaging using hydroacoustic sources" (Lamarche et al., 2019). The TAN1806 survey employed Kongsberg EM 2040 (200 kHz) and EM302 (30 kHz) MBES to map the Calypso HVF, operated synchronously. Both MBES systems are hull-mounted on RV Tangaroa. MBES compensation occurred at the start of the voyage at the NIWA test-patch sites in Palliser Bay and Palliser Bank, Aotearoa/New Zealand. Calibration was undertaken over multiple depths on seafloor zones with flat, homogenous sedimentology of silt or sand Lucieer and Lamarche, 2011;Ladroit et al., 2017). Lines were run for each MBES setting (frequency/pulse) anticipated. Backscatter correction curves were generated to compensate for beam pattern variations by: 1) removing all Kongsberg backscatter corrections, 2) modeling the backscatter angular variations, and 3) modeling the EM302 transducer beam pattern.
Sound velocity profiles (SVP) were taken at 12-h intervals during operations, 3-h intervals during compensation, at the start of new mapping areas, or when noticeable artefacts appeared in the data. MBES data were acquired via Kongsberg Seafloor Information System (SIS) software (v4.2.1) and cleaned in Caris HIPS/SIPS (v10.4). Vessel position, motion, SVP, and tidal corrections were merged with raw range/angle data to generate geolocated soundings. Manual cleaning removed remaining invalid soundings and data spikes. A bathymetric grid at 5 m resolution was produced in NZTM projection (EPSG: 2193) ( Figure 1B).

Active Fluid Expulsion Identification
Active fluid expulsion was identified in the EM302 MBES (30 kHz) water-column backscatter data using the QPS FM Mid-Water (v7.7.4) Feature Detect tool. Acoustic anomalies detected were analyzed in QPS Fledermaus (v7.7.4) ( Figure 2). Acoustic anomalies that were attached to the seafloor, "flares," were classed as potential fluid expulsion features or "seeps," indicative of gas bubbles emitting from the seafloor. Point clusters of potential seeps were exported to MATLAB and filtered based on a maximum distance between points (40 m), minimum number of points (200), and aspect ratio (1) of each cluster. A seafloor or "seep base" location was determined for each seep; locations were exported as a point shapefile.

Seafloor Backscatter
MBES seafloor backscatter data were imported into SonarScope software (IFREMER, France;Augustin (2016)). Bathymetric sounding validity flags from Caris HIPS/SIPS were applied to corresponding backscatter samples. Kongsberg backscatter corrections applied at time of recording were removed while compensation for local backscatter angular variation and a transmit beam pattern were applied to ensure precise data reduction (Lurton and Lamarche, 2015). Very high overlap of MBES lines during the survey allowed only beams outside nadir (-20°/20°) to be considered, with an average backscatter (decibel, dB) value calculated from multiple samples for each grid point. Backscatter mosaics of both frequencies (200 and 30 kHz) were produced at 2 m resolution ( Figure 3).
An averaging filter (3 × 3 moving window) was passed over the backscatter mosaic to interpolate small pixel-wide gaps from missing data. Both mosaics were clipped from −30 dB to −5 dB to remove large areas of no data and mosaic artefacts (outliers). Both mosaics were classified individually, following Jenks Natural Breaks method (Jenks, 1967) to highlight comparable and dissimilar backscatter classes within each mosaic ( Figure 4). Jenks method sorts values into natural classes by reducing the variance within classes while maximizing variance between classes (Jenks, 1967).

Segmentation of Seafloor Acoustic Data
Backscatter mosaics were segmented using Geographic Object Based Image Analysis (GEOBIA) to identify homogenous regions of seafloor for input into a random forest model. This was achieved using a multiresolution segmentation algorithm in eCognition software (v9.4). Three mosaic configurations were segmented, with bathymetry incorporated in each: 1) 200 kHz backscatter only, 2) 30 kHz backscatter only, and 3) equally weighted EM2040 and EM302 backscatter. Three regions were segmented: the full Calypso HVF dataset, the North Calypso HVF subset, and the South Calypso HVF subset. The individual Calypso HVF regions were analyzed separately to: 1) account for the large differences in acoustic response observed between the north and south ( Figure 3); 2) remove the influence of poor weather acquisition artefacts in the North Calypso HVF 30 kHz mosaic ( Figure 3B); and 3) account for the historic classification of the Calypso HVF-the North Calypso HVF region includes the data described by Duncan and Pantin (1969), while the South Calypso HVF includes the South flare areas, described by Stoffers et al. (1999) and Botz et al. (2002).
Each mosaic configuration and region were segmented using seven different scale parameters (SPs): 9, 14, 26, 46, 50, 75, and 100. The Estimation of Scale Parameter tool (ESP2) (Dragut et al., 2010) was used to determine the appropriate SP for analysis. The SP is a unitless term to describe the minimum size of the segments. SP levels of 9, 14, 26, and 46 were selected automatically using ESP2. Larger, manually selected SPs (50, 75, and 100) were also interrogated to explore differences between auto-and manually selected SPs. We selected segmentation settings that emphasize the influence of color (greyscale color heterogeneity, in the case of backscatter mosaics) rather than shape (shape 0.1, color 0.9). For eCognition, shape and color parameters influence each other (high influence of color equals low influence of shape, and vice versa). Compactness (border length/area) and smoothness of the boundary were equally weighted (0.5) with the aim of extracting backscatter areas with a high color contrast to surrounding seafloor (i.e. "hardgrounds"). Segments along the immediate  Very dark areas (black) in the 30 kHz mosaic show acquisition artefacts from aeration of the transducer due to rough sea conditions. Although original mosaic ranges shown include extreme outliers, color ramps are normalized here to ensure comparable grey tones. Mosaic distributions predominantly lie within −30-−5 dB.
Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 785396 6 perimeter of the survey area were removed to avoid edge effects created from backscatter processing.
Similar to the original Jenks classification of individual mosaics (pixel-based classification), segments were classified into five classes based on mean backscatter intensity values (derived from all pixels that contribute to a segment), following Jenks Natural Breaks method (Jenks, 1967) (objectbased classification). The five classes reflect the mean backscatter values for each segment from lowest to highest values and each class capturing internally homogenous areas of seafloor ( Figure 6).

Textural Analysis of Backscatter Layers
Statistical attributes were generated for each segment, including brightness (combined intensity of both mosaics), maximum difference (backscatter dB difference), and the mean, standard deviation, and skewness of both the 200 and 30 kHz mosaics ( Table 1). These values are calculated from all pixel backscatter dB levels that construct each individual image segment. Each segment produces one value for each measure.
Textural values were also generated for each segment in eCognition using the Haralick grey level co-occurrence matrix (GLCM) (Haralick et al., 1973). Texture metrics are commonly used to identify differences between seafloor substrates or associated habitats and can be used as input layers for seafloor classification models (Lucieer, 2008;Buscombe, 2017). Specific to hardground and fluid expulsion in a shallow hydrothermal system, GLCM texture variables provide high-resolution seafloor characteristics in areas where large topographic changes are limited (i.e. flat seafloor). Here, we used the invariant direction, which is the sum of values for four spatial directions (0°, 45°, 90°, and 135°), for the two mosaics individually and then combined. The GLCM analyses relationships between two pixels at a time, a reference and neighbor pixel, using an offset of 1 pixel (i.e. direct neighbor pixels). All backscatter derivative descriptions are provided in Table 1.
Each segmentation layer (for each mosaic configuration and SP) was exported as a shapefile with all textural and geomorphic values attached to each segment. The segment layers and a point shapefile of flare locations (from Section 3.4) were linked in ArcGIS. Segments that were co-located with flare locations were Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 785396 classed as "active" segments and those without flare locations were classed as "inactive." This process created attribute tables for each segmentation layer with each unique underlying textural and geomorphic values for all 'active' and "inactive" segments. Sixty-three joined attribute files were exported as a. csv file for further analysis.

Random Forest Classification
The segment attribute files exported in Section 3.6 were used as input data for random forest (RF) classification. The data encompass seven different SP (9,14,26,46,50,75, and 100), three layers from different frequencies (30 kHz mosaic, 200 kHz mosaic, and a combined mosaic (equally weighted) of the two frequencies) and three regions (the full Calypso HVF dataset, the North Calypso HVF subset, and the South Calypso HVF subset). Supervised RF models were run on each individual exported attribute file, executed in the "randomForest" R package (Liaw and Weiner, 2002). The RF algorithm trains an ensemble of decision trees by randomly selecting and testing a random subset of predictor variables at each node and then pooling all tree results (Brieman, 2001). The RF algorithm was chosen due to its ability to deal with unbalanced classes (as is the case with "active" and "inactive" seafloor segments (Mellor et al., 2015) and its robustness against over-fitting (Lucieer et al., 2012;Buscombe and Grams, 2018). The data sets were randomly split into two groups-training data (70%) for model development and test data (30%) for model evaluation. Binary "presence/absence" (i.e. "active" vs "inactive") classes were the target. Attribute values from each segment were the predictor variables: backscatter mean, standard deviation, and skewness; backscatter derivatives (GLCM textures), water depth, and bathymetric derivatives (BTM layers). Each model ran with 500 trees and the square root of the number of predictor variables as the subset at each node (rounded down). Error matrices and model performance statistics (Kappa, sensitivity, specificity, positive predictive value (PPV), and negative predictive value (NPV)) were calculated for each model. Error matrices compare known class data to predicted classes in a 2 × 2 table (true positive, true negative, false positive, and false negative). Model "goodness" thresholds were determined by Kappa values. This value allows for model comparison by quantifying the accuracy of each model above that of chance: 0 predictive power no better than chance, 0-0.2 poor, 0.2-0.4 fair, 0.4-0.6 moderate, 0.6-0.8 good and 0.8-1 very good, 1 perfect predictive power (Altman, 1991). Sensitivity conveys the ratio of true positive predictions (flare presence) to the total positive observations (i.e. number of correct seep predictions to the number of actual seeps). Specificity conveys to the ratio of true negative predictions (seep absence) to the total negative observations. PPV conveys the ratio of true presence predictions to the total positive predictions (i.e. number of correct seep predictions to the number of positive predictions). NPV denotes the ratio of true negative predictions to the total negative predictions. A "Mean Decrease Accuracy" measure is also generated for each model. This measure describes the relative importance of each predictor variable, by determining how much a model's accuracy degrades if a specific variable is excluded from the model. Variables with high importance drive the overall prediction outcome of an RF model.

Sediment Samples
Thirty-one Van Veen sediment grabs were collected over the survey region, 14 in North Calypso HVF, eight in South Calypso HVF, and nine in the remaining Calypso HVF region ( Figure 1B). A Van Veen sediment grab collected surface sediments, representing the top 30 cm of the seabed. Grain size distributions were extracted using a Beckman Coulter LS13 320 Diffraction Particle Size Analyzer (DPSA). Results were analyzed using Gratistat v8 (Blott and Pye, 2001). Calcium carbonate content was determined using a vacuumgasometric system (Jones and Kaiteris, 1983). Particle size distributions and D 50 (median) were measured and reported in log transformation units (ϕ) (Blott and Pye, 2001). Samples were classified into one of six descriptive sediment classes: medium silt, coarse silt, very coarse silt, very fine sand, fine sand, medium sand (Folk and Ward, 1957).

Towed Camera
Ten towed camera transects were sampled over three flares-ofinterest (FOI) (FOI-01, -02, and -03; Figure 1B white boxes), where high intensity water-column acoustic flare activity was observed during the survey. A towed video camera system recorded continuous digital color video footage (Sony HDR-CX550; Institute for Marine and Antarctic Studies). The video camera was weighted and mounted in a "fly frame" to ensure a consistent altitude above the seafloor. The camera was angled at 45°to capture gas bubble and fluid release from the seafloor, with lasers spaced at 15 cm for a scale reference. Additional GoPro cameras were attached to the towed camera system for some of the camera transects, affording multiple perspectives of the seafloor.
A USBL transponder was attached to the camera deployment line (∼0.5 m upline from camera attachment point) and linked to the RV Tangaroa Kongsberg HiPAP system. Transects were towed at ship speeds between 0.25 and 1.0 m s −1 , with the aim of maintaining an altitude of ∼1-2 m and a visible area of ∼1.5-2.5 m 2 . A total of 17 h of towed camera and 1,127 still images were obtained over the ten transects, capturing 12 km of seafloor.

Active Seafloor Fluid Expulsion Areas
A total of 3,222 individual seeps were identified over the mapped 115 km 2 Calypso HVF, between water depths of 111-292 m (Figures 1, 2). The seeps occupy ∼9 km 2 of seafloor, ∼8% of the Calypso HVF, with the remaining 92% of seafloor observed as inactive (at the time of survey). Several primary areas of active fluid expulsion were identified that agree with previous records of activity (Wright, 1992;Lamarche and Barnes, 2005;Hocking et al., 2010). These active areas occur either within or along the perimeter of ovoid or circular depressions that align with the visible NE-SW trending normal faults of Calypso 2, Pukehoko, Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 785396 8 and Nukuhou (Litchfield et al., 2013). The main North Calypso HVF activity is observed as a constellation of seeps that skirt the perimeter of two ∼1 km wide ovoid depressions that overlay the Calypso 2 normal fault ( Figure 2). The South Calypso HVF activity occurs within distinct depressions that occur between 1 and 1.5 km south of the Pukehoko North fault and between 0 and 1 km north of the Nukuhou fault. Fluid expulsion activity within the center of the Calypso HVF, the transition zone between north and south, lies between the Calypso 2 and Pukehoko North faults ( Figure 2).

Seafloor Backscatter Response
The two frequencies used to map the Calypso HVF produce mosaics that have visible differences in seafloor backscatter intensity ( Figure 3). Backscatter intensity values of the original 200 kHz mosaic range between −46 and +27 dB while the intensity values of the original 30 kHz mosaic range between −92 and 0 dB (Figure 3). These ranges included extreme outliers due to acquisition artefact ( Figure 3). Mosaic distributions were clipped to −30-−5 dB (see Section 3.4 and Figure 4). Differences between the two mosaics are most obvious in the North Calypso HVF when observed with the Jenks classification ( Figure 4). While very-high backscatter intensity patches (hardgrounds) are observed to largely correspond in both mosaics, the 200 kHz mosaic has an overall higher backscatter intensity than the 30 kHz mosaic in the vicinity of the North Calypso HVF seeps ( Figure 3C (200 kHz), d (30 kHz); Figure 4A). The 30 kHz mosaic also reveals more defined boundaries between hardgrounds and the surrounding seafloor while boundaries are more diffuse in the 200 kHz mosaic ( Figures 3E,F, 4D).
In the North Calypso HVF, an extensive assemblage of hardgrounds skirts the perimeter of two ovoid depressions, coinciding with the primary seep area; these patches are more conspicuous in the 30 kHz mosaic than the 200 kHz mosaic ( Figures 3C,D, 4D). The ovoid depressions show an inverse backscatter relationship-lower backscatter intensity in the 200 kHz (lower intensity) and higher in the 30 kHz (higher intensity). Two discrete very-high backscatter patches in the north-west ( Figures 3E,F) and west ( Figures 3G,H) largely correspond in both two mosaics.
In the South Calypso HVF, there are two large conspicuous patches ( Figures 3O,P, 4A) of very-high backscatter observable in both mosaics though more obvious in the 200 kHz mosaic. Here, a very-high intensity ring is visible in the 200 kHz but less obvious in the 30 kHz. Another two smaller patches of backscatter are observed. One is visible in the 200 kHz mosaic while almost absent in the 30 kHz ( Figures 3M,N), while the second shows the inverse relationship ( Figures 3A,B). Very-high backscatter intensity patches are also observed to correspond to fault scarps in the center (Figures 3K,L) and south ( Figures  3O,P) of the survey area. These linear shaped patterns align with the Pukehoko and Nukuhou normal faults (Figure 2).

Seafloor Sediment
The demarcation in seafloor backscatter observed between the two sites along the central scarp region (Figures 3A,B) is mirrored in the sediment grabs as a fining (decrease in median grain size) of sediments from north to south. The average D 50 (median ϕ) particle size for North Calypso HVF was 3 ϕ (125 μm) and for South Calypso HVF was 6 ϕ (16 μm). As D 50 particle size decreases, the difference in backscatter intensity between the two different frequency mosaics also decreases ( Figure 5A). North Calypso HVF sediment grab samples exhibit a range of particle sizes from very fine sand to very coarse silt (Table 2; Figures 4A,B). Of note are samples taken from within two circular depressions which are classified as coarse or very coarse silt due to their dominant silt fraction (64 and 75% weight) (Figures 5A,B). The South Calypso HVF sediment grab samples reveal a silt-dominant substrate ( Figures  4A,B), with all samples classified as medium silt ( Table 2).

Backscatter Segmentation
All 200 kHz-derived GEOBIA segments and most of the 30 kHz-derived and combined mosaic GEOBIA segments have irregular shapes for all SPs, aligning with observed natural features in both the bathymetry and backscatter (subsets shown in Figure 6 show SP46 for the NCHVF and SP75 for the SCHVF, the two best performing models for each region). Segments that do not align with natural features within the 30 kHz or combined segmentation are due to artefacts in the mosaics (i.e. predominantly artefacts that are a result of aeration of the transducer due to rough sea conditions).
As SP increased, segments that delineate very-high backscatter areas remained mostly stable while segments that delineate homogenous areas of seafloor increased in size ( Figure 6). Segments that delineate scarps also remain stable at higher SP levels ( Figures 6E,F,I,J). Though the ESP2 tool provides a repeatable and objective method for selecting SPs, we found the chosen SPs did not adequately identify boundaries of the hardground areas. Smaller SPs (SP9, SP14, and SP26), resulted in over segmentation, where too many segments were created without delineating conspicuous natural boundaries of features. These SPs resulted in highly granular segments, demarcating contiguous very-high backscatter areas into small segments.
The second-best performing model, NCHVF06 (combined mosaic, SP75, accuracy: 0.82) had a lower Kappa (0.62) though higher PPV (0.73) than the best performing model, being better at predicting seep segments. NCHVF05, the third "good" model, had highest overall accuracy (0.84) but lower Kappa (0.61) and lower PPV (0.62). The NCHVF12 model (30 kHz mosaic, SP100, Kappa: 0.52) had the highest PPV of 0.73 and NPV of 0.79. Of the SP levels selected by the ESP tool, only SP46 revealed a good model relationship.
The combined mosaic model with SP 75 (Table 4) was the best performing model for the South Calypso HVF (SCHVF) (SCHVF06; Kappa: 0.57, "moderate," accuracy: 0.81); this was also the model with highest PPV (0.61) (Figure 10). We found no significant difference between the combined model or the individual mosaic models with equivalent SP: SCHVF20 (30 kHz, SP75, Kappa: 0.56, accuracy: 0.85) and SCHVF13 (200 kHz, SP75, Kappa: 0.55, accuracy: 0.85), the second and third best performing models, respectively, for the South Calypso HVF. Of the SP levels selected by the ESP tool, none revealed a good model relationship.

Hardgrounds and Correctly Predicted Seep Presence
By applying a Jenks Natural Breaks classification to the updated distribution of mean backscatter values for each mosaic (moving from a pixel-based distribution to a segment-based distribution) we determined where high or very-high backscatter areas ("hardgrounds") are in the Calypso HVF and how they may correlate with water column acoustic flares ("active fluid expulsion").
Of the very-high 200 kHz backscatter segments, 0.3 km 2 (0.87% of the total North Calypso HVF area) were also true positives (i.e. correctly predicted active fluid expulsion segments) while 0.2 km 2 (0.67% of the total North Calypso HVF area) were false positives (i.e. incorrectly predicted active fluid expulsion) segments. Of the very-high 30 kHz backscatter segments, only 0.07 km 2 were true positives; 1.37 km 2 were false positives ( Figure 10C). The remainder comprised true negatives and false negatives.

Predictor Variable Importance
The nine most important predictor variables for the best performing model, NCHVF04 (combined 30 and 200 kHz,  Figure 7A). Skewness of 30 kHz, the distribution asymmetry of backscatter values from the pixels that construct a segment, has the highest relative importance to all other variables, meaning this variable would have the most negative impact on the model accuracy if removed. Segments with positive skew (segments that were made up of a majority of low mean dB values with a few very-high values in the tail, i.e. areas of soft sediment with flares) were aligned with areas that contained known seeps, (Figures 6D,H)

Seafloor Observations
Towed camera footage, though biased toward discrete seafloor areas of interest, provides additional information about active seafloor fluid expulsion that accords with model results. Though not spatially extensive, video footage confirms the assumption that very-high backscatter classes are due to indurated or lithified seafloor. Carbonate platforms and lithified seafloor observed in the towed camera aligned with high and very-high backscatter areas observed in the mosaics (Figure 8). This relationship was more obvious however in the 30 kHz mosaic and for FOI-02 and -03 ( Figures 8D,E , 9). No shell hash, mussel beds, or other potential reflectors traditionally also associated with fluid expulsion zones, and potential causes of high seafloor backscatter, where observed in the video footage. FOI-01 in the North Calypso HVF (Figures 8A,B) exhibits a high intensity of fluid expulsion coinciding with bioturbated sediment (Figure 9.1-3). Videos show shimmery fluid (due to density differences) escaping from the perimeters of anhydrite mounds (∼1 m diameter) covered in white filamentous bacteria (Figure 9.4-6); less intensive bubble expulsion was also observed. This fluid expulsion coincided with very-high backscatter in the 200 kHz mosaic ( Figure 8A) and broadly low backscatter with smaller patches of very-high backscatter in the 30 kHz mosaic ( Figure 8B). FOI-02 in the South Calypso HVF (Figures 8C,D) hosted a wide variety of substrates: sedimented regions with high-flux fluid expulsion (Figure 9.7); bioturbated sediment (as evidenced by active burrowing within surface sediment) with gas expulsion (Figure 9.8); hardground regions with low-flux fluid expulsion from the perimeters of, or in between, broken platforms ( Figure 9.9) or rocky piles of hardgrounds (Figure 9.10); occasional anemone gardens (Figure 9.9) and mussel beds (Figure 9.10); areas of extensive white filamentous bacterial mats with very-high-flux gas and liquid expulsion (Figure 9.11). Of note were the many instances of liquid and gas bubble expulsion at the same area of seabed, with one high flux gas and fluid expulsion area of rocky substrate blanketed by thick white bacterial mats (Figure 9.12). These hardground regions coincided with very-high backscatter intensity in both mosaics. Other very-high backscatter intensity patches in the South Calypso HVF survey area were not sampled with the towed camera.
FOI-03 in the North Calypso HVF ( Figures 8E,F), hosts fluid expulsion coinciding with a broad break of slope ( Figure 1B). This seafloor was characterized by: bioturbated sedimented regions with no fluid expulsion (Figure 9.13); white and lilac colored filamentous bacterial with and without gas bubble expulsion (Figure 9.14), with lilac colored bacteria corresponding with areas of minimal fluid expulsion (Figure 9.15); ∼0.5 m diameter anhydrite mounds with high flux liquid expulsion (observed in towed camera footage as shimmery film) (Figure 9.16); hardground regions with lower flux gas bubble expulsion, from the perimeters of, or in between, broken platforms (Figure 9.17) and rocky piles of hardgrounds; and occasional anemone gardens (Figure 9.18). Here, hardground areas corresponded to very-high backscatter intensity in the 30 kHz mosaic and also in the 200 kHz mosaic, though the overall very-high intensity signal of the 200 kHz mosaic obscured direct comparisons between veryhigh 200 kHz values and seafloor observations. More extensive and very-high backscatter intensity patches in the North Calypso HVF survey area were not sampled with the towed camera on this survey.

Quantitative Links Between Active Seepage and Seafloor Characteristics
We were able to link water-column acoustic flares (as an active fluid expulsion proxy) to specific backscatter and bathymetric derivatives of the seafloor using RF models at a shallow hydrothermal vent field. The best performing model (North Calypso HVF, combined mosaics, SP46, with overall accuracy: 0.75; PPV: 0.67; Kappa accuracy: 0.65) demonstrated the potential of using GEOBIA analysis, image segmentation, and RF modelling to predict regions of potential fluid expulsion based on seafloor characteristics ( Figure 10). Although the South Calypso models were not as accurate, they still performed well in highlighting areas where active fluid expulsion were observed ( Figure 10). Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 785396 15 We were unable however to link gas flares specifically to "hardgrounds," or very-high backscatter areas. Although areas of very-high backscatter areas in the 200 kHz mosaic were associated with predicted seep locations in the North, veryhigh backscatter dominated the whole mosaic most likely due to sand-dominant seafloor surface substrate overprinting discrete hardground areas, including areas that lacked seeps ( Figure 8A). This finding contradicted what was observed in the 30 kHz mosaic, with seeps associated with a range of backscatter values, from low to very-high. These results however, combined with those provided by the RF variable importance measure, provide us with initial seafloor indicators associated with active fluid expulsion. Naudts et al. (2008) observed that while fluid expulsion activity from shallow (<200 m water depth) seeps in the NW Black Sea and "hardground" backscatter were often correlated, higher fluid expulsion occurred adjacent to, or nearby, these hardgrounds. Along the Hikurangi Margin Naudts et al. (2010) noted that fluid expulsion from seeps in water deeper than 500 m escaped from sandy sediments away from highly reflective backscatter hardground platforms. Mitchell et al. (2018) however, supported the intuitive link between gas flares and potential hardgrounds when superimposing water column acoustic flares over anomalous seafloor backscatter intensities from ∼1,000 to 2000 m in the Gulf of Mexico. Thorsnes et al. (2019) showed on the Norwegian continental shelf (<400 m) that spatial relationships between hardgrounds and seeps are not always consistent, particularly when considering relict or dormant systems, or for fluid expulsion systems occurring from seafloor with more highly variable substrate types.

Implications for Predictive Modelling of Active Fluid Expulsion
Whilst we demonstrate from this example that areas of very-high seafloor backscatter are not good predictors of gas flares in this context, we suggest finer-scale textural indicators may be a useful tool to determine locations of active seeps. Skewness of the 30 kHz mosaic segments, whereby segments with positive skew were associated with known seeps, was the most influential variable in the best performing model and ranked highly in most of the other models. Positive skew describes the longtailed distribution of all the 30 kHz pixels that contribute to each segment, with a higher density of relatively lower pixel values with few much higher values (the positive tail). This corresponds to textures observable in the 30 kHz mosaic, for e.g. Figure 8B, where discrete patches of high to very-high backscatter values are scattered among surrounding lower backscatter values.
Towed camera footage shows a similar relationship-areas of rubble or broken hardgrounds with minimal to no fluid flux were interspersed with unconsolidated sediment that hosted the most white lines: segments of best performing models for North (SP46 combined) and South (SP75 combined). Numbers correspond to still images in Figure 9.
Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 785396 ebullient seeps (Figure 9). Although we found smaller SPs did not contribute to any of the best performing models, segment skewness may more aptly capture the fine-scale seafloor textures that are lost when using segmentation and GEOBIA methods, while still allowing for noise inherent in seafloor backscatter mosaics to be accounted for.

A Combined Frequency Approach to Improve Model Accuracy
Results from the Calypso HVF demonstrate that combining backscatter mosaics from multiple frequencies increased model accuracy compared to using individual mosaics. The combined models were the best performers for both the North and South Calypso HVFs. We suggest the improved predictive power, compared to a single frequency, is achieved by capturing information from both surface backscatter and notable acoustic returns in the immediate sub seafloor surface. Buscombe and Grams (2018) suggested that increasing dimensionality in a seafloor classification, by using multiple frequency mosaics, provides more points of understanding for a model to draw upon when determining substrate boundaries. Classification models using the single 30 kHz mosaics for North and South also outperformed their equivalent 200 kHz-based models at each of the segmentation scales, both in overall model accuracy and predictive power. In the North in particular, the 30 kHz mosaic better highlights areas of veryhigh to high backscatter than the 200 kHz, as there is difficulty differentiating between sandy substrate and hardgrounds in the higher frequency mosaic using backscatter alone ( Figure 4C). Sand grains (<4 ϕ), as observed in the sediment samples collected in the North Calypso HVF, are strong acoustic scatterers, comparable to hard substrate (Lurton and Lamarche, 2015). The lower frequency sounder was able to penetrate beneath the sandy overprinting veneer and better delineate the hardgrounds beneath. Work has been undertaken to estimate signal penetration depths, with low frequency systems having the greatest penetration potential: e.g. 4 m (Hillman et al., 2017) or 12 m (Schneider Von Deimling et al., 2013) penetration for a 12 kHz signal; up to 7 m penetration for 27-30 kHz (Mitchell,FIGURE 9 | Still images for three features of interest (FOI; see Figure 1B for FOI locations). Numbers correspond to approximate locations along towed camera transects in Figure 8. White dashed lines highlight where fluid was observed exiting the seafloor; white circles highlight observable bubbles.
Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 785396 1993), and decimeter penetration for a 95 kHz signal (Fonseca et al., 2002). This demonstrates the usefulness of employing a lower frequency sounder than would typically be considered ideal for these depth ranges (<200 m) to reveal sub-surface structures. The usefulness of using a lower frequency sounder, or combining two frequency sounders, for classifying the seafloor and interpreting sediment at depth has been demonstrated elsewhere. Hughes Clarke (2015) initially demonstrated that so long as a minimum separation of one octave between frequencies is maintained, two discrete MBES frequencies may enable frequency differences to be exploited and allow for improved discrimination between surface seafloor characteristics and subsurface information. Gaida et al. (2018) found that using a 100 and 400 kHz combination increased substrate discrimination in the Bedford Basin, Canada, compared with using a single frequency. Working on the same dataset, Brown et al. (2019), observed coarse sub-surface dredge spoils in a 100 kHz mosaic that were not observed in the 400 kHz mosaic. Feldens et al. (2018) found that lower frequency data were more sensitive to substrate changes, highlighting buried seafloor structures that were not detected in higher frequency data and highlighted the greatest contrast between seafloor facies. This allowed for additional sediment classes to be produced by combining surface acoustic information with subsurface acoustic information observed in the 100 kHz frequency backscatter data. Buscombe and Grams (2018) showed that multi-spectral data and spatially-aware machine learning algorithms that do not rely on assumptions of independence in substrate models may better determine substrate classifications. Janowski et al. (2018) however, highlighted that multi-spectral seafloor mosaics, though ideal, are not always required for seafloor discrimination. They showed a frequency pair (150 and 400 kHz) can provide enough additional frequency information to improve overall accuracy of RF models. Specific to shallow hydrothermal vent fields, penetrative backscatter differences may be useful in detecting emergent seeps or sub-surface sediment modification or classes that are difficult to detect from a conventional single frequency mapping approach, while spatially-aware machine learning algorithms may allow for the incorporation of broader geological influences such as fault scarps.

Factors Affecting Model Accuracy
We propose differences in model accuracy between the North and South Calypso HVF are the result of two different seep settings, both in surrounding substrate and seep style, and may imply two FIGURE 10 | Final model predictions for the best performing models of (A) North Calypso HVF (SP46) and (B) South Calypso HVF (SP75) with seep locations (black dots) overlain, alongside (C) segment area and very-high backscatter area confusion matrix results for the best performing NCHVF model.
Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 785396 distinct or diverging seep systems. We also speculate that a complex arrangement of precipitated minerals, including hydrothermal sulfates (e.g. anhydrites) as well as carbonates, may be producing the very-high seafloor backscatter anomalies at the seafloor surface of the two sites. Overprinting by terrigenous and volcanic sediments may also affect model accuracy due to hardground burial however sedimentation rates in the Bay of Plenty, in particular the Calypso site are relatively low (in the order of decimeters per 10,000 years) (Kohn and Glasby, 1978;Bostock et al., 2018). The South Calypso HVF seeps are represented by the veryhigh backscatter mounds surrounded by medium silt seafloor (predominantly low backscatter in the 200 kHz and medium in the 30 kHz mosaic). The North Calypso HVF in contrast exhibits a constellation of seeps over varying substrate types. The South Calypso HVF is nearer to the Whakatāne Graben axis than the North and noticeably, the South Calypso HVF seep clusters exhibit more ebullient liquid and gas flux in video footage than what we observed in the North. This nearness of the South Calypso HVF to the graben axis, its more active appearance, and more discrete emplacement style lacking extensive hardgrounds may imply a younger hydrothermal seep than the North. In addition, a higher abundance of gas at depth in the South Calypso HVF sediments may contribute to the higher intensity response observed in the 30 kHz mosaic. The relatively lower gas flux, compared to the South, and extensive sediment settling on most of the hardground platforms observed in the towed camera suggest the North hardground areas are less active and have been emplaced over a longer period of time (thousands of years) (Bowden et al., 2013). The combined extensive hardgrounds and active seeps in the North suggest a combined seep system of active hydrothermal venting and longer-lived carbonate hardgrounds. Largely, the emplaced hardgrounds closely resemble authigenic carbonate platforms observed at cold methane seeps, implying a gentler hydrothermal gradient with the potential to encourage microbial anaerobic oxidation of methane to occur. Procesi et al. (2019) pose that the Calypso HVF is a sedimenthosted geothermal system (SHGS). SHGSs are geological vent and seep systems that straddle both strictly hydrothermal and hydrocarbon arrangements; they emit both mantle derived CO 2 and biotic CH 4 . For Calypso HVF, Stoffers et al. (1999) and Kamenev et al. (1993) both reported the presence of hydrocarbons in vent fluid and bottom waters. Botz et al. (2002) measured 83.61% CO 2 and 6.196% CH 4 within the North Calypso HVF vent gas, with a temperature of 201°C. Values of 72.10% CO 2 and 9.898% CH 4 in vent gas and 51.79% CO 2 and 1.182% CH 4 in vent water, with a temperature of 181°C recorded at the South Calypso HVF. While the carbon isotope values observed between −3.4‰ and −5.5‰ of the reference standard for δ 13 C (Pee Dee Belemnite (PDB)) indicate a shallow magmatic source for the CO 2 gas, δ 13 C isotope values between −24.6‰ and -28.9‰ PDB at Calypso HVF also indicate CH 4 sourced from thermal maturation of organic matter (thermogenesis) within marine sediments (Botz et al., 2002). Sarano et al. (1989) and Kamenev et al. (1993) reported anhydrite mounds up to 8 m in height. Numerous smaller (∼1-2 m high) anhydrite mounds, with observable liquid escape (shimmery water-column features in the video footage) and associated filamentous bacteria, were observed during the 2018 voyage. No mounds were detected on the same scale as the previously observed mound (Figure 9). The small anhydrite mounds observed during the 2018 voyage may be a precursor to larger cones, if the vent pathway is maintained, or eroded versions of formerly larger features. Hocking et al. (2010) observed that the largest mound had visibly eroded over two decades, suggesting high temperature fluid expulsion may have waned or relocated in the North Calypso HVF since original observations were made by Sarano et al. (1989). Anhydrite forms as hot calcium-rich hydrothermal fluids mix with cold sulfaterich sea water (Sarano et al., 1989;Hannington et al., 2001;Lowell et al., 2003). The retrograde solubility for anhydrite is > 150°C; below this threshold, anhydrite dissolves (Mottl, 1983). Without a constant supply of hot vent fluid, anhydrite dissolves at shallow depths (<200 m) (Lowell and Yao, 2002;Lowell et al., 2003). Hardgrounds however exist on timescales many orders of magnitude longer than active fluid expulsion (Bayon et al., 2009;Liebetrau et al., 2010;Bowden et al., 2013); coexistent hardgrounds and anhydrite mounds suggest both long-lived and contemporary fluid expulsion. Consequently, modelling flares to hardgrounds may not capture such timescale differences.

Improvements and Future Applications
Fine tuning our approach may lead to more sophisticated means of locating and differentiating between active and relict fluid expulsion areas. Additional recommendations for future mapping might consider the following elements in each step of the modeling pathway.

Acquisition
A stratified random sampling regime, combined with more extensive video footage to observe first order indications of active and inactive fluid expulsion, would allow more robust associations between gas flares/fluid expulsion and seafloor characteristics. This would allow the unique seafloor textural differences associated with both to be captured and fully modeled. A single-pass survey using an integrated multi-spectral system for multi-frequency seafloor analysis (Brown et al., 2019) would ensure MBES calibration, correction, and data reduction issues could be reduced by collecting varying frequency information on a ping-by-ping basis. Further, repeat surveys to capture change detection in both the seafloor and water column acoustic backscatter would allow us to integrate short time frame fluctuations in a dynamic seep system into the model.

Processing
Seep identification using higher frequency (200 kHz) watercolumn acoustic data would allow us to interrogate more forms of fluid expulsion. We identified fluid expulsion using 30 kHz water-column acoustic data which favors gas expulsion. Additional fluid expulsion evidence may be found in the 200 kHz water-column data, e.g. liquid expulsion that lacks a gas fraction.
Image segmentation and pre-classification could also be used to remove obvious artefacts prior to segmentation.

Modelling
Sampling hardground platforms in both regions, to determine sediment mineralogy, combined with acoustic information, discrete camera sampling, and geochronology could be used to determine the genesis or timescales of the active fluid expulsion regions. Additionally, seep differences could be modeled to further understand seafloor modifications possible within of a sediment-hosted geothermal system.

CONCLUSION
We demonstrate that a dual frequency approach, within the bounds of a MBES survey, leads to improved seep prediction. Combining the penetration and resolution of two distinct MBES frequencies can provide added value for the interpretation of seafloor structure from backscatter data. Understanding the potential of different acoustic frequencies when targeting specific seafloor substrates is necessary to inform decisions of vessel choice or limitations when preparing for a survey of a key substrate type. We show that extensive and precise groundtruthing (through video, grabs, and cores) that enable strong associations between paired or multi-frequency seafloor maps are required to advance our understanding of penetration depth and influence of sounder choice on seafloor substrate classification. Understanding the limits imposed by chosen frequencies and sonar systems will increase these associations, while improving integration methods for paired frequency systems.
The main outcome of this study is that the extensive hardgrounds of the Calypso HVF are not directly linked with gas flares within this shallow hydrothermal system, while unconsolidated sediment nearby hosts higher flux expulsion. This aligns with the development and cessation of other submarine fluid expulsion system (Naudts et al., 2008;Judd and Hovland, 2009;Naudts et al., 2010). Seafloor fluid systems will refashion themselves on decadal timescales due to vent instability, clogging from deposited minerals, modification of the sediment, or rearrangement due to ongoing extensional processes and tectonic stresses (Desbruyères, 1998;Bowden et al., 2013;Levin et al., 2016). Once a seep self-seals, fluid will migrate to a more permeable pathway. Being able to remotely predict active and inactive regions of fluid expulsion will prove a useful tool in rapidly identifying seeps in legacy datasets, as well as textural metrics that will aid in locating nascent, senescent, or even extinct seeps during a survey, where important communities may persist even if they appear inactive. Understanding the seafloor qualities best able to predict active or inactive fluid expulsion will allow us to establish primary controls on fluid expulsion in a shallow marine environment and further understand seep succession stages and evolution.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
All authors give final approval of the article to be published. ES conceived and designed the project; processed and analyzed data; generated segmentation and random forest analysis output; substantially drafted, constructed, and revised the article and figures. GL was voyage leader of QUOI-TAN1806 voyage that produced the data and was science leader of public funded program that funded the research. GL and VL funded the voyage from project "Building capability for in situ quantitative characterization of the ocean water-column using acoustic multibeam backscatter data" (2017-2018) -Royal Society of New Zealand Catalyst Fund. JW, VL, and GL supervised the project; provided analytical support and guidance; critically reviewed and revised the article. SW acquired and processed data; provided technical and analytical support; revised the article. YL acquired and processed data; provided technical support; revised the article. EH acquired and processed data; provided technical support; revised the article. AP acquired and processed data; provided technical support; revised the article.