Quantification of the fine-scale distribution of Mn-nodules : insights 1 from AUV multibeam and optical imagery data fusion 2 3 4

Autonomous underwater vehicles (AUVs) offer unique possibilities for exploring the deep seafloor in high resolution over large areas. We highlight the results from AUV-based multibeam echosounder (MBES) bathymetry / backscatter and digital optical imagery from the DISCOL area acquired during research cruise SO242 in 2015. AUV bathymetry reveals a morphologically complex seafloor with rough terrain in seamount areas and low-relief variations in sedimentary abyssal plains which are covered in Mn-nodules. Backscatter provides valuable information about the seafloor type and particularly about the influence of Mn-nodules on the response of the transmitted acoustic signal. Primarily, Mn-nodule abundances were determined by means of automated nodule detection on AUV seafloor imagery and nodule metrics such as nodules m−2 were calculated automatically for each image allowing further spatial analysis within GIS in conjunction with the acoustic data. AUV-based backscatter was clustered using both raw data and corrected backscatter mosaics. 
 
In total, two unsupervised methods and one machine learning approach were utilized for backscatter classification and Mn-nodule predictive mapping. Bayesian statistical analysis was applied to the raw backscatter values resulting in six acoustic classes. In addition, Iterative Self-Organizing Data Analysis (ISODATA) clustering was applied to the backscatter mosaic and its statistics (mean, mode, 10th, and 90th quantiles) suggesting an optimum of six clusters as well. Part of the nodule metrics data was combined with bathymetry, bathymetric derivatives and backscatter statistics for predictive mapping of the Mn-nodule density using a Random Forest classifier. Results indicate that acoustic classes, predictions from Random Forest model and image-based nodule metrics show very similar spatial distribution patterns with acoustic classes hence capturing most of the fine-scale Mn-nodule variability. Backscatter classes reflect areas with homogeneous nodule density. A strong influence of mean backscatter, fine scale BPI and concavity of the bathymetry on nodule prediction is seen. These observations imply that nodule densities are generally affected by local micro-bathymetry in a way that is not yet fully understood. However, it can be concluded that the spatial occurrence of Mn-covered areas can be sufficiently analysed by means of acoustic classification and multivariate predictive mapping allowing to determine the spatial nodule density in a much more robust way than previously possible.


Introduction 1.1 Mn-nodules exploration
Research on Mn-nodules received increased attention in the last decade due to increasing prices for ores rich in Cu, Ni or Co, i.e. metal resources that are contained in Mn-nodules.In nature, the largest Mn-nodule occurrences are found in the deep sea, e.g. the equatorial Pacific between the Clarion and Clipperton fracture zone (CCZ), the Peru Basin as well as the Atlantic and Indian Ocean (Petersen et al., 2016 ).In the typically muddy sediments of the deep sea, Mn-nodules form an important hard substrate providing a habitat for deep sea sessile fauna such as sponges, corals and associated organisms (Vanreussel et al., 2016;Purser et al., 2016).Therefore, mapping Mn-nodule fields is a two-fold task, comprising not only the assessment of Mn-nodules and their density distribution for accurate resource assessment, but also the improved understanding of the natural habitat heterogeneity and its relation to the deep sea ecology.Knowledge about Mn-nodule habitats will support mitigation strategies for mining-induced impacts.Since an increasing number of countries move forward with exploitation plans for Mn-nodules in the CCZ, strategies for a detailed mapping of the deep sea Mn-nodule fields might become mandatory in order to proceed with licensing procedures prior to any mining activity.
Deep sea mining will cause substantial disturbances of the deep sea ecosystem since Mnnodules, the primary hard substrate, will be removed and massive re-sedimentation of the top 20 to 30cm of sediment of the mined area will occur (Bluhm et al., 1995, Vanreussel et al., 2016).Thus, efforts have been made to investigate the effects of potential mining disturbances in the past (e.g.Thiel et al., 2001) and currently during the project "Ecological Aspects of Deep Sea Mining" as part of the Joint Programming Initiative Healthy and Productive Seas and Oceans (JPI Oceans).To study in detail the potential effects of a deep sea disturbance by Mnnodule mining to benthic fauna, a plough-experiment was performed in 1989 in the Peru Basin as part of the DISturbance and reCOLonization project (DISCOL, www.discol.de).A plough of 8m width was towed 78 times over a 2nmi wide circular area (February-March 1989) to generate dense and less dense impact sub-areas.Photographic surveys, sediment and biological sampling before and after the disturbance (September 1989, March 1992, February 1996), showed that the plough marks were well visible even after 26 years and that the benthic fauna did not recover to its initial state.The data used in this study were collected during the SO242-1 cruise to the DISCOL area during summer 2015, 26 years after the DSICOL experiment.

The DISCOL study area
The DISCOL working area is situated 560 nmi SW of Guayaquil on the Pacific Oceanic Plate in the Peru Basin (Fig. 1A) in about 4150 m water depth.The larger DISCOL area ranges from 3800m to 4300m water depth (Fig. 1B) and is characterized by N-S oriented graben and horst structures with a deep N-S elongated basin with water depths down to 4300m.An 11 km wide seamount complex in the NE along with a second seamount complex to the SW and three higher mounds to the NW clearly show that the DISCOL area is not located on a flat and homogenous deep seafloor.
The ploughed DISCOL Experimental Area (DEA) itself is located on a relatively smooth, slightly elevated part of the seafloor with a central valley of about 20m depth that dips southward (Fig. 2A).When inspecting the bathymetry data generated by the autonomous underwater vehicle (AUV) in more detail, the central part of the area shows a 20m deep valley, the floor of which is comprised by low-relief N-S trending ridges giving the impression of a braided river system (Fig. 2A).Despite the rich morphological features in the study area, it does not contain steep slopes and represents a rather smooth seafloor (<5 degrees).

Acoustic mapping of Mn-nodules and study objectives
Acoustic mapping has proved to be a useful tool for supporting deep sea mineral resource assessments.The initial studies mentioned below, showed promising results for Mnnodule detection and quantification, however, progress in more detailed and meaningful method development and data processing capabilities has remained slow, mainly due to fluctuations in the global interest of deep sea mining.The majority of surveys performed for Mn-nodule mapping purposes rely on acoustic remote sensing and near-bottom photography (de Moustier, 1985).The applicability of acoustic methods is based on the clear acoustic contrast of at least 11 dB between the background deep sea soft sediment and the nodules (de Moustier 1985).Weydert (1985) found that the nodule size is proportional to the average backscatter strength for low frequency signals (<30 kHz).In addition, Weydert (1990) concluded that it is possible to map the percentage of seafloor covered by nodules based on backscatter measurements of sonar frequencies higher than 30 kHz , whereas for a frequency of 9 kHz it is possible to use the backscatter response to determine whether the nodule diameter is greater than 6 cm or smaller than 4 cm.Masson and Scanlon (1993) suggested that lower sonar frequencies produce a much weaker acoustic contrast between nodules and surrounding sediments for nodules of given size.They concluded that on a seafloor covered with mixed-size nodules larger nodules will have a greater impact on the backscattered energy than smaller ones.They also suggested that minor differences of nodule coverage will have a considerable effect in backscatter values.A more recent study by Chakrabotry et al. (1996) suggested that the nodule coverage is proportional to the backscatter strength and that for low frequency (15 kHz; wavelength ca. 10 cm) the main type of scattering is Rayleigh scattering (wavelength/10 < nodule size) for nodules and coherent scattering for fine sediments.During one of the first deep sea studies for acoustic mapping of Mn-nodules, de Moustier (1985) utilized a multi-beam echo-sounder (MBES) sonar combined with near-bottom acoustic measurements and photographs from a deep towed camera system to infer nodule coverage.He managed to obtain high agreement between relative backscatter intensity classes and three types of nodule coverage as interpreted from seafloor imagery (dense, intermediate and bare).At that time, his results highlighted the great potential of MBES technology in deep sea mineral prospecting.In more recent years Lee and Kim (2004) utilized side-scan sonar (SSS) to examine the relation of regional nodule abundance with geomorphology.According to their qualitative analysis, lower backscatter values are related with abyssal troughs whereas increased backscatter values are related to abyssal hills.Additionally, Ko et al. (2006) attempted to examine the relation between MBES bathymetry and slope with nodule density in the equatorial Pacific without identifying a solid pattern.Most recently, Okazaki and Tsune (2013) utilized AUV-based MBES, SSS and image data for Mn-nodule abundance assessment and its relation to deep sea micro-topography.More recent projects regarding resource assessment of Mn-nodules at large scales (0.1' by 0.1' grid cell size) have been based on various spatial modelling and decision making techniques (ISA, 2010).Most commonly, the kriging method has been applied on sparse ground truth data (obtained by physical box-corer sampling) while logistic regression and fuzzy logic algorithms were applied in multivariate data sets of Mn-nodule-related environmental variables such as sediment type, sea surface chlorophyll and Ca Compensation Depth (CCD) (Agterberg & Bohnam-Carter, 1999, Carranza & Hale, 2001).In this study we analyse AUV-based MBES and image data for quantitative mapping of Mnnodule densities in the Peru Basin.Particularly, we utilize local ground-truth information (Mnnodule measurements from AUV photographs) in order to investigate a) its relation to acoustic classification maps and b) its potential use for predictive mapping of Mn-nodules in wider areas where only hydro-acoustic information is available.Therefore, we apply two unsupervised methods (Bayesian probability and ISODATA) for seafloor acoustic classification and a machine learning algorithm (Random Forest) for Mn-nodule density predictions beyond the areas that were optically imaged using the AUV.By applying different algorithms for unsupervised classification, we aim at comparing their results against quantitative ground truth data of nodule metrics from automated analyses on AUV imagery.This way, we will assess the ability of classification methods in discriminating areas with distinct nodule densities.To our knowledge, this is the first time the Random Forest algorithm is applied for predictive mapping of Mn-nodule densities.Therefore, we examine its performance and the influence of various AUV MBES data on the Mn-nodule prediction results.

AUV MBES data acquisition and processing
The data in this study were collected using the AUV "Abyss" (built by HYDROID Inc.) from GEOMAR, during cruise SO242-1 where various AUV missions were flown.The AUV is equipped with a RESON Seabat 7125 MBES sensor with 200 kHz operating frequency, 256 beams with 1 by 2 degree opening angle along and across track, respectively.From the original PDS2000 sonar data, files backscatter snippet data were extracted into s7k format whereas bathymetry data were extracted into GSF format.Prior to exporting, MBES bathymetric data were filtered within the PDS2000 software.Bathymetry data from different AUV dive-missions were jointly used for interpolating one single grid of bathymetry and backscatter (Fig. 2).Latency and roll-related artefacts affected bathymetry in places due to a none-constant time delay for roll values creating uncorrectable artefacts in the resulting grid.Therefore, the bathymetry was smoothed by applying a Gaussian filter with a 10 m x 10 m rectangular window with 3 and 5 standard deviations as smoothing factors in SAGA GIS.Filtered bathymetry was visually inspected for artefacts using the hill-shade function in SAGA GIS, giving satisfactory results.Vertical differences between the smoothed grid with the originally processed surface were everywhere less than 1 m, highlighting that the filtering did not cause significant smoothing and removal of finer details.The filtered bathymetric grid was used for calculating a variety of derivatives listed in Table 2.
The MBES backscatter data were processed in two ways.First, the s7k/GSF pairs were automatically corrected (for radiometric and geometric bias) and mosaicked in QPS FMGT (Fig. 2B).In addition, backscatter mosaic statistics were calculated and exported as GEOTIF files using a 10 m x 10 m neighbourhood.The raw snippets data were exported prior to any processing using a combination of in-house conversion software and QPS DMagic for merging beam data with ray-traced easting and northing.The raw snippets data were transformed from 16-bit amplitude units to dB using the formula in Eq. (1): Backscatter (dB) = 20*log 10 (amplitude) (1) Raw backscatter data were processed by applying the Bayesian approach on certain beams as described in Alevizos et al. (2015 and2017) whereas the gridded data were analysed with Random Forest (RF) regression trees and ISODATA clustering (see section below).An overview of the software used to process and classify each type of dataset is presented in Table 1.

Seafloor imagery and automated image analysis
AUV surveys were undertaken for collecting close-up images from the seafloor using a camera system recently described by Kwasnitschka et al (2016).In this system the camera is mounted behind a dome port along with a 15mm fish-eye lens that produces extreme wide-angle images.This type of lens and dome port configuration induces significant distortions to the image which need to be corrected prior to any image analysis processing.Surveying at altitudes of 4-8m above the seafloor and using the novel state-of the-art LED flash system, the AUV collected several hundred-thousand seafloor images at a 1Hz interval.The respective AUV surveys were designed to cover a large part of the study area with a single-track dive pattern and also to focus on two selected areas running track lines 5m apart for dense 2D image mosaicking (Fig. 2A).Each image was individually georeferenced using the AUV navigation and altitude data.This way, each pixel of the AUV imagery is translated to an actual portion of the seafloor.For the automated image analyses (e.g.Mn-nodule counting), all images were smoothed by a Gaussian filter to remove noise and then converted to grayscale for computational speedup.Following, the images were corrected for inconsistent illumination due to the varying AUV altitude using the fSpice method described by Schoening, et al. (2012).The central (sharpest, best illuminated) region of each image was cropped and thresholded by an automatically tuned intensity limit before contours in the resulting binary images were detected and fused to blobs of pixels that served as nodule candidates.Each nodule candidate was finally fitted with an ellipsoid to account for potentially buried parts of the nodule.The sizes of these ellipsoids constitute the nodule size distribution within one image from which descriptive parameters were derived.This kind of automated image processing resulted in quantitative information such as: image area (square meters), number of nodules (n), percentage of seafloor covered by nodules (amount of nodule pixels divided by total amount of image pixels), and the threshold sizes (estimated 2D surface) of 1, 25, 50, 75 and 99 percent quantiles of the nodule size distribution (comparable to a particle size analysis).A detailed publication on the nodule delineation algorithm can be found in Schoenning et al. ( 2017), while the source code is available online as Open Source (https://doi.pangaea.de/10.1594/PANGAEA.875070)In this study, we considered the number of Mn-nodules per square meter as a normalized measure of nodule density in order to avoid overestimation of Mn-nodules due to multipledetections between overlapping images.This metric is derived from the ratio of the number of nodules detected to the area (m 2 ) of the image footprint (the size of the central 'good' part of the image).Therefore the results of the predictive mapping are presented with 6 m x 6 m resolution which is representative for the majority of image footprint sizes.

Seafloor classification and prediction methods
Three different approaches were applied for a predictive Mn-nodule mapping.The first approach is an unsupervised method based on Bayesian statistics applied on raw snippet data.It examines the within-beam backscatter variability in the entire area in order to estimate the optimum number of seafloor classes.The output acoustic classes can then be validated with available ground-truth data.The second approach, is based on the ISODATA algorithm (an unsupervised method as well), applied on gridded backscatter data.This algorithm can automatically adapt the number of classes to the data for given minimum and maximum values set by the user.Finally, a supervised machine learning method was applied on gridded bathymetric and backscatter data.This method requires a training set in order to model the complex relationship between the Mn-nodules occurrences and the bathymetry, bathymetric derivatives and backscatter information.The algorithm outputs a prediction grid for Mnnodule densities and also estimates the importance of each input variable in accurately predicting Mn-nodule densities.

Random Forests (MGET)
Table 1: Datasets and methods applied in this study.

Bayesian probability on beam backscatter
The raw backscatter data were classified by applying the Bayesian methodology developed and implemented by Simons and Snellen (2009) and Amiri-Simkooei ( 2009) and applied by Alevizos et al. (2015).In order to enhance the method's performance, strong outliers in the raw data were filtered by using a variance threshold set to 100 (i.e. 10 standard deviations).Thus, beams with a snippet data variance greater than 100 were disregarded from the classification process.The remaining snippet data were averaged for each beam for obtaining the mean relative backscatter intensity.The Bayesian method is based on the central limit theorem and the assumption that acoustic backscatter measurements of a homogeneous seafloor type would express normal distribution when derived from a certain incidence angle.Therefore all backscatter values were grouped per beam angle and their histograms were examined separately.At first, a number of Gaussian curves were fitted to each histogram and the goodness of fit was assessed by the χ 2 criterion.The minimum number of Gaussian curves that fitted well the overall distribution pattern of the histogram values (i.e.: χ 2 is less than 2), was considered as the optimum number of classes.Not all beam angles provided the same number of Gaussian curves; therefore it was important to identify those beam angles that gave consistent results about the number of classes.Usually the mid-range incidence angles provided the most consistent results (Alevizos et al., 2015) regarding the Gaussian fitting; hence beams from this range were utilized as reference in order to derive the optimum number of classes.Once the reference beams were identified, the mean and standard deviations of each Gaussian curve were used as conditions for classifying the backscatter values for the rest of the beams.
The Bayesian technique does not require the MBES to be calibrated and allows for class assignment per beam, thus maximizing the spatial resolution of the final map.The most important aspects of the Bayesian technique are the internal cluster validation based on χ 2 criterion and the increased geo-acoustic resolution, allowing for maximal acoustic discrimination of similar seafloor types (Alevizos et al., 2015).

ISODATA classification for grids
The ISODATA classification was applied to the backscatter mosaic and its derived statistics (Table 2) using the ISODATA algorithm implemented in SAGA GIS.ISODATA stands for Iterative Self-Organizing Data Analysis and has been applied in several marine mapping studies involving backscatter information (Diaz, 1999;Hühnerbach et al., 2008;Blondel and Gomez-Sichi 2009).The fundamentals of ISODATA processing are described in detail by Dunn (1977) and Memarsadeghi et al. (2007).A particular advantage of this method apart from its fast execution is that it estimates a suitable number of classes by dividing clusters with large standard deviations and by merging similar clusters at the same time (Diaz 1999).This is done automatically and the user only defines an empirical minimum and maximum number of classes.

Random Forest predictive mapping for grids
To exploit the full range of MBES gridded data and for comparison purposes, supervised classification was applied to the bathymetry, bathymetric derivatives and backscatter statistics (Table 2).Applying a machine learning algorithm was encouraged due to the abundant groundtruth data (nodule metrics from automated image analysis) and the high resolution of the various MBES layers.The Random Forest algorithm as implemented in the MGET toolbox for ArcGIS was used (http://mgel2011-kvm.env.duke.edu/mget).Initially developed by Breiman (2001) it has shown good results in marine predictive habitat mapping (Stephens and Diesing 2014, Lucieer et al., 2013, Che-Hasan et al., 2014).The algorithm requires a training data set with the response variable (here: nodule density from AUV imagery analysis results) and a set of explanatory variables (here: bathymetry, bathymetric derivatives, backscatter) as inputs in order to model the relationship between them.The training set provides the required "knowledge" about the response variable and its corresponding explanatory variable's values.At the next stage, an ensemble procedure based on several regression trees of random subsets of the explanatory variables is iteratively applied for classifying/predicting Mn-nodule density per grid-cell using a-priori information from the training sample.The prediction at a certain grid-cell is defined by the majority votes of all random subsets of trees (Gislason et al., 2006).During the iterative processing, the Random Forest will reserve randomly selected parts of the training sample for internal cross-validation of the results (out-of-bag sample).During each iteration, one explanatory variable is neglected and its importance score is calculated according to its contribution to the resulting prediction error.The variable importance calculation is considered one of the main advantages of the Random Forest algorithm.An important step prior to Random Forest application is data exploration.With data exploration it is possible to identify which explanatory variables are capable to discriminate patterns of nodule density in the study area better.A standard approach is to explore the probability density function of the response variable with each of the other gridded variables (e.g.slope, BPI, etc.).These plots give first indications about the distribution type of the response variable for a given explanatory variable.The explanatory variables presented in Table 2 were chosen as good descriptors of nodule density in the area based on the probability density functions of arbitrarily chosen classes of nodule density (Fig. 3, Appendix).The arbitrary classes where based on the quantiles method for classifying the nodule density histogram.It has to be noted that the arbitrary classes were used only for data exploration and not for the prediction of nodule densities.All descriptorgrids were resampled to 6 m x 6 m pixels in order to be compatible with the average effective area of the AUV images upon which nodule metrics were computed.An appropriate selection of training samples is fundamental for modelling the relationship between the response variable and the gridded descriptor data.Particularly, the training samples need to span the entire range of the study area capturing most of the data variability.They have to contain as diverse values as possible regarding both the nodule density and the corresponding gridded descriptor data.

Explanatory variables Description
From bathymetry Scale: 6 m cell size Depth AUV MBES, smoothed with Gaussian filter (5σ) Slope ArcGIS slope algorithm in percent units

BPI
Relative position of pixels compared to their neighbors.Inner radius 10m, outer radius 100 m (Iwashahi and Pike, 2007) SAGA GIS terrain analysis toolbox

LS factor
The integrated slope length and inclination, formula from Moore et al. (1991)

Automated nodule detection from AUV images
The automated nodule detection algorithm results for nodule density (number of nodules m -2 ) are shown in Fig. 4. The dense point cloud offers a detailed view of the nodule spatial distribution which can significantly enhance the interpretation of nodule density in conjunction with MBES bathymetry.In Fig. 4 the nodule density fluctuates in a pattern of alternating bands.By colorizing the seafloor surface and the bathymetric profile cross-section according to nodule density values, it can be seen that higher nodule densities appear on smooth slope features where the seafloor appears locally concave or terraced and also on the foot of these slopes which appear relatively lower compared to the surrounding area.By colouring the AUV bathymetry according to the nodule density it became clear that MBES derivatives may be useful for quantifying the nodule distribution in the entire study area.We thus calculated bathymetric derivatives such as BPI, concavity, slope and slope-related derivatives (LS factor, TRI) to be included in predicting nodule densities.

Bayesian acoustic classification of raw BS data
The Bayesian method identified six classes based on the analysis of beams with incidence angles between 38 and 42 degrees (Table 3).Despite the variance-based filtering, it was not possible to compensate for the remaining effects on beam incidence angles in the middle range and towards the nadir.We believe that these effects are responsible for the stripe-like classification at the outer part of the swath.The selection of six classes resulted from the agreement between two adjacent beams (Table 3) and the relative lower overlap of the Gaussian curves.The finally derived classes are ordinal; meaning that from class 1 to class 6 there is an increase in backscatter intensity.The spatial distribution of the acoustic classes expresses a gradient of high to low backscatter classes in the N-S direction (Fig. 5A).

ISODATA applied to BS data
The ISODATA algorithm was applied to the mean, mode, 10% and 90% quantiles of the backscatter mosaic.These datasets are considered more suitable than the raw backscatter data, as they hold a more realistic representation of backscatter spatial variability and they are slightly correlated (correlation coefficients: 0.5-0.9)with the mean backscatter.The ISODATA algorithm was set to produce an optimal number of clusters for different ranges of cluster amounts (minimum number of clusters from 2 to 5; maximum number of clusters from 6 to 10).The results for all possible pairs regarding the minimum and maximum clusters were divided, indicating five or six clusters as optimal.To have comparable results with the Bayesian method, six clusters were selected for further analyses.Although the algorithm does not output classes with ordering, the ISODATA classes were reclassified based on their nodule statistics to be comparable with Bayesian results (see discussion section).The classes show a decreasing amount of nodules from north to south with the nodule-free areas being sufficiently demarcated.

Random Forest predictions using bathymetry derivatives and BS data
The RF was performed in two steps: the training and the prediction step.First a sensitivity test was carried out using different percentages of training samples (Fig. 6B) and fitting models with 200 and 1000 trees.This test is essential for examining the optimal settings prior to applying a predictive model.It also helps in quantifying the stability of results (given the random character of the process) by running the model with optimal settings repeatedly.For quantifying the model accuracy we used the percentage of variance explained by the out-ofbag samples (RF algorithm output report) whereas for assessing the prediction results, calculation of R 2 was applied for measuring the correlation between the predicted and measured nodule density.According to the sensitivity analysis, a training set with 30% of the total amount of images with Mn-nodule statistics was sufficient to explain more than 70% of the variance of the out-of-bag sub-sample when training 200 trees.It was also found that this accuracy value is not improving significantly when increasing the training sample size (Fig. 6B).By maintaining the same amount of training samples (30% of the total images acquired, ca.2700 images) while using ten different parts of the data as training sample (ten-fold cross-validation), the model performance was relatively consistent (69-72%) regarding the out-ofbag variance explained (Table 4).These results refer to the Mn-nodules m -2 analyses.In addition we tested the predictability of the 2D size of nodules using the 50% and 75% quantiles of 2D sizes in square centimetres.The resulting out-of-bag variance explained was found to be much lower (35-40%), independently from the number of trees and the size of the training sample set.By using the results from the ten-fold cross-validation (or sensitivity test) we extracted the mean importance score of each bathymetry and backscatter parameter (Fig 6C).Considering the prediction of Mn-nodules m -2 , the mean backscatter data was found to be the most influencing variable which constantly scored first, followed by the BPI, bathymetry and concavity.After the sensitivity test an optimal model using 30% of all images as training data and growing 200 trees (1000 trees did not produce better results) was developed using the explanatory variables for prediction of nodule densities.The final results of the RF method express a gradient from higher to lower nodule densities from North to South (Fig. 6A).An independent subsample of nodule measurements was used for validating the prediction results.This validation sample consists of measurements selected at least six meters away from any training location, to avoid the introduction of autocorrelation effects on the validation process which could overestimate the performance of the model.The value of 6m was selected as the majority of images cover a 6 m x 6 m area on the seafloor.A comparison between the image-based Mn-nodule measurements and the averaged predicted values based on ten different RF runs show a good average correlation based on the R 2 coefficient (Table 4).This implies that there is a correlation between Mn-nodule density and MBES data, although there is some degree of uncertainty that remains in the prediction model (see Appendix

Discussion
Our results show that AUV imagery is capable to provide detailed information about Mnnodule densities hence assisting quantitative mapping of the Mn-nodule distribution on the seafloor.Consistency and repeatability of quantitative methods are fundamental factors in mapping studies and therefore automated image analysis is crucial in this regard.Expert assessments of several tens of thousands of images are practically not possible in a reasonable time frame and include a high rate of subjectivity.Thus, automated analysis of imagery is regarded as a very suitable method for quantitative mapping of Mn-nodules.This however comes at the cost that usually AUV image surveys are spatially restricted due to the low altitude above the seafloor.For larger scale quantitative mapping of nodule fields, AUV imagery data need to get spatially linked with AUV hydro-acoustic data supporting with data from all regions of interest at the seafloor.Results from image analysis can then be used as alternative information for acoustic class validation and predictive mapping.Although image analysis results do not constitute ground-truth information they are the best available data to correlate with acoustic classification and prediction results.By exploring the relationship between Mn-nodule data with bathymetry, bathymetric derivatives and acoustic backscatter, we aim to identify potential linkages that allow extrapolation of nodule information to larger areas to assess mineral resources, determine benthic habitats or learn about geological processes that might influence nodule growth.The following paragraphs discuss the performance of the applied classification and prediction methods highlighting the potential use of high resolution Mn-nodule density maps by considering various sources of errors induced throughout the data analyses.

Fine scale spatial variability of Mn-nodule density
Both, the unsupervised classifications (ISODATA, Bayesian) and the random forest prediction results are largely comparable to the nodule detection measurements map (Fig. 7).Hence, both classification and prediction data, and nodule measurements reflect a similar spatial distribution pattern of nodule densities.The Mn-nodule densities seen in the imagery highlight a pattern of alternating high and low density bands on bathymetric slope features.According to studies on the fine scale (tens of meters) distribution of Mn-nodules as summarized by Margolis and Burns (1976) higher nodule densities are related to hilltops, slopes and the foot of slopes.The authors particularly highlighted that e.g.nodule sizes vary significantly over short distances; unfortunately there were no methods to capture this variability sufficiently at the time of this study.The correlation to the bathymetry is supported by the variable importance plot of the RF model (Fig. 6C).This plot shows that both bathymetry and backscatter features contribute significantly to the prediction of the Mn-nodule densities with variables such as mean backscatter intensity, fine scale BPI, and concavity as good predictors.The predictive potential of these variables needs to be validated in future studies using MBES data from different study areas.
Both unsupervised acoustic classes and the Random Forest prediction suggest a gradient of decreasing nodule densities from north to south while the RF quantitative map (Fig. 6A) shows more gradual changes regarding the fine-scale spatial distribution of Mn-nodules.The northern part of the MBES survey is located very close to, and partly within, a seamount area.According to towed camera video footage these seamounts comprise ancient volcanoes that are now covered with deep sea fine sediments.In addition, a few pillow-basalt outcrops were found along with basalt slabs being exposed on the seamount slopes.Greater nodule densities can be observed from these images suggesting that accumulated nodules or exposed basalt rocks may be assigned to the same acoustic class that represents higher acoustic intensities.In the random forest prediction, high nodule densities could be confused with basalt rock as well (Fig. 6A, black arrows).Video data can be used in order to differentiate these seafloor types in the acoustic classes.Greater nodule densities in the vicinity of the seamounts area can be explained by the findings presented by Vineesh et al. (2009) and Sharma et al., (2013).These two studies propose that in the proximity of abyssal hills and slopes, abundant basalt fragments act as nodule nuclei that favour nodule development.Away from the seamount area, the nodule density variations follow a banded pattern of high and low density alternations with localized depressions representing nodule-free areas (Fig. 2B).The bandpattern variation is not fully understood by the datasets available in this study; however, it is assumed that it is the result of a combination of the deep sea benthic boundary layer hydrodynamics, local sediment movement and active tectonics that impacts pore fluid migration.It is not clear why and how the nodule-free areas are formed and why we observe moderate nodule densities in broad deep plains of the area.Margolis and Burns (1977) suggest that bathymetric valleys are more influenced by sedimentation hence not favouring nodule growth, but that hill tops and bathymetric slopes are covered by a greater amount of nodules due to a lower impact of local sedimentation.Whether this explanation is also true for the described study area remains speculative.In any case, backscatter data clearly indicate where areas of higher and lower Mn-nodule densities exist, allowing for future investigations of the underlying factors.

Assessing the Mn-nodule acoustic classification
To assess the performance of unsupervised classification methods in clustering homogeneous areas of Mn-nodules, we examined the within-and between-class variability of the Mnnodules densities (nodules m -2 ).The assessment is based on the descriptive statistics of nodule measurements from each class (Table 5) and box-plots of nodules m -2 from each class (Fig. 8).The box-plots assist to better illustrate the separation between classes as well.
To evaluate the separation of Mn-nodule densities that fall within different acoustic classes (Bayesian and ISODATA), we performed a Welch ANOVA along with a Games-Howell test for testing whether the mean values between the classes differ significantly.This test was selected, because the Levene's test (Martin & Bridgmon, 2012) indicated that there is no homogeneity between the class variances for both classification methods (p<<0.05).
Particularly the results of the Welch ANOVA for nodule populations belonging to the same Bayesian class (F(5,905)=700, p=<<0.05)and ISODATA (F(5, 2520)=810, p<<0.05)support the finding that the mean values of Mn-nodules densities differ significantly between the different classes.This finding supports that classification results effectively resolve acoustically homogenous areas of nodule patches which are statistically distinct to each other.Regarding the Bayesian classification results, the ordinal type of the classes can be noticed both in the statistics and the box-plots (Table 5, Fig. 8A).The mean and median values of nodules m -2 are increasing with increasing class number suggesting that higher backscatter values are related to higher nodule densities.Class 1 represents the lowest nodule densities but without including samples of zero nodules, this would make this class more distinguishable with an even lower mean value.Some class overlap can be observed in the box-plot for the Bayesian classes; the within-class standard deviation is increasing with acoustic class number, suggesting larger ambiguity for areas with increased nodule density.Classes resulted from the ISODATA clustering hold similar standard deviations suggesting a similar degree of within-class variability.Overall, Mn-nodule density classes express high within-class variability with almost 50% of within-class measurements spanning in a wider range of values causing class overlap (Fig. 8).This can be attributed to few factors such as inaccurate navigation between the different AUV deployments, shortcomings of the image-based nodule detection algorithm and noise in the backscatter data (see Appendix).However, it can be inferred from the box-plots for each unsupervised method that seafloor areas of homogeneous Mn-nodule density can be discriminated by classifying the MBES backscatter information only.No useful results were obtained for the 2D size of nodules (in cm 2 ) when examining their descriptive statistics and box-plots with acoustic classes.This might be explained by limited interfering between acoustic wavelength and the nodules radii.The high frequency (200 kHz) MBES signal results in ca. 8 mm pulse-wavelength for 1500 m s -1 sound speed in seawater.This wavelength is significantly shorter than the average nodule size in the study area (>3 cm) suggesting that the dominant backscattering is sensitive to nodule density and not to nodule size.Early acoustic studies on Mn-nodules were based on low frequency sonars; therefore there is little or no information about the acoustic backscatter of nodules at high MBES frequencies (> 100 kHz).However, results from this study are in agreement with findings of Weydert (1985) according to which, frequencies higher than 30 kHz are more suitable for mapping the nodule density than the nodule size.This can be attributed to the fact that high frequency signals are more susceptible to surface roughness which is caused by fluctuating nodule densities.Therefore it is suggested that backscatter would increase with increased nodule density given that seafloor roughness increases as more nodules occur per seafloor area.

Implications of acoustic mapping on Mn-nodule resource assessment and benthic habitat characterization
Obtaining high resolution seafloor acoustic classes and quantitative spatial predictions of the Mn-nodule density provides useful information for deep sea mining and impact management.The obvious application is a more realistic resource assessment (total tonnage of Mn-nodules per area) which can assist a better delineation of particular areas with mining interest on large and small scales.Resource assessment can be based on semi-quantitative information provided by acoustic classes that correspond to particular Mn-nodule densities or quantitative results from the RF predictive map.In addition, quantitative maps of Mn-nodule densities can be used to support extrapolations of benthic biota densities to seafloor areas where benthic information is not available.This is possible by considering the nodule substrate as surrogate for habitat mapping of certain biota.Surrogacy for mapping deep sea ecosystems has been incorporated in the study of Anderson et al. (2011); the authors point out, that geomorphic classes can be used for discriminating habitats in broad scales of tens to hundreds of kilometres.They also highlight that any surrogacy approach should be based on the correlation between the physical variables (e.g.bathymetry, backscatter) and the biological patterns that appear in the study area.In Vanreussel et al. (2016) and Amon et al. (2016) it is shown that seafloor covered with more Mn-nodules features higher epifaunal densities.This relation might be further evaluated to have a better and verified relationship between nodule and biota densities allowing estimating biota abundances in larger areas that have only been mapped acoustically.

Conclusions
AUV-based optical and acoustic mapping at high spatial resolution opens up new opportunities for mapping Mn-nodule fields.In this study, automated image analysis provided dense, quantitative information about Mn-nodules at fine scale.This information offers useful insights about the fine scale variability of Mn-nodule densities while it can be utilized for correlations with seafloor acoustic classes and predictive mapping.It was found that the Mn-nodule density within a 500 m x 500 m photo mosaic varies in a pattern of alternating bands (with denser and sparser amounts of nodules) according with smooth bathymetric slopes with a preference of increased nodule occurrence at concave seafloor morphologies.Areas with different nodule densities produced distinct backscatter classes that distinguished nodule populations with distinct mean density values.This suggests that Mn-nodule densities can be efficiently mapped with high resolution hydro-acoustic data.In addition, applying machine learning methodology showed great potential in quantitative predictive mapping of Mnnodules through modelling the complex relation between image-derived nodule metrics with bathymetric derivatives and backscatter statistics.In essence, by using a relatively small amount of AUV images (ca.2700) as the training set it was possible to obtain a 70% correlation between predicted and measured Mn-nodule densities.High quality and spatial resolution AUV hydro-acoustic and optical data can provide a fast and less costly mean for Mn-nodule mapping.This has three major implications in deep sea studies: 1) it raises questions about what causes the Mn-nodules to follow the fine scale bathymetric morphology, 2) it assists in better resource assessment of Mn-nodules and provides the information needed for planning the optimal mining path and 3) it provides more accurate information about Mn-nodule substrate as a benthic habitat, hence it can be utilized for better understanding the deep sea ecology and ecological impact of potential Mn-nodule mining.2).Error sources in quantitative Mn-nodule mapping A few error sources need to be considered when performing seafloor classification and nodule density estimates with optical and acoustic data acquired during multiple AUV deployments.1) Noisy backscatter data: Since the Bayesian approach uses the raw backscatter data, any final classification is susceptible to the effects of noise.Hence, beam incidence angles less than 20 degrees were discarded due to extreme nadir noise effects.The ISODATA classification was based on the backscatter mosaic and its statistics which are also affected mainly by nadir specular noise.It is thus strongly recommended that backscatter data are properly corrected for geometric and sensor-related effects during pre-processing and grids are also filtered/smoothed before the final classification.2) AUV navigation: As exact underwater navigation in 4 km water depth is generally a difficult task, relative misalignments of data from different deployments are very common.Differences in absolute positioning between two deployments can easily amount to 100 m.Thus correlating image based nodule densities from one deployment with backscatter values from another dive might introduce correlation errors that also impact predictability.Although the large scale spatial pattern of classes is well defined, these misalignments can slightly alter the position of class boundaries causing disagreement with the nodule density measurements in places.A correct and verified re-navigation of all AUV-tracks is important for all subsequent analyses.This was done during this study, but slight misalignments remain.3) Nodule sediment blanketing: The effect of Mn-nodules being blanketed by sediment needs to be considered as a source of error here as the individual nodule size and thus the seafloor coverage might be underestimated by automated annotation.Apart from natural sedimentation, the re-deposition of the plume cloud caused by ploughing during the first disturbance experiment (conducted in 1989), has covered certain parts of the nodule field which might lead to a lower nodule densities in those areas.This effect can artificially reduce the correlation between acoustic classes and Mn-nodule densities given that backscatter is not affected by sediment blanketing.

Fig. 1 :Fig. 2 :
Fig. 1: A) The DISCOL area location in the Peru Basin (red star).B) Ship-based, shaded bathymetry of the wider DISCOL area with 40 m pixel size.The black rectangle represents the boundaries of the AUV MBES dataset used in this study (Fig.2).

Fig. 4 :
Fig. 4: A) Points with nodule measurements derived from automated nodule detection, draped on AUV bathymetry, showing Mn-nodules per square meter from perspective view, B) Longitudinal section of bathymetric profile from same area highlighting the local scale morpho-bathymetry of Mnnodule fields.

Fig. 5 :
Fig. 5: A) Bayesian classification map based on AUV backscatter beam data, B) ISODATA classification map based on AUV backscatter neighbourhood statistics (mean, mode, 10 th Q and 90 th Q, see Table2).

Fig. 6 :
Fig. 6: A) Random forests prediction map of Mn-nodules densities, Sensitivity analysis results: B) Percentage of training sample size and performance of RF model in terms of percentage of variance explained (out-of-bag).C) Importance scores of MBES explanatory variables, based on average percentage increase of mean prediction error from ten model runs.

Fig. 7 :
Fig. 7: Inter-comparison of quantitative methods results from the same coverage area (Rectangle made by dense black lines in Fig. 2 A): A) Mn-nodules per image-point (automated nodule-detection from optical images), B) ISODATA classes (10m cell size), C) Bayesian classes (6m cell size), D) RF Mn-nodule density prediction map (6m cell size).

Fig. 8 :
Fig. 8: Box-plots of nodule densities grouped by acoustic class to illustrate the between-class variability.A) Variation of measurements, from samples belonging A) to the same Bayesian classes and B) same ISODATA classes.Blue rectangle bottom and top represent the 25% and 75% percentiles respectively whereas the red line indicates the median value.The whiskers extend to the minimum and maximum value of the samples that are not considered outliers (i.e.: they are no more than ±2.7σ apart).Outliers are marked with red crosses.

Table 2 :
Description of MBES features (bathymetric derivatives and backscatter statistics) that are used as explanatory variables in random forests predictions.

Table 3 :
The nodule-free areas holding lowest backscatter values are captured clearly.Averaged central dB values of the Gaussians derived from reference beam angles on both sides of the AUV MBES

Table 4 :
). RF model performance for minimum optimal settings of training sample and number of trees regarding prediction of Mn-nodule densities.

Table 5 :
Descriptive statistics highlighting the within-class variability of Mn-nodules for both classification methods.