Continuous, High-Resolution Mapping of Coastal Seaﬂoor Sediment Distribution

: Seaﬂoor topography and grain size distribution are pivotal features in marine and coastal environments, able to inﬂuence benthic community structure and ecological processes at many spatial scales. Accordingly, there is a strong interest in multiple research disciplines to obtain seaﬂoor geological and/or habitat maps. The aim of this study was to provide a novel, automatic and simple model to obtain high-resolution seaﬂoor maps, using backscatter and bathymetric multibeam system data. For this purpose, we calibrated a linear regression model relating grain size distribution values, extracted from samples collected in a 16 km 2 area near Bagnoli–Coroglio (southern Italy), against backscatter and depth-derived covariates. The linear model achieved excellent goodness-of-ﬁt and predictive accuracy, yielding detailed, spatially explicit predictions of grain size. We also showed that a ground-truth sample size as large as 40% of that considered in this study was sufﬁcient to calibrate analogous regression models in different areas. Regardless of some limitations (i.e., inability to predict rocky outcrops and/or seagrass meadows), our modeling approach proved to be a ﬂexible tool whose main advantage is the rendering of a continuous map for sediment size, in lieu of categorical mapping approaches which usually report sharp boundaries or rely on a few sediment classes. Conceptualization, S.I., M.I. and R.T.; methodology, S.I., M.I. and M.D.F.; software, S.I., M.I. and G.D.M.; validation, M.D.F.; formal analysis, S.I. and M.I.; investigation, S.I. G.D.M.; resources, M.S. and R.T.; data curation, S.I. M.I.; writing—original S.I. M.I.; writing—review and M.D.F. G.D.M.; visualization, S.I. G.D.M.; M.S.


Introduction
Studying the spatial and temporal distribution of seabed sediments is a pivotal topic in marine and coastal environmental research, ranging from biological and ecological research to engineering applications. It has been shown that seafloor geology, in particular topography and grain size distribution, can influence benthic community structure and ecological processes at many spatial scales (e.g., [1][2][3][4][5][6][7]). At the same time, characterizing grain size distribution is also important in harbour and coastal management [8][9][10], and can be essential when dealing with pollution and environmental remediation [11][12][13][14].
Consequently, there is a strong interest in multiple research disciplines to undertake seafloor geological and/or habitat mapping, and recent studies have presented different modeling approaches for the classification of seabed sediments [15][16][17][18]. Such classification approaches have been performed manually based on the experience of operators [4,19,20], while supervised or unsupervised modeling approaches have recently been proposed [10,[21][22][23][24][25][26]. Manual segmentation and classification can be problematic due to subtle variations that may occur in the large volume of data being collected during modern surveys which might present reproducibility issues between different operators [27]. Supervised approaches, alternatively, rely on statistical methods based on ground-truth data, such as grab samples or dives [28], while unsupervised approaches often rely on some form of statistical clustering [29]. Supervised methods are the preferred approach in marine geological and ecological research [30] but, despite their importance, studying seabed sediments is challenging, given that obtaining representative samples of ground-truth data is usually complex, time-consuming and expensive [31]. Automated approaches relying on remote sensing, along with several statistical tools and dedicated software, can address these issues and are commonly employed nowadays (e.g., [27,[32][33][34][35][36][37][38][39]).
Accordingly, a widespread remote sensing instrument used in studying seafloor sediment characteristics is the multibeam echosounder (MBES), which is currently becoming a fundamental tool for the exploration of both coastal and deep-water environments, although other instruments, such as side-scan sonar, are still widely used [5]. MBES is able to return both bathymetric and acoustic information on seabed characteristics, with applications rapidly increasing in seafloor mapping. MBES backscatter allows the application of quantitative and qualitative analysis to obtain important information for mapping seafloor habitats [40][41][42][43], such as surface roughness [44][45][46][47], seafloor classification [36,[48][49][50][51], substrate type and benthic biota [2,35,52]. Backscatter intensity is related to sediment properties [47,53,54], depending on two main components: volume scattering from sediment inhomogeneities and interface scattering from bottom roughness. The scattering from the sediment volume is created by fluctuations in sediment density or sound velocity. The scattering from the sediment surface is controlled by the impedance difference between the overlying water and the sediments [55,56]. Thus, fine sediments generally exhibit low backscatter intensity due to low sediment bulk density and low acoustic impedance contrast at the water-sediment interface, whereas coarse sediments generally result in higher backscatter strength (e.g., [44,57]). Furthermore, it has been shown that backscatter intensity decreases with mean grain size in sandy sediments [47]. Since the relationship between acoustic backscatter, sediment grain size and benthic community structure has been demonstrated [58][59][60], an accurate grain size classification deriving from the acoustic backscatter becomes a priority from ecological, geological and management perspectives. At the same time, seafloor geology, depth, and topography, along with sediment distribution, are significant drivers of benthic community structure and ecological processes, being widely used in the development of benthic habitat maps [4,5,10,31].
Even though many methods are already available for categorical seabed classification based on MBES data, it has been shown that all of them may have low accuracy [30]. Such lack of accuracy could be due to the classification system, which is often an oversimplification of reality that fits into a rigid scheme failing to detect gradual changes in sediment size and substrate type [30,61]. Some of these issues might be circumvented by means of continuous rather than categorical mapping and, in recent years, continuous modeling approaches have been proposed [18,62,63]. Yet, there is still a need to develop new approaches that might render accurate and detailed benthic habitat maps for the management of coastal environments, which are also those with higher threats from anthropogenic activities [64].
The aim of this study was to provide a novel, automatic and simple model to obtain high-resolution seafloor maps that could provide accurate grain size distribution in coastal environments. For this purpose, we regressed grain size observations gathered from 50 ground-truth samples collected in a 16 km 2 area near Bagnoli-Coroglio (Southern Italy) against backscatter and depth-derived covariates. We compared spatially explicit model predictions in the area with an expert-based seafloor map already available [65]. Additionally, the same model was also projected onto a test area of 18 km 2 , i.e., Lampedusa Island, although here no ground-truth samples were available, except for a few remotely operated vehicle images, along with a former map of the seafloor features based on a mixed expert-based and unsupervised approach [35]. Lastly, we estimated the minimum number of ground-truth samples per km 2 that were needed to recalibrate the model and/or validate the results for new areas.

Data Acquisition
The bathymetric data were collected during two different surveys, i.e., "Abbaco 2017" and "Lampedusa 2015". The Teledyne SeaBat 7125, a multi-frequency MBES set at 400 kHz, providing high-resolution data and collecting backscatter signal for acoustic mosaic production, was used for both surveys [35,65,66]. An SBG EkinoxU inertial system and a Trimble BX982 dual antenna differential global positioning system (DGPS) were used. Sound velocity corrections were provided by a Valeport mini-SVS probe and a Valeport Mini SVP sound velocity profiler. Data were acquired and processed using the PDS 4.1 (Teledyne©inc) software package. Tide data were applied (from https://www. mareografico.it, accessed on 1 February 2022) to the logged files to set up the real depth, and data de-spiking was carried out in order to produce a high-resolution image of the seafloor (chart datum mean sea level, ellipsoid WGS84). The backscatter signal (i.e., snippet data) was processed with an FMGeocoder Toolbox (FMGT) Fledermaus 7.6 version [67]. These data were corrected for receiver gain, transmit power, transmit pulse width, spherical spreading, attenuation in the water column, insonified area, beam pattern, speckle noise, and for angular dependence and local slope [35,65,68].

Bagnoli-Corogolio Calibration Area
The data acquisition of this study area was carried out in 2017 by CNR-ISMAR Naples (formerly CNR-IAMC), within the framework of the ABBaCo Project (Italian project 'Restauro Ambientale e Balneabilità del SIN Bagnoli-Coroglio'-pilot experiments for the environmental restoration and bathing possibility of the Bagnoli-Coroglio coastal area) coordinated by the SZN-Anthon Dohrn, Naples. The project was designed to produce new approaches addressing the remediation of contaminated sites and restoration of marine habitats [66,69]. The study area (Figure 1) was located in the eastern Gulf of Pozzuoli (Naples, southern Italy) and is an integral component of the Campi Flegrei volcanic system, located in the Campania region [66].

Data Acquisition
The bathymetric data were collected during two different surveys, i.e., "Abbaco 2017" and "Lampedusa 2015". The Teledyne SeaBat 7125, a multi-frequency MBES set at 400 kHz, providing high-resolution data and collecting backscatter signal for acoustic mosaic production, was used for both surveys [35,65,66]. An SBG EkinoxU inertial system and a Trimble BX982 dual antenna differential global positioning system (DGPS) were used. Sound velocity corrections were provided by a Valeport mini-SVS probe and a Valeport Mini SVP sound velocity profiler. Data were acquired and processed using the PDS 4.1 (Teledyne©inc) software package. Tide data were applied (from https://www.mareografico.it, accessed on 1 February 2022) to the logged files to set up the real depth, and data de-spiking was carried out in order to produce a high-resolution image of the seafloor (chart datum mean sea level, ellipsoid WGS84). The backscatter signal (i.e., snippet data) was processed with an FMGeocoder Toolbox (FMGT) Fledermaus 7.6 version [67]. These data were corrected for receiver gain, transmit power, transmit pulse width, spherical spreading, attenuation in the water column, insonified area, beam pattern, speckle noise, and for angular dependence and local slope [35,65,68].

Bagnoli-Corogolio Calibration Area
The data acquisition of this study area was carried out in 2017 by CNR-ISMAR Naples (formerly CNR-IAMC), within the framework of the ABBaCo Project (Italian project 'Restauro Ambientale e Balneabilità del SIN Bagnoli-Coroglio'-pilot experiments for the environmental restoration and bathing possibility of the Bagnoli-Coroglio coastal area) coordinated by the SZN-Anthon Dohrn, Naples. The project was designed to produce new approaches addressing the remediation of contaminated sites and restoration of marine habitats [66,69]. The study area ( Figure 1) was located in the eastern Gulf of Pozzuoli (Naples, southern Italy) and is an integral component of the Campi Flegrei volcanic system, located in the Campania region [66].  The marine area of the Bagnoli coastal zone is part of the caldera collapse that originated from the eruption of the Neapolitan Yellow tuff [70]. The infralittoral setting is characterized by relatively coarse-grained sediments, mainly sand and silty-sand, representing the grain size distribution in equilibrium with a shallow-water, high-energy environment. The sediment constituents are mainly formed by bioclasts and volcaniclasts with a sandy-silt matrix that becomes more abundant towards the distal sector of the gulf [65]. The Bagnoli area has undergone intense industrial activity during most of the last century, from the 1910s until the 1990s. Human-related activities altered the original morphology of the coast, e.g., for the construction of infrastructures on both land and at sea. Steel production ultimately ended in 1990, and the industrial facilities were completely dismantled in the early 2000s [70]. Since then, a series of surveys were carried out with the aim of providing a characterization of soil on land and marine sediments offshore, addressed to the environmental recovery of the area (e.g., [11,13,71,72]). With the ABBaCo project, new results were obtained aimed at the acquisition of new and updated data [8,65,66,[73][74][75][76]. Similarly, a high-resolution shaded image derived from bathymetric data ( Figure 1) and acoustic mosaic from backscatter data ( Figure 2) were obtained, covering about 16 km 2 in a water depth range of 1.5-115 m bsl. The marine area of the Bagnoli coastal zone is part of the caldera collapse that originated from the eruption of the Neapolitan Yellow tuff [70]. The infralittoral setting is characterized by relatively coarse-grained sediments, mainly sand and silty-sand, representing the grain size distribution in equilibrium with a shallow-water, high-energy environment. The sediment constituents are mainly formed by bioclasts and volcaniclasts with a sandy-silt matrix that becomes more abundant towards the distal sector of the gulf [65]. The Bagnoli area has undergone intense industrial activity during most of the last century, from the 1910s until the 1990s. Human-related activities altered the original morphology of the coast, e.g., for the construction of infrastructures on both land and at sea. Steel production ultimately ended in 1990, and the industrial facilities were completely dismantled in the early 2000s [70]. Since then, a series of surveys were carried out with the aim of providing a characterization of soil on land and marine sediments offshore, addressed to the environmental recovery of the area (e.g., [11,13,71,72]). With the ABBaCo project, new results were obtained aimed at the acquisition of new and updated data [8,65,66,[73][74][75][76]. Similarly, a high-resolution shaded image derived from bathymetric data ( Figure 1) and acoustic mosaic from backscatter data ( Figure 2) were obtained, covering about 16 km 2 in a water depth range of 1.5-115 m bsl. The inner shelf is characterized by the presence of stepped terraced surfaces located, in the western sector, at water depths of 10, 25, and 35 m, S-SW oriented. These terraces are expressions of erosional features or toplap/offlap surfaces of local prograding wedges and are likely the result of the dynamic equilibrium between seafloor erosion at the base of storm waves and the sediment supply from the coastline during the Holocene sea-level rise [77]. Stepped terraced surfaces, mainly W oriented, are also present around Nisida The inner shelf is characterized by the presence of stepped terraced surfaces located, in the western sector, at water depths of 10, 25, and 35 m, S-SW oriented. These terraces are expressions of erosional features or toplap/offlap surfaces of local prograding wedges and are likely the result of the dynamic equilibrium between seafloor erosion at the base of storm waves and the sediment supply from the coastline during the Holocene sea-level rise [77]. Stepped terraced surfaces, mainly W oriented, are also present around Nisida island, where the morphology displays a rugged seafloor with a gentle slope (<1 • ) within 30 m of depth. Two slope breaks bound the island (~10 • at 30 m depth and~15 • at 45 m), connecting the inner continental shelf to the deeper part of the bay (see Figure 4 in [66]). During the ABBaCo project, 50 grab samples were collected and processed for sedimentological and grain-size analysis. The gravel/sand fraction (4000−0.063 mm) was analysed using dry sieving at phi interval, while the fine-grained fraction (0.063−0.0040 mm) was processed with a SYMPATEC laser particle size analyzer [74].

Lampedusa Test Area
Lampedusa is the largest island of the Pelagie Archipelago, followed by Linosa island and the Lampione islet (Sicily, Italy; Figure 3). The island is located in the southern Mediterranean Sea, inside the Sicily Channel, arising from the northern edge of the African continental plate [78,79]. island, where the morphology displays a rugged seafloor with a gentle slope (<1°) within 30 m of depth. Two slope breaks bound the island (~10° at 30 m depth and ~15° at 45 m), connecting the inner continental shelf to the deeper part of the bay (see Figure 4 in [66]). During the ABBaCo project, 50 grab samples were collected and processed for sedimentological and grain-size analysis. The gravel/sand fraction (4000−0.063 mm) was analysed using dry sieving at phi interval, while the fine-grained fraction (0.063−0.0040 mm) was processed with a SYMPATEC laser particle size analyzer [74].

Lampedusa Test Area
Lampedusa is the largest island of the Pelagie Archipelago, followed by Linosa island and the Lampione islet (Sicily, Italy; Figure 3). The island is located in the southern Mediterranean Sea, inside the Sicily Channel, arising from the northern edge of the African continental plate [78,79]. The island (surface area 20 km 2 ; maximum elevation 133 m asl) shows a sub-planar surface inclined toward the south-west and is carved by deep valleys that connect with the narrow inlets along the coast. The island displays rocky cliffs up to 120 m high, which outline a rugged coastline, with tens of rocky highs at a short distance from the coast, mostly occurring along the north-western sector [79,80]. The Lampedusa shoreline can be divided into two main sectors: in the north (from Capo Ponente to Capo Grecale) they are indented, with dominant coastal features varying from steep to sub-vertical cliffs, along with small promontories and bays. In the south, the coast gradually slopes down and includes pocket beaches in the coves. A geophysical survey of Lampedusa was carried out in 2015 by CNR-ISMAR Napoli (formerly CNR-IAMC), within the framework of the project "CAmBiA-Contabilità Ambientale e Bilancio Ambientale" to assess the conservation status and map the distribution of Posidonia oceanica (L.) Delile (Potamogetonaceae) meadows. The aim of the survey was the characterization of the seabed of Lampedusa Island with MBES (both for bathymetric and backscatter data) along with video-camera inspections as ground-truth data (see Figure 3 in [79]) to assess the distribution and the conservation of the P. oceanica meadows. The high-resolution bathymetric map of Lampedusa The island (surface area 20 km 2 ; maximum elevation 133 m asl) shows a sub-planar surface inclined toward the south-west and is carved by deep valleys that connect with the narrow inlets along the coast. The island displays rocky cliffs up to 120 m high, which outline a rugged coastline, with tens of rocky highs at a short distance from the coast, mostly occurring along the north-western sector [79,80]. The Lampedusa shoreline can be divided into two main sectors: in the north (from Capo Ponente to Capo Grecale) they are indented, with dominant coastal features varying from steep to sub-vertical cliffs, along with small promontories and bays. In the south, the coast gradually slopes down and includes pocket beaches in the coves. A geophysical survey of Lampedusa was carried out in 2015 by CNR-ISMAR Napoli (formerly CNR-IAMC), within the framework of the project "CAmBiA-Contabilità Ambientale e Bilancio Ambientale" to assess the conservation status and map the distribution of Posidonia oceanica (L.) Delile (Potamogetonaceae) meadows. The aim of the survey was the characterization of the seabed of Lampedusa Island with MBES (both for bathymetric and backscatter data) along with video-camera inspections as ground-truth data (see Figure 3 in [79]) to assess the distribution and the conservation of the P. oceanica meadows. The high-resolution bathymetric map of Lampedusa (depth range of 2-50 m bsl) showed a rugged seafloor all around the island within the first 20-40 m depth, due to extensive rocky outcrops, localized talus deposits and relict morphologies of the seabed, such as P. oceanica meadows on 'matte' facies (as defined by [31,79,81] vertical scarps 10-20 m high are present, such as in the eastern and northern shallow water sectors. Along the northern coast, the seabed is steeper and dips down to over 50 m bsl (see [79]). Conversely, along the southern sector of the island, the seabed slopes with gradually decreasing gradient in the depth range of 10-50 m. As in Bagnoli, in this study area the snippet data of the MBES was also processed, obtaining an acoustic mosaic at 2.5 m resolution ( Figure 4) that was used in [35] for automatic image classification.
(depth range of 2-50 m bsl) showed a rugged seafloor all around the island within the first 20-40 m depth, due to extensive rocky outcrops, localized talus deposits and relict morphologies of the seabed, such as P. oceanica meadows on 'matte' facies (as defined by [31,79,81]). Locally, vertical scarps 10-20 m high are present, such as in the eastern and northern shallow water sectors. Along the northern coast, the seabed is steeper and dips down to over 50 m bsl (see [79]). Conversely, along the southern sector of the island, the seabed slopes with gradually decreasing gradient in the depth range of 10-50 m. As in Bagnoli, in this study area the snippet data of the MBES was also processed, obtaining an acoustic mosaic at 2.5 m resolution (Figure 4) that was used in [35] for automatic image classification.

Modeling and Mapping Approach
We used the granulometric distribution values derived from the 50 ground-truth samples collected in Bagnoli to calibrate our model. In detail, the different percentages of sediments belonging to each particle size class (from the coarse to the finest) were used to calculate the mean value of the phi scale using the 'gran.stats' function of the 'rysgran' package [82]. Phi is a logarithmic parameter for grain size classification whose values range from −9 (i.e., boulders) to 9 (i.e., clay) and is calculated as: The package 'rysgran' is based on the size scales adopted by [83,84] while the phi scale was introduced by [85].
We used the values of phi as the dependent variable in our modeling approach using backscatter and depth data as independent variables. In order to take into account the effects of seafloor morphology, slope, aspect, terrain ruggedness index (i.e., TRI), the topographic position index (i.e., TPI), and roughness were computed from the depth data using the function 'terrain' from the package 'raster' [86].
The set of independent variables was checked for multicollinearity according to a variance inflation factor (i.e., VIF) using the function 'vifstep' in the package 'usdm' [87]. Specifically, we dropped all variables with VIF ≥ 5 as a threshold indicator of multicollinearity among predictors [88].

Modeling and Mapping Approach
We used the granulometric distribution values derived from the 50 ground-truth samples collected in Bagnoli to calibrate our model. In detail, the different percentages of sediments belonging to each particle size class (from the coarse to the finest) were used to calculate the mean value of the phi scale using the 'gran.stats' function of the 'rysgran' package [82]. Phi is a logarithmic parameter for grain size classification whose values range from −9 (i.e., boulders) to 9 (i.e., clay) and is calculated as: The package 'rysgran' is based on the size scales adopted by [83,84] while the phi scale was introduced by [85].
We used the values of phi as the dependent variable in our modeling approach using backscatter and depth data as independent variables. In order to take into account the effects of seafloor morphology, slope, aspect, terrain ruggedness index (i.e., TRI), the topographic position index (i.e., TPI), and roughness were computed from the depth data using the function 'terrain' from the package 'raster' [86].
The set of independent variables was checked for multicollinearity according to a variance inflation factor (i.e., VIF) using the function 'vifstep' in the package 'usdm' [87]. Specifically, we dropped all variables with VIF ≥ 5 as a threshold indicator of multicollinearity among predictors [88].
The remaining predictors were used in a linear regression model using phi as the dependent variable. We used a backward model selection using function 'drop1' in package 'stats' [89], recursively dropping all non-significant predictors. The remaining predictors were further sub-selected according to the Akaike information criterion (i.e., AIC) using Remote Sens. 2022, 14, 1268 7 of 18 the function 'stepAIC' in the package 'MASS' [90]. The AIC is routinely adopted in model selection applications, as the optimal model can be identified according to the lowest value of AIC [91]. Final models' goodness-of-fit was evaluated calculating the adjusted coefficient of determination (R 2 adj .). Moreover, predictive performance of the best models was assessed through a 10-fold cross-validation approach, which was repeated five times using the function 'train' in the package 'caret' [92]. At each cross-validation run, we calculated the values of the predictive coefficient of determination (R 2 pred .) and rootmean-square error (RMSE), then averaging the results. The residuals from the best model according to the aforementioned criteria were inspected to assess possible violations of the linear regression assumptions. We also took into consideration the possible presence of spatial autocorrelation in the model residuals, which we tested using the function 'correlog' in the package 'ncf' [93], using 5 m increments for the uniformly distributed distance classes and 1000 permutations to assess the level of significance.
The model was used to produce a 2.5 m resolution map for phi values over the Bagnoli area and was also projected onto a test area (i.e., Lampedusa) where ground-truth data were not available. The final maps were produced by projecting the best regression model on both Bagnoli and Lampedusa areas, also filtering all those pixels reporting phi values >9 and <−9, which could result from minor anomalies in either the backscatter or the depth data. In order to render a more homogenous final map, we computed mean values for each pixel using a moving window approach by a 5 × 5 mean filter, using thefunction 'focal' in the package 'raster' [86].
In both study areas, we evaluated whether the models were extrapolated beyond the calibration range by means of multivariate environmental similarity surface (i.e., MESS), using the function 'mess' in the package 'dismo' [94]. In detail, we highlighted in the final output all those areas which had MESS ≤ −10, which can be considered areas of strict extrapolation [95].
Finally, we performed a sensitivity analysis to establish how many ground-truth samples would be needed in order to calibrate and/or validate the model outputs when applied to other areas beside the ones we studied. In order to do so, we took random subsamples from our original dataset, each with increasing sample size compared to the original one (i.e., from 10% to 90% of our original 50 ground-truth data). For each subsample, we recalibrated the model and compared the coefficients to the full model. Each subsample was randomly repeated 50 times. All statistical analyses were performed using R 4.1.2 [89].

Ground-Truth Grain Size Classification from Bagnoli-Coroglio
The results from the grain size classification analyses are presented in supplementary material Table S1. The majority of the samples could be generally classified as sand (76%), ranging from very coarse sand (i.e., CS, 2%) to very fine sand (i.e., VFS, 14%). The remainder of the samples (24%) could be classified as silt, ranging from coarse silt (i.e., CSi, 10%) to fine silt (i.e., FSi, 6%). In terms of phi, the computed mean values showed an average of 2.67 ± 1.80 standard deviations. The whole distribution of phi values used for the model calibration can be seen in Figure 5.

Model Calibration
After checking for multicollinearity among independent variables, we were left with a set of five covariates, i.e., depth, backscatter, aspect, TPI, and roughness. Starting from a full model, we recursively dropped all non-significant predictors until we were left with only depth and backscatter as significant covariates. In addition, these two variables were included in the best model after backward AIC selection. The final model achieved a high goodness-of-fit, reporting a R 2 adj .= 0.858, also showing excellent predictive performances (i.e., R 2 pred .= 0.865 and RMSE = 0.647, Table 1). The statistical analyses showed that both depth and backscatter were excellent predictors for phi, albeit backscatter remarkably Remote Sens. 2022, 14, 1268 8 of 18 outperformed depth when the covariates were considered singularly ( Table 1). The final model showed no deviation from linear regression assumptions in the residuals. Spatial autocorrelation was negligible, with a mean value of Moran's I of −0.027 and only 10% of 633 distance classes results significant (p < 0.05).

Final Maps for Bagnoli and Lampedusa
The final maps for Bagnoli and Lampedusa are shown in Figures 6 and 7, respectively. For Bagnoli, only 0.33% of the total pixels in the image could be considered as strictly extrapolated according to MESS analysis, while in Lampedusa the corresponding figure was 1.28%. The mean value of phi for the overall predicted area of Bagnoli was 3.57 ± 1.91, while for Lampedusa the value was 2.84 ± 1.17.

Final Maps for Bagnoli and Lampedusa
The final maps for Bagnoli and Lampedusa are shown in Figures 6 and 7, respectively. For Bagnoli, only 0.33% of the total pixels in the image could be considered as strictly extrapolated according to MESS analysis, while in Lampedusa the corresponding figure was 1.28%. The mean value of phi for the overall predicted area of Bagnoli was 3.57 ± 1.91, while for Lampedusa the value was 2.84 ± 1.17.

Final Maps for Bagnoli and Lampedusa
The final maps for Bagnoli and Lampedusa are shown in Figures 6 and 7, respectively. For Bagnoli, only 0.33% of the total pixels in the image could be considered as strictly extrapolated according to MESS analysis, while in Lampedusa the corresponding figure was 1.28%. The mean value of phi for the overall predicted area of Bagnoli was 3.57 ± 1.91, while for Lampedusa the value was 2.84 ± 1.17.  The results from the sensitivity analysis can be seen in Figure 8. The results showed that, even with a sample size as low as 40% of our original dataset, all model parameters (i.e., intercept and coefficients for depth and backscatter) converged on the results of the full model, albeit with some variability resulting from the 50 random replicates with each subsample. Thus, considering that the area of Bagnoli was 16 km 2 and that we took 50 ground-truth samples, we could have reached similar results with only 20 ground-truth samples. Consequently, 1.3 ground-truth samples per km 2 would be a viable trade-off for validating and/or calibrating our modeling approach in other areas. Accordingly, for the area of Lampedusa, for instance, the results from our model could be validated by 24 ground-truth samples, as the overall surveyed area was approximately 18 km 2 . Figure 8. Results from the sensitivity analysis. Each point represents the mean and the standard deviation of 50 replications of random subsamples of the original dataset, each with increasing sample size compared to the original one (i.e., from 10% to 90% of the original 50 ground-truth data). For each subsample, the model was recalibrated and compared to the intercept and coefficients of the full model.

Bagnoli-Coroglio Calibration Area
In the Bagnoli-Coroglio area, the seafloor map showed good correspondence with the sedimentological map preliminarily recognized by expert-based interpretation of the backscatter mosaic published in [65]. The coarse-grained sediments were concentrated in two sectors: the area between Pozzuoli and Bagnoli and the southern sector of the Bagnoli-Coroglio offshore, including Nisida island and the area at its south. In particular, the northern sector was characterized by relatively high acoustic reflection (Figure 2) that corresponded to "gravelly coarse sand" in [65]. These sediments were mostly concentrated in two areas, one below a slide deposit that contributes a supply of coarse material from the coast, and the other corresponding to a sub-circular depression about 0.1 km 2 wide where coarse sand settled. Such depression was described in more detail by [8]. Figure 6 showed the same gravelly coarse sand distribution (blue area), both from north and south sectors, and also showed coarse sand entrapped by bottom currents near the piers of the disused industrial site of Bagnoli [8].
In the southern part of the first study area, offshore of the Coroglio tuffaceous cliff, a series of morphologic depressions and ridges have favored a current-driven separation between the coarse deposits (depressions) from the fine-grained material (ridges), as shown by [8]. These morphological structures are clearly visible in Figure 6, supporting the outputs of the applied method. The central part of the study area, off the pier of the former industrial site, was classified as fine sand/very fine sand in [65], based on both backscatter values and seafloor sampling. Nevertheless, the acoustic mosaic showed some patches of intermediate to high backscatter, likely due to the occurrence of coarser shallow water deposits, which would be hard to discriminate visually by an expert. In this specific case, the seafloor map differs from that of [65], because the modeling approach presented here was able to recognize the heterogeneous distribution of medium sand to fine sand, with a predominance of the latter. In contrast, expert-based classification, or even automatic image classifications, might coerce aggregations or render sharp boundaries between the facies, preventing the realization of a more realistic distribution between the granulometric classes. Thus, even when close to reality, polygon classification can be limited in terms of resolution and detail compared to the outputs produced in this study. The same can be said for silty sediments, as this classification creates the possibility of distinguishing between coarse silt to very fine silt, avoiding creating a single class of silt and clay, as was done in [65]. Silt and clay were detected in the south-western distal part of the study area corresponding to the Bagnoli Valley [77], a morphostructural depression connecting the Bagnoli foreshore with the outer continental shelf of Pozzuoli Bay, which delivers coastal sediments towards deep water areas. In [66], a flow analysis, based on seafloor morphology data, was carried out to derive the stream paths towards the Valley, that showed a preferential direction of the sediment transport that fits with what is shown in Figure 6.
Finally, the presence of bedrock in the southern sector of Nisida island was recognized by experts from the morphology data published in [65]. In this sector, linear regression classified this sub-outcropping structure as an alternation of fine sand and coarse sand. Admittedly, as there were no grabs done in this sector and considering that the model was calibrated on loose sediments, rocky outcrops or very large boulders could have been misidentified.
Yet, the MESS sectors of the Bagnoli study area were only 0.33% of the total pixels and, analysing the map, they corresponded mostly to backscatter noise or to coarser sediments, which were underrepresented in the calibration data and might have led to strict extrapolation in some minor parts of the map.

Lampedusa Island Test Area
The sediment distribution map obtained from the linear regression model for Lampedusa, where no ground-truth was available, is in concordance with the classification carried out in [35]. In the latter study, a preliminary benthoscape classification of the island was carried out with the use of remote sensing object-based image analysis (RSOBIA), that segmented and classified the acoustic mosaic obtained from snippet data processing [96]. The map in [35] was presented as preliminary because no grab sampling was available and the seabed was mainly classified only by experts based on its acoustics facies pattern (i.e., fine sediments exhibiting low backscatter, and coarse sediments corresponding to high backscatter), and the results of RSOBIA segmentation. In general, the seafloor map produced here by the model was characterized mostly by fine sand, representing the predominant grain size. In detail, the southern sector of the island was mostly characterized by a distribution of fine and very fine sand, with channels of medium coarse sand. Posidonia oceanica meadows were identified with video-inspections and added in benthoscape classification as overprint layers [35,79]. Well-developed meadows of this plant cause a decrease in the backscatter values, shifting the classification toward fine sands, as already shown in [31]. So, in the case of meadows or any other plant-dominated seafloors, visual inspections or scuba diving as ground-truth information can outperform the linear modeling approach. The sediment distribution of the southern flank showed a classic depositional model of coastal environments, with the onshore coarse grain-size (blue areas in Figure 7) and finer material offshore. Coarser sand was present along the coastline and in all north-western sectors of Lampedusa, i.e., from Capo Ponente to P. Muro Vecchio. As shown by ROV images (see Video 1 in Figure 4 published in [35,79]), very coarse heterogeneous materials (well-rounded pebbles and boulders) in poorly sorted calcareous gravels were detected, forming a terraced surface at the foot of the submerged cliffs. This formation was interpreted as a foot-cliff deposit generated by the wave-related winnowing of talus, likely derived by progressive cliff failures [97]. As already showed for the Bagnoli-Coroglio site, rocky outcrops can be misclassified by the linear regression model, due to underrepresentation in the calibration ground-truth data.
Accordingly, the MESS areas were more present in sectors with higher reflectivity, as for Bagnoli, because there were no calibration samples for this class. The same reasoning can be followed for very fine sediments, such as clay without silt (e.g., Cala Creta in Figures 4 and 7), which were also underrepresented in the calibration data from Bagnoli-Coroglio, where the finest sediments were very fine silt grouped with clay. Admittedly, the linear model predictions in the Lampedusa area need to be validated by grab samples that, according to our sensitivity analysis, could be as few as 24 samples.

Comparisons with Other Approaches
Supervised methods for seafloor grain size and/or habitat classification based on remote sensing data have been shown to provide accurate results, especially when paired with robust ground-truth sampling [21,27,98]. Yet, several supervised approaches resulted in comparatively low accuracy and high error rates when dealing with a morphologically simple seabed characterized by a high sedimentary complexity [30]. Despite state-of-theart statistical approaches, such low map accuracy could be attributed to an inadequate classification system, low discriminatory power of the available predictors, and the spatial complexity of the survey site compared to the positioning accuracy of the grab samples [30]. Our approach can be seen as different and innovative compared to other studies, as it used a quantitative model for the spatial distribution of sediment fractions, which did not oversimplify the final output into a few categories. On the contrary, as the final model outputs a map for phi values, it can be seen as an extremely detailed map that clearly shows gradients and boundaries among sediments, allowing for more precise habitat mapping or geological study of the seafloor in coastal environments. Such a detailed map can be easily transformed into a simpler map of substrate types, as for instance those defined by the European Nature Information System (EUNIS) hierarchical classification [99] or other marine habitat classification schemes [61] After cross-validation, our model showed excellent predictive performance (R 2 pred. = 0.865). Such a level of accuracy exceeds an 85% target accuracy that has usually been considered ideal for thematic mapping based on remotely sensed data [100]. Such a level of accuracy can be surpassed by other modeling approaches, such as fully connected conditional random field [101], but this model was calibrated on a few substrate classes rather than a continuous variable, such as phi.
Another typical classification model for seafloor data based on MBES is random forest [27,[102][103][104]. This type of model has been shown to offer accuracy as high as 83% [103]. Moreover, when applied to seafloor data, random forest is used as a classification approach for few categories, such as those proposed by EUNIS, showing limitations and underperformance when substrate complexity is present [30]. Admittedly, random forest regression for continuous data is widely used in ecological research [105,106], yet, in continental shelf research, random forest is mostly used as a supervised classification approach, although there are recent examples where it has been employed for continuous modeling [18,62,63].
Other statistical methods rely on the use of proprietary software compared to free software, such as R or QGIS. Some approaches rely on object-based image analysis, for example, RSOBIA, which is implemented in ArcGIS [30,35,96,107]. The segmentation implemented in RSOBIA requires the user to define both the number and the minimum size of the clusters, which are then classified by a k-means clustering approach. Another ArcGIS tool is Terrain Attribute Selection for Spatial Ecology [108], which derives six variables from bathymetry that can be used along with backscatter data in a multiresolution segmentation algorithm requiring expert judgment to define scale and input variables. Both the aforementioned approaches in ArcGIS can offer good predictive performance, but they require proprietary software, some degree of expert handling for fine tuning the model and, finally, they both suffer from the same problems as random forest when dealing with complex seafloors [30].
Other approaches to seafloor mapping have been applied as well, including both expert-based and supervised classification. The authors of [109] focused on soft-sediment physical features to delineate distinct acoustic facies that presumably represented different lithological types of seabed sediments. The different facies were defined in ArcGIS according to similar backscatter characteristics, which were later validated by sediment samples analysed for particle size.
In [110], the applicability of supervised learning techniques for benthic habitat characterization using angular backscatter response was demonstrated. The study involved characterization of acoustic backscatter data from multibeam systems, using four different supervised learning methods to generate benthic habitat maps.
In [5], the authors combined a suite of geomorphometric and textural analytical techniques to map specific types of seafloor morphologies and compositions. The seabed was classified into five elementary morphological zones, using morphometric derivatives, bathymetric position index and geomorphometric mapping, to produce morphologic and seabed composition maps of the predominant habitats. Finally, ref. [111] used ensemble modeling to map seafloor sediment distributions and shifts with minimal ground-truth data combined with hydroacoustic datasets in the Sylt outer reef.

Limits and Future Developments
Despite several, advantages, our modeling approach was unable to correctly identify rocky outcrops and P. oceanica meadows. The model was also constrained by the calibration data, which was underrepresented in both very fine and very coarse sediments. Moreover, as our survey areas were both limited to shallow water, the model needs to be recalibrated when waters are deeper than 110 m. In its current stage of development, our model is best suited to all those shallow-water areas with loose sediments intermingled with a few hard substrates, regardless of terrain, that can be either smooth or rough [43]. Future developments of the model will be able to detect the presence of consolidated sediments and function in deeper waters by expanding the calibration dataset.

Conclusions
The modeling approach presented in this research can be seen as an accurate method for high resolution and continuous mapping of loose sediments in the continental shelf, which can be used in several applications from habitat mapping to substrate shifting when the analysis is repeated across time. The method relies on free software and can be validated and/or recalibrated using 1.3 ground-truth samples per km 2 of survey area. The results from the Bagnoli-Coroglio area showed excellent accuracy and negligible areas of extrapolation, while the map for Lampedusa was largely in concordance with expert-based previous seafloor mapping [35,79], although there is a need to validate the results by grab samples. Regardless of these limitations, our modeling approach is a flexible tool whose main advantage is the rendering of a continuous map for sediment size that does not show sharp boundaries or few sediment classes. These features can be pivotal in ecological, geological and management applications in marine coastal ecosystems.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/rs14051268/s1, Table S1: Results from the granumoletric analyses from the 50 grab samples in the Bagnoli-Coroglio calibration area. The results include the percentages from all granulometric classes (according to Udden & Wentworth) and the corresponding phi values (according to Krumbein). Data Availability Statement: The entire data for the grab samples is available as Supplementary Materials. The raster files for backscatter and depth for both the calibration and test areas can be requested at Sara Innangi sara.innangi@cnr.it.