**Assessing Fine-Scale Distribution and Volume of Mediterranean Algal Reefs through Terrain Analysis of Multibeam Bathymetric Data. A Case Study in the Southern Adriatic Continental Shelf**

**Fabio Marchese 1,2,\*, Valentina Alice Bracchi 1,2, Giulia Lisi 1, Daniela Basso 1,2, Cesare Corselli 1,2 and Alessandra Savini 1,2,3**


Received: 21 November 2019; Accepted: 2 January 2020; Published: 4 January 2020

**Abstract:** In the Mediterranean Sea, crustose coralline algae form endemic algal reefs known as Coralligenous (C) build-ups. The high degree of complexity that C can reach through time creates notable environmental heterogeneity making C a major hotspot of biodiversity for the Mediterranean basin. C build-up can variably modify the submarine environment by affecting the evolution of submerged landforms, although its role is still far from being systematically defined. Our work proposes a new, ad-hoc semi-automated, GIS-based methodology to map the distribution of C build-ups in shallow coastal waters using high-resolution bathymetric data, collected on a sector of the southern Apulian continental shelf (Southern Adriatic Sea, Italy). Our results quantitatively define the 3D distribution of C in terms of area and volume, estimating more than 103,000 build-ups, covering an area of roughly 305,200 m2, for a total volume of 315,700 m3. Our work firstly combines acoustic survey techniques and geomorphometric analysis to develop innovative approaches for eco-geomorphological studies. The obtained results can contribute to a better definition of the ocean carbon budget, and to the monitoring of local anthropogenic impacts (e.g., bottom trawling damage) and global changes, like ocean warming and acidification. These can affect the structural complexity and total volume of carbonate deposits characterizing the Mediterranean benthic environment.

**Keywords:** geomorphometry; seafloor mapping; spatial analysis; algal reefs; bioconstruction volume; multibeam bathymetry

#### **1. Introduction**

Calcareous algal reefs are 3D biogenic-carbonate structures mostly built by calcareous red algae (CCA) that typify the seascape, from roughly 10 m of water depth down to those areas where reef-forming algae can thrive in dim light conditions [1–4]. In the warmer oceans, CCA may contribute to the complex of *coral reefs* or, on a small scale, may form independent bioherms. In temperate waters, like the Mediterranean Sea, CCA form endemic algal reefs known as Coralligenous (C) buildups, which represent the second most important hot-spot of biodiversity in the Mediterranean Sea after *Posidonia oceanica* meadows [4]. C is able to produce large deposits of biogenic calcium carbonate [5] according to climate change, oceanographic conditions, accommodation space, terrigenous supply variations and human-caused impact [6–10].

C geomorphological expression plays a crucial key-role in shaping the Mediterranean seascape and affecting its evolution through geological times, as it constitutes the most massive reefs along the Mediterranean continental shelf [11]. Nevertheless, a precise distribution of C is still uncomplete [11–14]. C morphotypes have been categorized mainly in two groups, banks and rims [4,15,16] that is based upon the nature of the substrate. Among banks, several different terms have been used to distinguish diverse local morphotypes (heads, blocks, patches, or banks [17]; vertical pillar, [18]; algal reefs, [19,20]; mound-shaped algal banks, [21]; minute reef aggregation and superficial layer formations, [22]; columns and ridges, [14]). A first attempt to use shape geometry descriptor to extract C from backscatter data was proposed by [23], but only recently a new categorization of C morphotypes, based on such a quantitative technique, has been proposed [11]. Moreover, C is responsible for the production and storage of a large amount of carbonate, but its contribution to the ocean carbon budget is still poorly quantified.

We apply geographic information system (GIS)-based analysis on digital elevation models (DEMs), obtained from high-resolution bathymetry data collected by multibeam echosounder (MBES) to investigate the variety of seaforms created by C build-ups, reaching scales and resolutions that are comparable to subaerial geomorphology studies [24].

We propose a new method for C feature extraction, designed upon the principals of geomorphometry and using established algorithms for surface analysis [25,26], applied to a well-studied C area located offshore Puglia Adriatic coasts (Italy) where a diverse type of C build-ups occur [11,14]. Our work focused on: (1) the development of a technique for automated identification and extraction of C build-ups; (2) the characterization of the spatial and morphological distribution of C; (3) the quantification of the total volume of selected C build-ups in relation to the field's subsurface geology, as determined from sub-bottom profiler (SBP) data.

Exploring and imaging submarine landforms requires the development of innovative approaches in using methods and techniques originally developed for the terrestrial environment (i.e., geomorphometry). In particular, we focus on providing quantitative information on C habitat as carbonate deposits on the Mediterranean shelf, which can be used as a baseline to which other C areas may be compared quantitatively.

#### **2. Data and Methods**

#### *2.1. Data Acquisition and Processing*

High-resolution acoustics data ground-truthed by video records were obtained by ship-based research surveys performed within the framework of BioMAP project (BIOcostruzioni Marine in Puglia, -P.O. FESR 2007/2013) that promoted actions in order to map and monitor C habitats along the Apulian shelf (southern Adriatic margin and northern Ionian margin, Mediterranean Sea) [27]. The survey has been carried out using a small research boat (Calafuria ISSEL, property of CoNISMa) in March and April 2013, acquiring 16 nautical miles of MBES survey lines, covering the seafloor over an area of 1.7 km2, at depths ranging from 5 to 40 m and using a Teledyne RESON SeaBat 8125 (455 kHz) MBES (Figure 1). Teledyne PDS2000TM software was used to acquire and process bathymetric data producing very high-resolution DEMs (0.5 m × 0.5 m). The correct position of the acquired data has been assured by a Hemisphere Crescent R-Series dGPS, and an IXSEA OCTANS motion sensor and gyro. Water sound velocity was obtained by a Valeport MiniSVS sound velocity sensor and by a Teledyne Reson SVP15 sound velocity profiler.

A set of high-resolution SBP lines have been also acquired on board the same research boat in an area NW of the Bari port (Figure 1), using a parametric Innomar 2000 (SES 2000), with an operating frequency of 10 kHz and 5 kHz. SBP data were processed using ISE Post-Processing Software, converted into SGY format and analyzed using PETREL (Schlumberger).

**Figure 1.** Location of the study area, bathymetric model for inset from [28]. Shaded relief map of the high-resolution multibeam bathymetry data acquired over the study area. Red track lines indicate the SBP survey and the trawled camera route. A, B, C, D, E, F, G areas are magnifications for the other figures in this paper.

The study area has also been surveyed by video inspections to ground-truth remote dataset, using an underwater trawled camera (Quasi-Stellar © Elettronica Enne) (Figure 2).

**Figure 2.** Examples of underwater images of Coralligenous build-ups made by the trawled camera. Following the classification proposed by [11]: (**A**) Discrete Relief, (**B**) Hybrid Bank and (**C**) Tabular Banks.

#### *2.2. Geomorphometric Analysis and Extraction of C Build-ups*

The entire spatial dataset (DEMs, survey routes and video tracks) was integrated into ArcGIS™ and ad-hoc geomorphometric analysis was performed with SAGA (System for Automated Geoscientific Analysis [29]). In particular, basic terrain parameters were extracted from selected areas (Figure 3) and the Topographic Position Index (TPI) proposed by [30] has been performed at the finest possible scale (1 inner radius and 5 outer radii), according to the DEM resolution (i.e., 0.5 × 0.5 m grid cell size) [30–32]. TPI is actually only one of a vast array of morphometric parameters based on neighboring areas that can be useful in DEM analysis [32]. In particular, positive TPI values well represent locations that

are higher than the average of their surroundings, as defined by the neighborhood (ridges), whereas negative TPI values represent locations that are lower than their surroundings (valleys). TPI values near zero are either flat areas (where the slope is near zero) or areas of constant slope (where the slope of the point is significantly greater than zero). We considered as a minimum TPI value 0.3 and therefore, not all cells under this value were considered as C builds-up. This approach facilitates the extraction of C build-ups from surrounding seafloor and reduces the occurrences of DEM artifacts. TPI scale (1–5) and value (0.3) were selected through a trial-and-error method in order to maintain the high-resolution of the extraction on which it depends the volume computation. Resulting rasters were re-classified and converted to obtain polygonal vector format. The remaining artifacts were manually deleted.

**Figure 3.** Topographic Position Index (TPI) maps, the color scale applied to the map indicate TPI value (low value in blue and high value in red) for the entire bathymetric dataset: (**A**–**E**) refer to the boxes in Figure 1, (**F**) histogram of the distribution for TPI value.

For each C polygon the shape index (SI) [33–35] was calculated that measures the complexity of patch shape compared to a circle shape. [11] demonstrates how SI is a useful tool to describe a seafloor landscape characterized by distinct C morphotypes.

The approach here proposed in performing the geomorphometric analysis was designed to specifically detect C build-ups and determine their 3D extension over the surrounding seafloor.

#### *2.3. Acoustics Profiles Analysis*

SBP dataset was analyzed using PETREL (Schlumberger) software in order to map the thickness of the C build-ups (Figure 4) along all survey lines and to explore the seismostratigraphic relationships between C build-ups and the substrate over which C has settled down. The level of details offered by our SBP data allowed us to assume that the depth of the substrate over which C growth at the time of its settlement is likely the same as the present-day seafloor that surrounds C build-ups. We, therefore, considered the seabed surrounding C build-ups as the base surface from which we extracted C build-ups as described in detail within the next section (Section 2.4).

**Figure 4.** (**A**) Thickness map of the C build-ups calculated using SBP interpretation. (**B**–**D**) profiles represent examples of C echotype in the area (black arrows).

#### *2.4. Volume Estimate*

The volume calculation was developed by the creation of a reference surface without C build-ups that was created for each DEM through ESRI ArcGIS geostatistical analysis toolbox. The interpolation function, used for the creation of the reference surface, was the natural neighbor [36]. A natural neighbor interpolation algorithm is a weighted average interpolation technique not affected by anisotropy or variation in the data density because the selection of the neighbors is based on the configuration of the data [37]. It does not infer trends and will not produce peaks, pits, ridges or valleys that are not already represented by the input samples [37]. The surface passes through the input samples and is smooth everywhere except at the locations of the input samples. It adapts locally to the structure of the input data, requiring no input from the user pertaining to search radius, sample count, or shape. We decided to use the natural neighbors interpolator with the aim of leaving a coarser morphology, thus avoiding smoothing effects given by other methodologies such as Spline or Kriging [38]. The volume computational was performed by the comparison of the analyzed DEMs with the corresponding reference surface in order to calculate the volume of all detected C build-ups. ArcGIS provides the cut/fill tool that summarizes the areas and volumes of change from a cut-and-fill operation, if the input raster surfaces are coincident and with the same cell size. The attribute table of the output raster presents all changes in surface volumes generated by the cut/fill computation. Positive values for the volume difference indicate regions of the original raster surface that have been cut (material removed). Negative values indicate areas that have been filled (material added).

#### **3. Results**

#### *3.1. Coralligenous Extraction Accuracy*

The model extracted 172.548 polygons but only 103,494 positive morphologies (Figure 5) were finally related to C build-ups by visual inspection of the hillshade bathymetry and validation from video observations collected within the study area. Most of these artifacts have been attributed to: (a) bad roll correction that creates false positive elongated structures (Figure 5b,d), (b) occurrence of *Posidonia oceanica* dead matter on the shallow area that results in reliefs with typical sub-circular shape [39,40] (Figure 5e), and (c) artifacts concentration on DEM discontinues on its boundaries.

**Figure 5.** Result from geomorphometric analysis. (**A**–**F**) refer to the boxes in Figure 1. Red polygons represent coralligenous build-ups extracted, polygons in yellow are examples of artifacts filtered from the study.

#### *3.2. Acoustics Profile Characterization*

The 38 SBP profiles analyzed in our study are mostly characterized by an indistinct echotype [41], where few discontinuous reflectors can be detected, in agreement with that previously observed by [11]. When associated with the occurrence of C build-ups, the first return from the seafloor is imaged by an indistinct, low-amplitude, and highly transparent reflector, whereas lateral continuity and overall amplitude increase crossing the surrounding seafloor. A similar seismic structure for algal reefs was described in [22], where the underlying acoustic basement was also visible. Interestingly, in our surveyed area, C build-ups are in contact with an erosional truncation that marks the top of the underneath sequence (Figure 6). The erosional truncation can be associated to aerial erosion of past highstand marine deposits, or to the ravinement surface formed at the onset of the post-glacial sea-level rise, that on the South Adriatic shelf margin enhanced transgressive erosion of older deposits [42]. Both the reflectors (the seafloor and the erosional truncation) were digitalized and exported to calculate the thickness of C frameworks for all survey lines (Figure 6). In the entire area, C shows a width between 0.1 m to 2.6 m with an average of 1.04 m (Figure 6).

**Figure 6.** (**A**) SBP dataset with a selected profile (in green) crossing red polygons representing C build-ups obtained from the geomorphometric analysis. (**B**) Selected SBP profile showing in detail the echotype related to C occurrences and the erosional truncation that marks the top of the underneath sequence (blue horizon). Blue and red X represent intersect profiles.

#### *3.3. Area and Volume*

TPI parameters extracted the distribution of the C build-ups with high resolution in terms of perimeter boundary. Volume calculation is heavily dependent on the reference surface, consequently dependent on the interpolation algorithm (a large data gap could not be filled in with a good resolution by interpolation). Volume calculation was performed following two different approaches, according to the different C morphotypes described by [11] and mapped in our survey area. In order to not force the interpolation and therefore the volume calculation, especially for tabular bank morphotypes, SI of C morphotypes was used as a threshold in selecting the most appropriate approach. In particular: (Table 1, Figure 7):


**Table 1.** Classification of C polygons modified from [11] based on SI and results in terms of area and volume.


The *cut and fill* tool was also used to analyze the performance of the segmentation process and the interpolation algorithm. The segmentation process directly affects the extraction of C build-ups (expected to be positive) as it will cut the DEM and consequently how the surrounding cells will influence the interpolation algorithm for the reference surface creation. In this way, with the extraction of positive structures from the DEM, positive balance was expected after the cut and fill tool for a correct volume computation. In our analysis, only in the deepest DEM, 0.3% presents a negative volume balance. This could be attributed to a small group of cells on the border of DEM that has been influenced by the surrounding lack of data.

**Figure 7.** Distribution of C features extracted with TPI colored by morphotypes, following the classification proposed in [11]. (**A**–**E**) refer to the boxes in Figure 1, they are details from the entire dataset and show: Red polygons for discrete relief, yellow polygons for hybrid bank and blue polygons for tabular bank. (**F**) Scatterplot indicates the direct correlation between area and SI, Spearman's ρ = 0.96.

#### **4. Discussion**

#### *4.1. TPI-Based Feature Extraction*

The results of this work illustrate the importance of combining high-resolution acoustic techniques and geomorphometric analysis in order to have a preliminary quantitative characterization of C habitat distribution. To our knowledge, this is the first attempt to quantify the volume of Mediterranean algal reefs. Conventionally, segmentation of MBES data sets into defined submarine landforms is carried out manually despite the process might be highly subjective, slow and potentially inaccurate [43,44]. On the contrary, geomorphometric techniques (e.g., terrain attributes, feature extraction, automated classification) can objectively characterize seabed terrain from the coastal zone to the deep sea [45–52], although a common and standardized technique for automated seabed classification has never been developed [53]. Only recently other works successfully develop object-oriented methods applying object-based image analysis (OBIA) or considering a complete set of remote data in order to precisely characterize targeted landforms on the seabed that will document the extension of biodiversity hotspots [52,54–56]. In this study, we considered a C build-up as a discrete feature in both two and three-dimensional spaces in order to use a geomorphometric parameter for object-oriented modeling and extraction from the seafloor. The method developed here (Figure 8) intentionally chose

a geomorphometric algorithm commonly available in existing terrain analysis software to achieve maximum applicability.

**Figure 8.** Conceptual model of the workflow for the extraction of volume and area of C features.

Variability of C morphotypes [11] poses several challenges to their automated extraction from DEM. The 3D complexity reached by C tabular banks [11] hampers the automatic detection of their edges, and consequently the creation of the hypothetical base surface for volume calculation. Since C build-ups raise from the surrounding seafloor (up to 4 m in height), their detection could be performed by slope analysis that is in general efficiently used to operate an accurate segmentation of isolated, small-scale features [47,48]. Nevertheless, the slope does not efficiently work for the tabular banks that usually have a lateral extension of more than 1 km2. In this case, although the slope is efficient in detecting the boundary of C (determined by a high difference in elevation) it is not able to include within the segmentation process the inner part of tabular banks, where the high tridimensional complexity makes difficult the creation of a continuous polygon. On the other hand, geomorphometric parameters that have a higher performance in detecting the 3D complexity, such as the rugosity index (i.e.,: vector ruggedness measure [57,58] or terrain ruggedness index [59]) are able to better define the entire distribution of C tabular banks, although they fail in giving an accurate estimation of the dimension of smaller morphotypes like discrete relief. TPI represents, therefore, a good compromise for C morphotypes detection, since it measures the relative topographic position of the central point as the difference between the elevation at this point and the mean elevation within a predetermined neighborhood.

The time-consuming operation in manual filtering the erroneous polygonal could be definitely avoided using a high quality larger multibeam coverage without artifacts. The accuracy of the detection of C build-ups can be appreciated by the results of the cut and fill tool explained in the previous section.

#### *4.2. Volume Computation*

Volume computation is obviously constrained by the substrate depth from which C build-ups started to grow. Our approach considered the seafloor surrounding all detected C build-ups, as the original substrate over which C build-ups started to develop. The analysis of SBP dataset and the associated resolution supports our assumption since seismic profiles clearly show the continuity of the erosional truncation, detected at the base of C build-ups, with the surrounding seafloor bare of build-ups, below less than 2 meters of mobile sediment (visible also on videos, Figure 2). We considered two different approaches in computing the volume of C build-ups, according to the two different groups of C morphotypes (Table 1). Our volume computation strategy took into consideration the influence of the 3D geometry of C build-ups on the performance of the computational process, in order to not force the interpolation algorithms. For smaller data gaps (associated with C build-ups with SI ≤ 2) we used the results coming from the cut and fill tool. The volume of C tabular banks (with SI > 2, and therefore polygons with larger areas), that represent big gaps to solve by interpolation, was calculated using the average of C thickness value obtained from SBP data analysis. The analysis of the acoustic profiles was critical in supporting our volume computation strategy for tabular banks. Despite the considerable relief that C can reach in the survey area (up to 4 m if calculated after TPI extraction), SBP data show a maximum thickness of 2.5–3 m (only a few cells up to 4 m, Figure 4) decreasing to 1–1.5 m in the inner part of the C tabular banks.

#### *4.3. Coralligenous Growth Rate and Age*

Since present-day C depth distribution roughly matches those areas of the seafloor that were emerged during the last glacial maximum, C build-ups likely started to occur during the last transgression [60]. The C growth rate has been estimated by radiocarbon dating with a mean growth rate of 0.19 mm yr<sup>−</sup>1, and a range of 0.11 to 0.26 mm yr–1 for C build-ups sampled between 10–35 m of water depth [61,62]. The results of this work suggest that a tridimensional quantitative characterization of C habitat could be related to the C growth rate. Using observed thickness values (Table 1) and assuming a prevailing vertical growth direction, we obtained respectively an age of at least 1.7 kyrs BP for discrete relief and hybrid banks, and 7.8 kyrs BP for tabular banks. Both these results are in agreement with the general estimated Mediterranean C build-ups age [4,60–62]. Finally, we calculate the C volume production *per* year that results in 22.04 m3 yr–1 for discrete relief and hybrid banks and 35.38 m<sup>3</sup> yr–1 for the tabular banks morphotypes.

#### **5. Conclusions**

We demonstrate that our ad-hoc semi-automated, GIS-based methodology is useful to portray the three-dimensional complexity of C build-ups on a sector of the southern Apulian continental shelf and to quantify the extent to which C builders shaped the submarine environment, affecting the evolution of submerged landforms.

Notably, the maps and volume computation obtained from our work represents the first quantitative data supporting an estimation of the contribution of C build-ups to the ocean carbon budget. Our dataset also yields new details about the genesis of C build-ups (i.e., their potential relationships with drowned continental shelf morphologies), by imaging the development of C tabular banks in the study area.

Finally, despite the lack of a comprehensive investigation on the role of quaternary geomorphic processes on the C distribution and growing processes, the obtained results will support any future exploration of the relationship between framework morphologies and changes observed in the associated biological community, which still represent a major gap in knowledge in eco-geomorphological research of Mediterranean algal reefs.

**Author Contributions:** F.M. collected the entire dataset, performed and designed the methods. F.M. wrote and coordinated the manuscript with the significant contribution of A.S. V.A.B. and D.B. supervised the Coralligenous sections. F.M. and G.L. processed the SBP data. C.C. supervised the BIOMap Project. A.S. supervised the research. All authors contributed substantially to the discussion of the results and revision of the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was financed by the BIOMAP project (P.O FESR 2007/2013, Puglia Region, Italy).

**Acknowledgments:** We are grateful to Michele Panza for the assistant during Calafuria Issel surveys (2013). SBP lines are property of Acquedotto Pugliese SPA and have been acquired in the framework of a commercial survey executed by the CoNISMa Local Research Unit of Milano-Bicocca. We thank Acquedotto Pugliese SPA for the scientific use of these data. We also thank three anonymous referees for their valuable comments on sections of this manuscript. FM and VAB are funded through a post-doctoral fellowship in Earth Sciences by the University of Milano-Bicocca. This article is also an outcome of Project MIUR—Dipartimenti di Eccellenza 2018–2022.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
