**Benthic Habitat Morphodynamics-Using Remote Sensing to Quantify Storm-Induced Changes in Nearshore Bathymetry and Surface Sediment Texture at Assateague National Seashore**

#### **Arthur Trembanis 1,\*, Alimjan Abla 2, Ken Haulsee <sup>1</sup> and Carter DuVal <sup>1</sup>**


Received: 25 August 2019; Accepted: 16 October 2019; Published: 18 October 2019

**Abstract:** This study utilizes repeated geoacoustic mapping to quantify the morphodynamic response of the nearshore to storm-induced changes. The aim of this study was to quantitatively map the nearshore zone of Assateague Island National Seashore (ASIS) to determine what changes in bottom geomorphology and benthic habitats are attributable to storm events including hurricane Sandy and the passage of hurricane Joaquin. Specifically, (1) the entire domain of the National Parks Service offshore area was mapped with side-scan sonar and multibeam bathymetry at a resolution comparable to that of the existing pre-storm survey, (2) a subset of the benthic stations were resampled that represented all sediment strata previously identified, and (3) newly obtained data were compared to that from the pre-storm survey to determined changes that could be attributed to specific storms such as Sandy and Joaquin. Capturing event specific dynamics requires rapid response surveys in close temporal association of the before and after period. The time-lapse between the pre-storm surveys for Sandy and our study meant that only a time and storm integrated signature for that storm could be obtained whereas with hurricane Joaquin we could identify impacts to the habitat type and geomorphology more directly related to that particular storm. This storm impacts study provides for the National Park Service direct documentation of storm-related changes in sediments and marine habitats on multiple scales: From large scale, side-scan sonar maps and interpretation of acoustic bottom types, to characterize as fully as possible habitats from 1 to 10 m up to many kilometer scales, as well as from point benthic samples within each sediment stratum and these results can help guide management of the island resources.

**Keywords:** side-scan sonar; swath bathymetry; habitat monitoring; hurricane Sandy; hurricane Joaquin; climate change; shoreline detection; remote sensing

#### **1. Introduction**

The overall goal of this study was to map the morphodynamic changes to the nearshore zone of Assateague Island National Seashore (ASIS) to determine what changes in bottom sediments and benthic habitat occurred are attributable to storm impacts such as Superstorm Sandy and hurricane Joaquin. Specifically, we (1) mapped the full survey area with side-scan sonar and multibeam bathymetry at a resolution comparable to that of the pre-storm survey, (2) resampled a subset of the benthic stations that represented all sediment strata previously identified, and (3) compared newly obtained data to that from the pre-storm survey to examine storm impacts. This survey design aimed

to document storm-related changes in sediments and marine habitats on multiple scales: From large scale, side-scan sonar maps and interpretation of acoustic bottom types, to characterize as fully as possible habitats from 1 to 10 m up to many kilometer scales, as well as from point benthic samples within each sediment stratum.

#### *1.1. Background on Acoustic Mapping*

Human reliance upon, and interference with benthic ecosystems necessitates an understanding of the spatial extent, structure, and function of these unique ecosystems [1,2]. This project work utilized advanced geoacoustic techniques for remote benthic habitat mapping and employed both traditional technologies for wide area coverage combined with point sampling and robotic systems including an autonomous underwater vehicle (AUV) and a remotely operated vehicle (ROV) for gathering data in an unprecedented level of resolution over targeted regions of interest. Knowledge gained from this study thus allows coastal and estuary environmental stakeholders responsible for the Assateague Island National Seashore (ASIS) to better manage resources to protect these environments by providing information on the location of habitats, species' distributions and associations to varying sedimentary and morphological regimes. The new benthic habitat and geomorphologic datasets will serve as a useful reference for an array of mapping and management efforts particularly providing insights into temporal dynamics from before and after storms.

Estuarine and coastal regions of the US face multiple marine spatial planning issues and challenges, including the effects of episodic storms, climate change, and other stressors [3]. In particular, there has been a strong interest in and multiple efforts to examine the effects of hurricane Sandy on the seabed within the Mid-Atlantic Bight [4–6]. The anticipated impacts of offshore development on the seabed and to the ecology and, more generally, initiatives in marine spatial planning [3] point to an ever-increasing need for both base-line mapping and monitoring efforts directed at biological habitats and geological features of the littoral zone.

In this study, we realized the need to fuse disparate data streams and recognized the most effective means of analyzing and understanding them is through data visualization. By repeated benthic habitat and morphology mapping, this field effort and the subsequent data analysis provides a critical extension to the previous study of the ASIS site (2011–2012). We recognize that together these data sets collected over multiple years represent a knowledge base much greater than the sum of its parts, and it is this fact that motivated the project.

Recent benthic mapping on similar prior projects to this project [7–9] have demonstrated that, in general, shallow coastal ecosystems are far more complex and heterogeneous in bottom type than typically anticipated. Furthermore, conventional sampling by direct methods (i.e., bottom grabs or dredges) is also generally recognized as severely limited by several factors including gear bias, spatial averaging and limited sample density. It is only from the careful fusion of multiple technologies (i.e., ship-based, ROV, AUV, grab samples, and dredge) that one is able to compile the complete picture of the nature and composition of the bottom [8–11].

Accurate bathymetric and backscatter data are essential for hydrographic charts, dredging, and habitat mapping. Previously, the primary method for seafloor mapping was the use of single-beam bathymetry. It is an accurate technique for collecting seafloor data, yet there are a number of limitations. In single-beam surveys, a discrete footprint of the seabed is ensonified leaving large gaps between adjacent survey tracks that are not surveyed. Single-beam surveys typically cover only 5-10% of the total seafloor area leaving large portions unsampled. Therefore, the maps generated from single-beam sonar samples may lack high-enough resolution to detect small and ecologically important habitat features [12]. In a recent study [12], the authors found that a full coverage swath sonar system resulted in a 90% classification accuracy compared to 70–80% accuracy from a single beam system approach. Using conventional methods like RoxAnn (e.g., [13,14]), found that the sampling is limited to a single sonar point footprint directly below the vessel, and maps necessitate heavy interpolation between track lines. In shallow settings, single-beam systems therefore require progressively more and more tightly

spaced lines as the local depth decreases. In contrast, multi-beam echosounders, interferometric sonar, and side-scan sonar systems characterize a wide swath on either side of the survey vessel. Surfaced towed acoustic sensors suffer from positioning errors associated with the ambiguity of the towfish location relative to the ship global positioning system (GPS). Surface vessel mounted sensors also suffer from reduced resolution (larger footprint) as water depth increases. Wide swath acoustic methods afford the potential of effectively mapping large areas of the seafloor at a resolution appropriate for habitat mapping and correlation with biological communities. Seafloor echoes can be interpreted in terms of hardness and roughness and correlated with the sediment and faunal characteristics.

Several reviews and studies of coastal sedimentary systems (e.g., [15–17]) have repeatedly demonstrated that dramatic spatiotemporal heterogeneity dominates seabed settings with respect to composition and morphology and that this heterogeneity is more the norm than the exception. Further exacerbating the problem, seabed heterogeneity often exists at spatial scales on the order of 1–10's m, well below typical single-beam survey spacing and gridding schemes. Understanding and being able to map this spatial patchiness at sufficient coverage and resolution scales is a critical component of habitat mapping efforts and an important motivation behind this project.

#### *1.2. Study Design and Objectives*

This study was designed to characterize the physical and biological components of the nearshore benthic zone of Assateague Island National Seashore (ASIS), including delineations of habitat types using the CMECS system. We characterized benthic sedimentary habitats by analyzing Young grab samples for grain size distribution, and total organic carbon (TOC). The study methods were modeled off of the 2011 ASIS benthic habitat assessment projects undertaken by Versar, Inc. (Columbia, MD) and the Maryland Geological Survey (MGS), originally contracted by the National Park Service (NPS). A Natural Resource Condition Assessment undertaken by ASIS and the University of Maryland [18] ranked the Atlantic subtidal habitats as reaching 99% of their reference condition (as determined by historical observations and data), albeit with very limited confidence or knowledge of indicators for the reference condition. Thus, data from the assessments by Versar, Inc. and the Maryland Geological Survey serve as a crucial pre-Sandy (2011) point of comparison for the possible effects of the physical stressors of hurricane Sandy on sediment distribution. While this project was one of four simultaneous benthic mapping regimes funded by the NPS in the aftermath of Sandy, this study covered the largest habitat area, focused on the Atlantic subtidal zone, and was unique in its availability of pre-Sandy (2011) reference data.

In this paper we set out to: (1) Provide an updated subtidal resource inventory to the National Park Service by determining benthic distribution utilizing geoacoustic surveys compelled with grab samples for grain size analysis; (2) analyze pre-Sandy (2011) and post-Sandy (2014-2015) habitats with a uniform statistical approach and compared them to determine what changes occurred through the study period; (3) conduct targeted acoustic mapping of areas of particular interest before and after storms to provide an unprecedented higher resolution of these target locations.

#### *1.3. Study Site*

Assateague Island is a 58-km long barrier island stretching from the Ocean City inlet in Maryland, down past Chincoteague Island in northern Virginia. This dynamic island is characterized by barrier beach and back-bay environments that are subject to the constant eroding and re-shaping forces of the Atlantic Ocean. The pre-storm and post-storm surveys of Assateague Island subtidal resources discussed here include the majority of the length of the island (Figure 1). We sampled benthic sediments up to 1 km offshore in accordance with the pre-storm survey and NPS jurisdictions (Figure 2), whereas the USGS surveyed areas further offshore of Assateague Island. Historically the nearshore subtidal habitats at ASIS are comprised mostly of fine sand, with patches and swaths of mud and coarser sediments. Much of the observed heterogeneity in sediments is due to the relict sediments left behind by the characteristic landward 'rolling-over' of barrier islands.

**Figure 1.** Pre-Sandy (2011) survey sites (n = 125) sampled in October 2011 on the beach (n = 15) and subtidally up to 1 km offshore (n = 110) of Assateague Island National Seashore, Maryland.

**Figure 2.** Post-storm (n = 120 grabs) sampled in October 2014 (n = 60 grabs, white) and 2015 (n = 60 grabs, yellow), up to 1 km offshore of Assateague Island National Seashore survey sites, Maryland.

#### **2. Materials and Methods**

#### *2.1. Acoustic Mapping: Side-Scan Sonar Workflow*

The side-scan sonar data collected was processed using Chesapeake Technology, Inc., Mountain View, CA, USA SonarWiz 6.0 sonar processing software. Sonar data processing followed industry standard procedures for side-scan data, focusing on gain corrections to remove the variation in across track brightness inherent in side-scan sonar performance in order to achieve the most representative mosaic of the seafloor. Standardized gain settings were chosen to make the sonar data as internally consistent as possible between different daily mission datasets, while also attempting to replicate the products produced from the Versar sonar data collected off Assateague in 2011. Raw Edgetech sonar

files (.jsf) were imported using a time-varying gain, to compensate for signal loss on the outer swath of the sonar files. Once the files were imported, bottom-tracking corrections were applied if automatic bottom tracking had failed. Once verified or corrected, an empirical gain normalization correction was applied. This gain correction makes the data set consistent from a file-to-file basis, removing gain biases caused by variations in backscatter intensity from differing sediments between each file. This produces an interiorly consistent mosaic. From this work flow (Figure 3), several mapping products were created. A standard georeferenced image file (geotiff) and Google Earth georeferenced file (.kmz) are generated. Vessel navigation or track-lines files were exported as shapefiles (.shp). Each gain-corrected file was exported as a waterfall image (.jpg). Target reports were generated (in .pdf format), highlighting unique features or objects on the seafloor. Finally, coverage reports (.xls) list the metadata for each file (start/end latitude and longitude, length, time, etc.), as well as the total area and line-line overlap achieved for the complete mission.

**Figure 3.** Side-scan sonar post processing workflow used by the University of Delaware, Assateague Island National Seashore mapping project: From EdgeTech raw data to SonarWiz to shapefiles and Google Earth .kml data products.

#### *2.2. Bathy Workflow*

The bathymetric data collected was processed using Chesapeake Technology, Inc., Mountain View, CA, USA SonarWiz 6.0 sonar processing software. Sonar data processing followed industry standard procedures for bathymetric data, focusing on injecting tidal offsets, sound velocity offsets, and configuring the vessel files in order to achieve the most representative mosaic of the seafloor. The bathymetric processing workflow (Figure 4) was standardized to keep the bathymetric representations as internally consistent as possible. The first step was to configure a vessel offset file for the vessel that

was recording the sonar data. All of the research vessels measurements were captured including the inertial measurement unit (I.M.U.) to the echosounder and the I.M.U to the antennas were recorded and saved in a Sonarwiz6 vessel file. After the vessel file is created, the raw sonar data (.jsf) was imported into a new Sonarwiz6 project that was set-up with an internally consistent naming convention based off the date of the survey. Before the data was merged, the sound velocity profiles, tide file, and the patch test values were imported into the project. After the echosounder was installed and before surveying operations began, a patch test was completed to ensure that each installation of the echosounder provided internally consistent data. The process accounts for offsets in roll, pitch, time latency, and heading. These patch test values were then entered into the vessel offsets file. In the field, sound velocity profiles were collected every hour on a Castaway CTD. The sound velocity files (.csv) were then imported into the Sonarwiz6 project. A tide file from the days that we surveyed was downloaded from the National Oceanic and Atmospheric Administration (NOAA) Tides and Current website in the form of .csv and then imported into the Sonarwiz6 project. After the patch test values, sound velocity profiles, and the tide file were imported into the project, the data was merged. This merging process injects all of the above information onto the XYZ point cloud that was recorded during the survey to constrain the positional error, while providing the most replicable data products that are internally consistent. Gridded mosaics, at a resolution of 0.75 m, were then exported as a .grd file using the CUBE algorithm, which is now the industry standard for providing the most accurate gridded solutions. The backscatter from the bathy was also processed in Sonarwiz6 and was exported as a .grd file. In all, a geotiff, .kmz, .xyz, and .grd of both the bathymetry data and the backscatter data were exported as bathymetric data products.

**Figure 4.** Schematic of workflow used by the University of Delaware mapping team for multibeam bathymetry data from EdgeTech 6205 acquisition to gridded data products. Data collected as part of a submerged mapping study for Assateague Island National Seashore in 2014–2015.

Based on the accuracy and resolution of our navigation system (a Coda F190R with dual head RTK and OmniSTAR enabled GNSS receivers) and standard patch test performed during our survey along with confirmation from overlapping passes we estimate that the horizontal and vertical uncertainty is approximately 25–30 cm respectively. Previous studies conducted using this same positioning system, sonar, and vessel in a further offshore site with the exact same sensors and configurations yielded repeatable target acquisition of within 1 m horizontally and was reported by [19]. Another study by the U.S. Geologic Survey (USGS) utilizing the same inertial measurement unit (IMU) and satellite correction service also yielded similar uncertainty values of general less than 50 cm for 95% of the

hydrographic soundings [20]. Therefore, it is reasonable to set our achievable positioning reliability and overall total propagated error at approximately 50 cm horizontally and vertically. Our baseline inertial measurement accuracy was approximately 5cm horizontal and vertical and was aided by either OmniSTAR corrections or RTK base stations when available. Even without the satellite or RTK corrections our positioning would have been accuracy to a meter or two and thus much less than the scale of the changes seen across the seabed that range from 20 to 180 m of variability.

#### *2.3. Acoustic Seabed Classification Map Preparation*

#### 2.3.1. Repeated Seabed Mapping Surveys

Our three, 2014–2015 acoustic surveys allowed us to make both inter and intra-annual comparisons for storm-related changes in bottom sediments. The first survey resulted in side-scan sonar mosaics and sediment acoustic class shapefiles from June 2014. These data include 5 side-scan sonar mosaics and sediment class shapefiles that were generated from the side-scan sonar data collected during 20–25 June 2014. The second survey was the most broadly defined and resulted in side-scan sonar mosaics and sediment class shapefiles for May 2015. These data include 12 side-scan sonar mosaics and sediment class shapefiles that were generated from the side-scan sonar data collected during 12–21 May 2015.

Finally, we conducted a set of spot selected surveys in October 2015 producing side-scan sonar mosaics and sediment class shapefiles. This survey was spatially focused on areas of change following hurricane Joaquin. Data were collected during 13–16 October 2015.

#### 2.3.2. Acoustic Seabed Classification

First, all the available datasets were imported into ArcGIS platform under coordinate system - WGS\_1984\_UTM\_Zone\_18N (EPSG:32618). Then, manually digitized acoustic classes were matched against side-scan sonar mosaics and grain size data of grab samples to make sure the side-scan sonar mosaics were segmented correctly. After that, sediment class boundaries were matched between neighboring sediment classes to correct any inconsistencies in the interpretation of side-scan sonar mosaics, (i.e. to keep the class assignments consistent from one survey day to the next). Minor discrepancies with the side-scan sonar data interpretation and inconsistencies in sediment class boundaries were corrected manually using ArcGIS Edit tool.

After fixing all the discrepancies related with side-scan sonar data interpretation, all the same sediment classes in different shapefiles were merged into a single shapefile using the "Merge" tool in Arc Toolbox (Figure 5). Overlapping areas of the same sediment classes between different shapefiles were dissolved into one using the "Dissolve" tool in Arc Toolbox. At the end, 4 single shapefiles representing 4 different sediment classes namely Coarse Sand, Medium Sand, Fine Sand, and Mud, were merged into one single shapefile named Sediment Class\_201406\_201505.shp. A new field named "Area" was added to the table of contents of the newly generated shapefile and areas in m<sup>2</sup> were calculated for the sediment classes using "Calculate Geometry" method in ArcGIS.

The acoustic backscatter data collected in the field surveys was used to generate mosaic sonar images as outlined in the section above and then digitized using ArcGIS to manually outline and build categorical maps of the Assateague Island survey domain. The same acoustic class types as used in the 2011 precursor study by Versar/MGS were adopted here for consistency and the previous acoustic class maps were used along with the newly acquired and created side-scan sonar mosaics to guide the human operator in the segmentation process. All segmentation effort was conducted by the same team member for consistency and then reviewed by the principle investigator. The backscatter segments were then assigned to the acoustic class and were assigned similar colors again adopting the same nomenclature and color scheme as used in the precursor study to aid in comparison. Note that up to this point the classes are simply acoustic classes and require some external ground-truthing to further relate the classes to actual seabed types. The next step of relating the acoustic classes to seabed types comes through the correlation of class maps to the benthic grab samples and grain size analysis samples.

Figures 3 and 4 illustrate the workflows used in data processing both the side-scan sonar and phase measuring bathymetric sonar.

**Figure 5.** Flow-chart of data processing steps for manually segmented side-scan sonar-sediment classification into a habitat map for Assateague Island National Seashore, developed by the University of Delaware 2014–2015.

#### 2.3.3. Ground Truthing

Essential for any geoacoustic marine habitat mapping initiative is the need to gather ground truthing observations to compare to the acoustically derived maps. Throughout our previous habitat mapping campaigns we have utilized benthic grab samples for ground truthing [6,9]. Grab samples were analyzed for grain size distribution and benthic faunal composition using standard approaches (see below and references therein for details). We used a modified Young grab sampler for recovering surface sediment ground truthing samples. The available grain size data included grain size analysis results from 98 grab samples taken during June of 2011 and 146 grab samples taken during October of 2014 and October of 2015 for benthic samples plus additional ones for ground truthing only. However, only the ones taken during 2014 and 2015 were used to check the side-scan sonar data interpretation in this study. Additionally, a small portable YSI Castaway CTD (http://www.ysi.com/castaway) was used to gather vertical profiles of salinity, temperature, and sound speed for use in environmental characterization and swath bathymetry data processing.

#### **3. Results**

#### *3.1. Acoustic Mapping Products and Findings?*

In total, over 65 km<sup>2</sup> of geoacoustic mapping was accomplished along the entire length of Assateague Island from as near to the shore as was possible out to at least 1km offshore. Additional sonar lines were acquired in order to connect the nearshore ASIS park boundary to the regional USGS mapping effort (Figure 6). The mapping area covered depths from 2 m–12 m extending from the beach to ~2 nm offshore, with final map resolution 50 cm/pixel for side-scan sonar and 1 m2/pixel for bathymetry products. The field survey design as developed and agreed to by all project partners conducting mapping at this and other post-Sandy NPS sites was to achieve 100% side-scan sonar coverage and as much bathymetric coverage as possible given the constraints imposed by such shallow water survey operations. In this study we were able to achieve approximately 50% swath bathy coverage (Figure 7) and 100% side-scan sonar coverage (Figure 8).

**Figure 6.** Overview of study site area showing park boundary (green) and recent USGS coverage along with UD mapping team gap filling coverage (red).

**Figure 7.** Bathymetry mosaic (1m/pixel) from 2014–2015 acoustic surveys (550 kHz) for the entire study area. Notable features captured in the bathymetry included pronounced shore-attached oblique bars.

**Figure 8.** Side-scan mosaic from 2014–2015 acoustic surveys for the entire study area.

Resulting swath bathymetry data for the study side was gridded to a resolution of 1m/pixel. Reliable bathymetry data was obtained for approximately half the total sonar swath width resulting in 50% seafloor coverage in depths ranging from 2 to 12 m. Early sampling in 2014 around Fishing Point was conducted using spacing resulting in 100% bathymetry coverage (200% side-scan sonar) but the survey rate was unsustainable for the entire study site and it was agreed upon by all project partners to expand the line spacing to achieve 100% side-scan and 50% bathymetry. This survey strategy is also in keeping with the same approach used by the USGS in their larger regional coverage provided for a meaningful connection of this study to that regional effort. The resulting swath bathymetric grids extend what was only single-beam echosounder data from the precursor 2011 project and thus provide an enhanced baseline for AINS resource managers. Notable features captured in the bathymetry included pronounced shore-attached oblique bars.

A complete set of side-scan sonar mosaics with a resolution of 50 cm/pixel was generated for the entire length of the island and included extensions to fill gaps between the ASIS park boundary and the USGS regional survey effort. Gain settings and post-processing were set to generate an internally consistent sonar backscatter map (Figure 8) that allowed for consistent segmentation and interpretation with respect to sediment class. The convention used in the backscatter mosaics is for bright yellow color for high-amplitude returns (i.e. coarse sand and gravel) and dark brown/black for low amplitude returns (i.e. mud). Bright high-amplitude returns in the backscatter mosaic were clearly associated with the raised features of shore-attached oblique sand ridges noted in the bathymetry (Figure 7). Mud patches were most notable in the adjoining troughs associated with the sand ridges or from relict marsh outcrops exposed by scour and barrier island migration.

As is shown in Figure 9, the sediment classification map includes four different sediment classes were recognized and shaded with different colors.

**Figure 9.** Sediment classification from 2014–2015 for the entire study area. Four sediment classes (coarse and medium sand, mud, and fine sand) were recognized and shaded with different colors. Mapping survey conducted 2014–2015.

Segmentation and acoustic classification of the side-scan sonar backscatter mosaic was performed via manual heads-up segmentation utilizing the same class types from the 2011 precursor study but informed by the sediment grab samples collected in this study. The GIS shapefile vectors for each acoustic class type are shown in Figure 9 using the same color and naming convention from 2011 for ease of comparison. The distribution of class type (Table 1) followed similar overall trends as documented in 2011 with medium and coarse sands associated with positive relief features like the large shore-attached oblique ridges and mud associated with the adjoining troughs. Notably absent in the 2014–2015 data were any mappable patches of gravel. The percentage area distribution of each class type are summarized in Table 1.


**Table 1.** Surficial sediment distribution by class type based on the 2014–2015 habitat mapping expedition.

The following three figures illustrate an example of an area with pronounced spatiotemporal changes in acoustic class texture between the 2011 precursor map (Figure 10) and repeated surveys with the EdgeTech 6205 in 2014 (Figure 11) and 2015 (Figure 12). The same backscatter amplitude convention is used but the gain settings and sonar frequencies are slightly different between the 2011 and 2014/15 surveys. Nevertheless, the general trend shows a migration of a set of coarse sand patches and the burial of a previously exposed mud patch.

**Figure 10.** Side-scan mosaic and classification shapefiles based on 2011 survey data.

**Figure 11.** May 2015 side-scan mosaic and sediment classification shapefile compared to 2011 from Figure 9.

**Figure 12.** October 2015 side-scan mosaic and sediment classification shapefile compared to May 2015 (Figure 11) and 2011 (Figure 10).

#### *3.2. Temporal Changes in Seafloor Sediment Type*

Select portions of the survey domain were mapped both in May 2015 and October 2015 capturing before and after the passage of hurricane Joaquin. Sediment class maps were determined in the same manner and using the same grain size informed class types in order to examine changes in the before and after storm surveys. Sediment types were coded from 1 to 4 from smallest to largest grain size with 1 as mud, 2 as fine sand, 3 as medium sand and 4 as coarse sand and these codes were manually inserted into the table of contents of sediment class shapefiles. Sediment class shapefiles were then converted into raster files with sediment class type code being the value field. In order to understand if the seafloor sediments are coarsening or fining from May 2015 to October 2015, a sediment property change map was generated by subtracting sediment class raster data of May 2015 from sediment class raster data of October 2015 (Figure 13). This map shows the pattern and intensity of change in sediment properties from May 2015 to October 2015 as a magnitude of sediment property change index (Figure 13). For example, positive sediment property change index means the seabed became coarser, negative sediment property change index means the seabed became finer and zero value means sediment property was unchanged. Higher absolute value means higher intensity of change, for example +3 area represents a change from mud to coarse sand −3 area represents a change from coarse sand to mud.

**Figure 13.** Sediment property change map for submerged resources of Assateague Island National Seashore between May 2015 and October 2015. Positive sediment property change index means seafloor became coarser, negative sediment property change index means seafloor became finer and zero value means sediment texture did not change.

Comparison of sediment class boundaries between May 2015 and October 2015 suggest a major shift (mostly south to eastward) in different sediment class boundaries. The distance of the shift in the fine-coarse sediment class boundary is about 20–40 m in the northern portion of the study area (Figures 14 and 15) and up to 80 m in the southern portion (Figures 16–19). The shift in mud-coarse sediment class boundary however, was measured to be up to 180 m in the southern part (Figure 17).

**Figure 14.** Side-scan sonar mosaic interpretation, Assateague Island National Seashore, Maryland (**A**) side-scan sonar mosaics collected during 13 May 2015; (**B**) interpretation of side-scan sonar data collected in 13 May 2015; (**C**) side-scan sonar mosaics collected in October 2015; (**D**) interpretation of side-scan sonar mosaics collected in October 2015. Note that red polygons represent the same areas mapped during both 13 May 2015 and October 2015.

**Figure 15.** Temporal changes in sediment class boundaries, Assateague Island National Seashore, Maryland. The underlying image on the bottom is the side-scan sonar mosaic collected during May 2015. Shaded polygons on the top represent sediment classes inferred from the side-scan sonar mosaics collected during October 2015. Note the southeastward shift (18–50 m) in the boundary of fine-coarse sand sediment classes.

**Figure 16.** Side-scan sonar mosaic interpretation, Assateague Island National Seashore, Maryland. (**A**) Side-scan sonar mosaics collected during 13 May 2015; (**B**) Interpretation of side-scan sonar data collected in 13 May 2015; (**C**) side-scan sonar mosaics collected in October 2015; (**D**) Interpretation of side-scan sonar mosaics collected in October 2015. Note that red polygons represent the same areas mapped during both 13 May 2015 and October 2015.

**Figure 17.** Temporal changes in sediment class boundaries, Assateague Island National Seashore, Maryland. The profile in the bottom is the side-scan sonar mosaics collected during May 2015. Shaded polygons on the top represent sediment classes inferred from the side-scan sonar mosaics collected during October 2015. Please note the northward shift (180 m) in the boundary of mud-coarse sand sediment class and Southeastward shift (37–40m) in the boundary of fine–coarse sand sediment class.

**Figure 18.** Side-scan sonar mosaic interpretation, Assateague Island National Seashore, Maryland. (**A**) Side-scan sonar mosaics collected during 14 May 2015; (**B**) interpretation of side-scan sonar data collected in 14 May 2015; (**C**) side-scan sonar mosaics collected in October 2015; (**D**) interpretation of side-scan sonar mosaics collected in October 2015. Note that red polygons represent the same areas mapped during both 14 May 2015 and October 2015.

**Figure 19.** Temporal changes in sediment class boundaries, Assateague Island National Seashore, Maryland. The profile in the bottom is the side-scan sonar mosaics collected during May 2015. Shaded polygons on the top represent sediment classes inferred from the side-scan sonar mosaics collected during October 2015. Please note the west-south-eastward shift (38–57 m) in the boundaries of fine–coarse sand sediment classes.

#### *3.3. Temporal Changes in Bathymetry*

Storms and extreme weathering events [21,22] or human induced disturbances such as dredging [23] have the potential to change the nearshore bathymetry. These changes can have serious impacts to the marine flora and fauna that most marine environment protection efforts are focused on.

The purpose of the bathymetric change comparison is to see if it is possible to determine significant changes in the bathymetry after the passage of hurricane Joaquin, which passed through the area in late September and early October of 2015. A total of 18 bathymetric cross sections were constructed using Interpolate Line tool in ArcScene across areas with overlapping survey coverage from before (May 2015) and after (October 2015) the passage of hurricane Joaquin. Nine of the cross sections (A1, B1 to J1) were made from May 2015 bathymetry, and the other nine cross section (A2, B2 to J2) were constructed from October 2015 bathymetry data along the same trackline as their pairs in May 2015. These bathymetry profiles were then exported as Excel spreadsheets. Sediment class type was added to each of the points along the bathymetry profiles based on the attributes from the acoustic sediment class maps discussed previously. Next, the Excel spreadsheets of the bathymetry profiles were imported into MATLAB R2017 (Natick, USA) to generate comparison charts between pairs of bathymetry profiles. Sediment types were coded from 1 to 4: 1 as mud, 2 as fine sand, 3 as medium sand and 4 as coarse sand. Bathymetry and sediment class comparison charts were generated by subtracting sediment type values of May 2015 profiles from October 2015 profiles, in order to understand if the seafloor sediments are coarsening or fining with the lateral shift in bathymetry (Figures 20 and 21).

In the areas where cross sections were constructed more or less parallel to the direction of sediment class boundary shift, southwestward lateral shift of sand bars (up to 60 m) were recognized in bathymetry cross section profiles (Figures 22 and 23). In the areas where cross sections were made more or less normal to the direction of sediment class boundary shift, vertical shifts of up to 65 cm are recognized (Figure 24). The largest vertical change was recognized in October 2015 profile J2 (Figure 24), which suggests axial incision of about 65 cm. This is the closest profile to the shoreline. Other profiles parallel to profile J2 suggests that vertical change in bathymetry seems to decrease with increasing distance from the shoreline.

**Figure 20.** Bathymetry profile change versus sediment class boundary shift after Hurricane Joaquin. Curved bathymetry profiles A1 (May 2015) and A2 (October 2015) are color coded according to the sediment types in the chart located in the top. Horizontal scale are meter in both panels, and vertical are meters in the top panel and size class index in the bottom panel. Sediment classes were indexed to 4, 1 representing the finest and 4 representing the coarsest. The chart in the bottom was generated by subtracting sediment type values of A1 (May 2015) from A2 (October 2015) to understand if the seafloor sediments are coarsening or fining with the lateral shift in bathymetry. Please note the coarsening of sediments in the top of the sand bars and fining of sediments in the trough. Inset scale is in meters on both axes.

**Figure 21.** Bathymetry profile change versus sediment class boundary shift after Hurricane Joaquin, Assateague Island National Seashore. Curved bathymetry profiles C1 (May 2015) and C2 (October 2015) in Figure 25 are color coded according to the sediment types in the chart located in the top. Horizontal scale are meters in both panels, and vertical are meters in the top panel and size class index in the bottom panel. Sediment classes were indexed from 1 to 4, 1 representing the finest and 4 representing the coarsest. The chart in the bottom was generated by subtracting sediment type values of C1 (May 2015) from C2 (October 2015) to understand if the seafloor sediments are coarsening or fining with the lateral shift in bathymetry. Please note the coarsening of sediments in the top of the sand bars and fining of sediments in the trough. Inset scale is in meters on both axes.

Apart from the axial incision recognized in October 2015 profile J2 (Figure 24), comparison of May and October profiles suggests that seabed changed from fairly smooth in May 2015 to incised by shore oblique/perpendicular bars in October 2015 (Figures 22 and 23). These are clear signs of hurricane effect on the seafloor. Decrease in vertical offset towards the sea may suggest an inverse relationship between hurricane effect on the seafloor and the distance from the shoreline. A Comparison between bathymetry change and sediment class boundary shift in this area suggests that sediments in the top of sand bars went coarser after the hurricane, while the ones in the trough of sand bars went finer (Figures 20 and 21).

#### **4. Discussion**

This nearshore benthic morphodynamic study at Assateague Island addresses a large gap in knowledge about the spatiotemporal dynamics of subtidal habitats and geomorphology in response to storm events (Figure 25). This knowledge will inform subsequent general management plans to best safeguard subtidal resources from anthropogenic and climate change stressors. This study also provides relevant data and recommendations for future assessments at ASIS.

**Figure 25.** Timeline of Storm events and seafloor mapping studies and resulting benthic habitat and morphology products.

Although storms have varied morphologic impacts on different benthic substrates (i.e. coarse sand versus mud), the ASIS nearshore zone experienced large lateral shifts on the order of 30–70 m of sediment type contacts resulting in changes from fine to coarse sediment type with additional vertical depth changes from erosion and migration of between 50 cm to as much as 2 m associated with migration of nearshore oblique sand bars as a result the passage of storms.

We are not able to directly link the geomorphology and habitat changes outlined here to Hurricane Sandy, excluding the possibility of impacts from other storms throughout the overall study period. This is largely due to the study timeline (Figure 26) and constraints of sampling logistics around unexpected storm events. It is more realistic to characterize the sum of geomorphology and habitat shifts between the two surveys as integrated changes over a multi-year scale. Despite the inherent uncertainty in the study, this benthic morphodynamic study is a valuable resource to inform management of publicly-valued and historically under-studied habitats.

**Figure 26.** Significant wave height (Hs) in meters (m) at offshore NDBC Buoy 44009 (38.461 N 74.703 W) between 2011–2016 including the passage of hurricanes Sandy and Joaquin near the study site.

Spanning four years and a total of 22 storms—including four hurricanes between the two overall site surveys is the prime challenge facing our interpretation of these results and direct comparison of pre- and post-Sandy (2014–2015) benthic assemblage data. Typically, benthic storm-response studies must be conducted fairly soon before and after an event, in order to minimize any effects from additional weather events [6]. However, in this project we were unable to begin benthic mapping until almost two years after Hurricane Sandy made landfall in the region. While it is in the interest of the National Park Service who manages Assateague island to gain any and all knowledge of the impacts of hurricane Sandy on subtidal habitats and resources along ASIS, it is more reasonable to look at the post-Sandy (2014–2015) data as the results of continuous storm action, of varying degrees of magnitude, over the entire study period from October, 2011 through October, 2015.

The passage and impact of hurricane Joaquin in September 2015 afforded us a unique opportunity to conduct targeted post-storm surveys in areas that we had mapped in May 2015. This provided the basis for examining storm induced changes. The resulting shifts illustrated in the sediment type contacts in the side-scan sonar (Figures 11 and 12) and in the seabed morphology revealed in the bathymetry (Figures 21 and 23) indicate significant shifts in the seabed ranging between 20–180 m horizontally with bathymetric changes of as much as 2 m with the passage of shore oblique bars. In general, the shifts were in a South to Southwest direction consistent with the impacts from strong northeasterly winds and waves that accompanied the storm.

In terms of recommendation for future research efforts we would recommend the development and implementation of a targeted rapid response mapping program. The development of a short-term program for immediate storm-response surveys during the hurricane season at Assateague—while likely constrained to smaller spatial scales—could provide more accurate links between storms and benthic morphodynamics. Such a program would necessarily work on a reduced scale, targeting a select number of sites that are determined to undergo the most change from storm impacts.

**Author Contributions:** The authors contributions were as follows. A.T. was responsible for overall project conceptualization, resources, writing of the draft and revisions along with project administration and funding acquisition. A.A. was responsible for GIS analysis on bathymetric and sediment texture change, assisted in writing the draft and preparing figures. K.H. was involved in field data collection, data analysis and data curation. C.D. assisted in data acquisition, data analysis and provided input and guidance to the draft and revisions.

**Funding:** This study was funded by the National Park Service under CESU Task P14AC00380 and modification 14-0001.

**Acknowledgments:** Pre-storm data was provided by Versar, Inc. We thank the Assateague Island National Seashore staff and project partners for their input, including Monique Lafrance-Bartley, Bill Hulslander, Charles Roman, Mark Finkbeiner, Sara Stevens and Robin Baranowski. Field and lab work was made possible by Captains Kevin Beam and Evan Falgowski, Ellie Rothermel, Hannah Rusch, Danielle Ferraro, Stephanie Dohner, Tim Pileguard, Trevor Metz, Anna Schutschkow, Annie Daw, Jason Button, Leah Morgan, Drew Friedrichs, Rachel Dalafave, Alec Halbruner, Meghan Owings, Rachel Oidtman, and Alex Itin. We thank the Environmental Protection Agency Region III's Ocean and Dredge Disposal Program Team, especially Renee Searfoss and Sherilyn Morgan Lau, for the loan of a Young grab sampler, as well as the University of Delaware's Chris Sommerfield and Adam Marsh for the loans of additional laboratory and field equipment.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2019 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/).

## *Article* **Application of Sentinel-2 Multispectral Data for Habitat Mapping of Pacific Islands: Palau Republic (Micronesia, Pacific Ocean)**

**Francesco Immordino 1, Mattia Barsanti 2, Elena Candigliota 1, Silvia Cocito 2, Ivana Delbono <sup>2</sup> and Andrea Peirano 2,\***


Received: 31 July 2019; Accepted: 31 August 2019; Published: 12 September 2019

**Abstract:** Sustainable and ecosystem-based marine spatial planning is a priority of Pacific Island countries basing their economy on marine resources. The urgency of management coral reef systems and associated coastal environments, threatened by the effects of climate change, require a detailed habitat mapping of the present status and a future monitoring of changes over time. Here, we present a remote sensing study using free available Sentinel-2 imagery for mapping at large scale the most sensible and high value habitats (corals, seagrasses, mangroves) of Palau Republic (Micronesia, Pacific Ocean), carried out without any sea truth validation. Remote sensing 'supervised' and 'unsupervised' classification methods applied to 2017 Sentinel-2 imagery with 10 m resolution together with comparisons with free ancillary data on web platform and available scientific literature were used to map mangrove, coral, and seagrass communities in the Palau Archipelago. This paper addresses the challenge of multispectral benthic mapping estimation using commercial software for preprocessing steps (ERDAS ATCOR) and for benthic classification (ENVI) on the base of satellite image analysis. The accuracy of the methods was tested comparing results with reference NOAA (National Oceanic and Atmospheric Administration, Silver Spring, MD, USA) habitat maps achieved through Ikonos and Quickbird imagery interpretation and sea-truth validations. Results showed how the proposed approach allowed an overall good classification of marine habitats, namely a good concordance of mangroves cover around Palau Archipelago with previous literature and a good identification of coastal habitats in two sites (barrier reef and coastal reef) with an accuracy of 39.8–56.8%, suitable for survey and monitoring of most sensible habitats in tropical remote islands.

**Keywords:** Sentinel-2; Remote Sensing; habitat mapping; mangroves; coral reefs; climate change; vulnerable habitats

#### **1. Introduction**

Coral reefs, seagrasses, and mangroves are threatened worldwide by climate change, whose main effects are sea temperature increase and ocean acidification [1–3]. In addition, the frequency of discrete extreme warming events (heat waves) threatening global biodiversity has increased, with projections indicating they will become more frequent [4]. Moreover, coastal coral ecosystems are threatened by additional anthropogenic pressures as overfishing, urbanization, and tourism [5]. Consequences are producing concerns for the maintenance of ecosystem functioning and the associated flow of services that coral reefs provide.

Having a deep cultural heritage for ocean conservation, the Pacific Ocean countries are strong advocates of a 'Blue Economy Strategy' and the sustainable use of ocean resources for economic growth, and are turning towards an increased reliance on green tourism. Palau Republic, among others, has emerged as a global leader in ocean conservation, receiving the 2012 Future Policy Award for developing the world's best policies to protect oceans and coasts. People living in Pacific Islands depend on healthy coastal ecosystems for their survival. Ecosystems such as coral reefs, mangroves, and seagrass beds favor coastal protection, provide food, building materials, and they represent the principal economic incomes for fishing and tourism industries [5]. Hermatypic corals are the most sensitive organisms to the synergic effects of warming, hurricane destruction, and ocean acidification [1]. Future climate scenery predicts that over the coming decades coral mortality may reach up to 60% in the areas where shallow coral reefs are present [6], driving to the elimination of most warm water coral reefs by 2040–2050 [1]. Hence, the preservation of global reef biodiversity, the monitoring of climate change effects on reef and islands together with the development of global strategies to reduce greenhouse gas emissions are among the major management issues to counter the effects of climate change [7].

The need to effectively manage coral reef systems and their associated coastal environments, such as mangroves and seagrasses, requires the ability to document their present status and monitor changes over time. Benthic habitat mapping of coastal ecosystems is an important and essential mean to provide marine resource assessments for coastal management and ecological analysis. Habitat mapping by remote sensing allows large scale environmental patterns and is highly cost-effective compared to the sampling of physical areas achievable by field survey, which provide accurate data but at highly detailed scales [8–13]. RS imagery were used in conjunction with state-of-the-art RS algorithms to map reefs geomorphology and habitat distribution [14–17].

Today, widely available orthorectified satellite imagery (Google Earth) and rapid development and cost reduction in GIS, multibeam echosounders, Lidar, and RS technologies are making large-scale morphometric quantification of reefs feasible [18–20].

In 2001, the "Millennium Coral Reef Mapping Project—Understanding, Classifying and Mapping Coral Reef Structures Worldwide Using High Resolution Remote Sensing Spaceborne Images" examined ~1500 images to design a thematically rich (966 classes) geomorphological classification scheme, used to interpret and map every single reef of the planet. Distributed as Geographic Information System (GIS) layers in late 2003, the map products have been used for a variety of applications, from establishment of marine protected areas in Papua New Guinea and Eastern Caribbean, reef condition assessment in the Caribbean, morphometric analyses of Maldivian reefs, and geochemical budgets in French Polynesian atolls [21].

Sentinel mission is one of the last RS programs dedicated to the study of marine and coastal environment. A constellation of two satellites, Sentinel-2 A and B, was launched in 2015 by European Space Agency (ESA) as part of the Copernicus Program. The Sentinel-2 satellites carries an innovative wide swath high-resolution multi-spectral imagery with 13 bands (443–2190 nm) in the visible, near infrared, and short wave infrared part of the spectrum. The combination of high spatial resolution of 10–20–60 m, novel spectral capabilities, a swath width of 290 km, a global coverage of land surfaces from 56◦ S to 84◦ N and frequent revisit times (5–15 days) provides unprecedented views of Earth. Sentinel mission not only offers continuity of services for the moderate resolution multispectral Spot XS and Landsat Thematic Mapper series sensors, but it also has several technical improvements that may lead to enhanced capability in coral reef mapping applications [22–24].

The aim of this paper is to enhance the effectiveness of RS repositories as a powerful tool for coastal resources assessment and management of remote Pacific islands, where mapping data are missing or lacking at present. A low-cost methodological approach for a preliminary habitat classification useful for coral coastal management purposes is proposed for Palau Archipelago (Micronesia). Different available resources such as free satellite images from the Copernicus Open Access HUB services and ancillary information as Google Earth imagery, scientific reports and literature were used to create maps of relevant habitats without field data validation. Results on mangrove forest distribution around the islands and on habitat classes in two sites (barrier reef and coastal reef) are compared with available detailed habitat maps of the Palau coastal area, and pros and cons of this approach are discussed.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The Republic of Palau is an island country located in the western Pacific Ocean; the country contains approximately 340 islands, forming the western chain of the Caroline Islands in Micronesia, and has an area of around 466 square kilometres (Figure 1).

**Figure 1.** The main island of the Palau Republic (Micronesia—West Pacific). Numbers indicate the barrier reef zone (site no. 1) and coastal reef zone (site no. 2) where the shallow-water benthic habitat mapping was performed (basemap and inset Google © 2019).

The Palauan coral reef ecosystem has the most diverse flora and fauna of Micronesia. Palau has some of the most extensive seagrass beds in all of Micronesia, hosting 10 species of seagrass [25]. Seagrasses are valuable habitats that provide important ecological components of coastal ecosystems worldwide. Moreover, one of the most significant ecosystems in Palau are the mangrove forests, the transition zone between terrestrial and marine ecosystems. The most extensive areas of mangrove forests occur along the west coast of the main island of Palau Republic (Babeldaob), covering approximately 80 percent of the shoreline. In the mangrove forest of Palau, there are 18 mangrove trees and associated plant species, the most diverse in Micronesia. These habitats provide a number of ecological functions, from nurseries for juvenile fish to food and shelter for invertebrates and rare, protected species as sea turtles, crocodiles, and dugongs. Mangroves are ecologically important also because they help stabilize coastal areas by trapping and holding sediments washed down from inland areas and local watersheds. Moreover, the islands of the Pacific region lie in one of most seismically active regions of the world and for the coastal population mangrove forests are often the first line of defence against such natural calamities.

#### *2.2. Image Processing and Shallow-Water Benthic Classification*

Two Sentinel-2 2017 satellite images (4 July 2017) were downloaded from the open source portal of the European ESA-Copernicus project (https://scihub.copernicus.eu/). Sentinel-2 product used in this work provides orthorectified Top-Of-Atmosphere (TOA) reflectance (Level 1-C), with sub-pixel multispectral registration. Cloud and land/water masks are included in the product. SENTINEL-2 products are available to users in SENTINEL-SAFE (Standard Archive Format for Europe) format, including image data in JPEG2000 format, quality indicators, auxiliary data, and metadata (https://sentinel.esa.int/web/sentinel/user-guides/SENTINEL-2-msi/data-formats). The SAFE format has been designed to act as a common format for archiving and conveying data within ESA Earth Observation archiving facilities. The SAFE format wraps a folder containing image data in a binary data format and product metadata in XML. This flexibility allows the format to be scalable enough to represent all levels of SENTINEL products. The SENTINEL-2 product refers to a directory folder that contains a collection of information and includes: manifest.safe file which holds the general product information in XML; preview image in JPEG2000 format; subfolders for measurement datasets including image data (granules/tiles) in GML-JPEG2000 format; subfolders for datastrip level information; subfolder with auxiliary data (e.g., International Earth Rotation & Reference Systems (IERS) bulletin); HTML previews.

The ERDAS ATCOR radiative transfer model was used for atmospheric corrections (Figure 2): it eliminates atmospheric and topographic effects in satellite imagery and extracts physical surface properties, such as surface reflectance, emissivity, and temperature. ATCOR Workflow employs a database containing the result of radiative transfer calculations based on MODTRAN® 5.

**Figure 2.** Atmospheric correction steps: ERDAS ATCOR radiative transfer model used for Sentinel-2 L1C product. \* Tropical corresponds to a water vapor column of 4.11 cm at sea level, while \*\* Maritime represents the aerosol conditions in areas close to the sea.

The processed multispectral satellite images showed the presence of calm, clear waters, and—considering the shallow depths—water penetrating correction and sun glint corrections procedures were not applied [26,27]. To this end, standard atmospheric corrections have been carried out, considering the type of atmospheric column, the SENTINEL-2 image acquisition period and the sensor instrument parameters that the ERDAS software automatically recognizes by reading the metadata at the beginning of the acquisition procedure.

All Sentinel-2 bands were processed and 're-sampled' in order to get the highest resolution possible (10 m); then the ENVI image analysis software was used for classification procedures.

Two separate approaches were used to classify inland mangrove cover and habitat classes related to coral platform; the analysis processing chain is illustrated in Figure 3.

**Figure 3.** Workflow diagram illustrating the steps followed for the habitat mapping of Palau.

First, for the 'coastline emergent vegetation' or mangrove classification, all Sentinel-2 13 bands were processed and 'resampled', building a metafile, in order to get the highest resolution possible (10 m). Then Band 11 (wavelength range 1542–1685 nm) was particularly enhanced in order to distinguish the land, habitat of mangroves, and sea (Figure 4).

**Figure 4.** Sentinel-2 satellite image on the Republic of Palau on 4 July 2017, with a 10 m spatial resolution. Enhanced Band 11 was used, showing the land (white), the mangrove cover (grey), the sea (black).

The spectral signatures of mangroves and land vegetation were extracted from the Sentinel-2 image and here reported to show their spectral separability (Figure 5).

**Figure 5.** Spectral signatures of mangroves and terrestrial vegetation in Sentinel-2 spectral bands.

Hence, the ENVI 'supervised classification' through 'region of interest' (ROI) and 'maximum likelihood' methods were applied in order to separate the land from mangroves and sea. This approach allowed to compute the total shoreline length and the mangroves cover for the whole Palau Archipelago (Figure 6).

**Figure 6.** Sentinel-2 satellite image (Band 11) on the Republic of Palau on 4 July 2017, with a 10 m spatial resolution: the ROIs in the process of a 'supervised classification' (on the left); the classification result (on the right) where land is orange, the mangrove cover is green, the sea is blue. Clouds are isolated and labelled as unclassified.

For the habitat classification related to coral platforms, we restricted the analyses to two sites, located in the western side of the main Palau island (Figure 1) and representatives of the two main

morphological zones of the island: the barrier reef (site no. 1) and the coastal reef (site no. 2). The ENVI algorithm that showed better results was the 'unsupervised isodata classification' performed for a maximum of 20 classes on Sentinel-2 using bands 2-3-4-8 with 10 m of resolution. Unsupervised classification yields an output image in which a number of classes are identified and each pixel is assigned to a class [28]. This classification often results in too many land cover classes, particularly for heterogeneous land cover types, and classes often need to be combined to create a meaningful map; the classification is useful when there is no pre-existing field data or detailed aerial photographs for the image area, and the user cannot accurately specify training areas of known cover type. In marine habitat classification, the unsupervised isodata classification is considered the most appropriate for an approach with no sea truth validation [29–31].

The spectral signatures of the reef coverage were extracted from Sentinel-2 image and here reported on a plot (Figure 7) to show the spectral separability among them. To check the reliability of atmospheric corrections and subsequent classification of the Sentinel-2 image, the ROI of the coverage classes present in the coralligenous environment was extracted.

**Figure 7.** Spectral signatures of reef coverage (the barrier reef zone, site no. 1) in ERDAS software.

The image interpretation for the coastal benthic classification was based on two main components: the geomorphologic description of the seabed and the biological relevance of the associations living in the zone (Table 1).

We referred to the NOAA (National Oceanic and Atmospheric Administration, Silver Spring, MD, USA) Shallow-Water Benthic Habitats Manual [32] to identify priority habitats in terms of ecological function and coral presence. Hence, on the barrier reef, the following five main habitat classes were used for RS classification:

(1) The 'fore reef' is the outward part of the reef barrier. On its underwater cliff, all coral biodiversity is concentrated along the first 20–40 m of depth. It represents the most important reservoir for coral maintenance and survival; (2) The 'reef crest' is the part of the outer barrier reef more exposed to open-ocean waves. Its associated subclasses are the living coral and one algal ridge formed mainly by coralline algae. It is the most important area of the barrier for the defence of the coastline and it is the most sensitive to mortality due to low tides, elevated seawater temperatures, and storms; (3) The 'back reef' is the area of the barrier reef formed by a coral platform and is limited towards the coastline by the lagoon. It may be very large and in its shallow areas is characterized by an eroded platform and rubbles with associated subclasses of coralline algae, massive corals, and algal turf. It is the area where fragments of corals of the reef crest damaged by open-ocean waves may survive. The slope of the back reef limited by the lagoon is normally formed by sand and coral rubbles; (4) The 'patch reef' includes as subclasses the scattered living coral formations like coral knob, coral head, irregular little islands of aggregated corals surrounded by sand, found on the back reef, on coral platform, or dispersed in the

lagoon. The patch reef is important because it includes isolated living reef that represent a 'reservoir' both in term of larvae and individuals of different coral species; (5) The 'lagoon' is the area between the back reef and the coastline. Here coral, sand, and rubble are abundant. However, the lagoon represents a limit for the RS investigations when the depth may exceed 30 m. In this deep area, indicated as 'deep lagoon', the coral cover may be important but should be investigated with other methodologies (transects, multibeam, side scan sonar). The species of coral inhabiting the lagoon and the deep lagoon could be an important reservoir of larvae for the surrounding communities.


**Table 1.** Geomorphological description and biological importance of the main classes/subclasses used for the Palau habitat mapping with Sentinel-2 imagery.

On the coastal reef, the following three main habitat classes were used:

(1) The 'coastal fringing reef' or 'coastal reef crest' is the seaward fringe of the coastal reef flat. Its associated subclasses are living corals in good health; this habitat class is subjected to erosion by waves and it is the most important indicator of the coastal reef erosion together with blue holes; (2) The 'coastal reef flat' includes the reef platform formed by dead coral surrounding the coast of the major island and it can reach the coastline as a rock substrate. Its subclasses include seagrasses and algae (mainly macroalgae) in the area where fine sediments are deposited. Its importance is related to the coastal defence from erosion; (3) Shoreline and emergent vegetation is formed from the emergent vegetation habitat composed primarily of red mangroves, generally found in areas which are sheltered from high-energy waves.

Finally, the classes 'unknown' or 'unclassified' represent the uninterpretable areas due to cloud shadow, water depth, or other interference.

#### *2.3. Validation of Image Interpretation*

To validate the Sentinel-2 image interpretations, the habitat mapping results were compared with available Palau NOAA maps (https://products.coastalscience.noaa.gov/collections/benthic/e102palau/). These maps are the final results of a 4-year project [33]: NOAA collected 2002 and 2004 multispectral Ikonos and supplemental Quickbird satellite imagery of Palau Archipelago with spatial resolution of 4 m (raw multispectral) and 1 m (pan sharpened). Color balanced imagery proved suitable for visual

extraction of the habitat classes. NOAA image process included atmosphere correction, deglinting, color balancing, orthorectification, correction for water column effects and pan sharpening. Collection constraints were set to control environmental effects such as glare, glint and other interferences that would limit visualization of benthic features. NOAA multiple collects were conducted to mosaic multiple scenes up to a maximum of 10% cloud cover. These images were used by NOAA to manually interpret and delineate geomorphologic features, cover types, and habitat boundaries. NOAA maps were produced following the schemes for habitat classification prepared by coral reef biologists and mapping experts. Ground validation information was used to investigate uncertainties on the photo-interpreter behalf during the decision-making process of the manual delineation of zones or structures. To test the interpretation accuracy, sea-truth validation was performed by NOAA through 623 benthic habitat characterizations conducted in four areas of Palau. Results showed an overall accuracy in habitat identification from 88.4 to 97.3% and with a thematic accuracy of the habitat classification schemes ranging from 0.79 Tau for detailed cover (79.9% overall accuracy) to 0.96 Tau for major structure (97.3% overall accuracy).

In the present work, main habitat classes cover derived from Sentinel-2 data were compared with NOAA and used as reference data for accuracy evaluation of the proposed method.

#### **3. Results**

The shoreline length and the mangrove cover were calculated for the whole Palau Archipelago. A comparison between NOAA results [28] on 2003–2006 Ikonos imagery and the present work on 2017 Sentinel-2 images is shown in Table 2. Shoreline and mangroves area values show relevant differences in the considered period, equal to 18.8% and 25.0%, respectively. The shoreline length of the Republic of Palau is calculated by NOAA, derived from 2003–2006 Ikonos Imagery, as 1021 km by visual interpretation and manual delineation of satellite imagery.

**Table 2.** Length of shoreline (km) and mangrove areas (ha) on Palau Archipelago computed in the present work (ENEA) and by NOAA (Analytical Laboratories of Hawaii, 2007). Differences are in percentages (%).


Regarding the mangrove habitat, thanks to the Band 11 of the Sentinel-2 imagery, the mangrove cover for the whole Palau Archipelago was determined with a resolution of 10 m and with a great accuracy. The optimum classification is mathematically confirmed by the Jeffries-Matusita distances among classes: the output is shown in Table 3. Good coefficients range between 1.999 and 2.000 [34].

**Table 3.** ROI matrix separability, for the SENTINEL-2 imagery, using Jeffries-Matusita algorithm.


As well as mangroves, also for coastal benthic habitats, Jeffries-Matusita (JM) spectral separability coefficients (Table S1, Supplementary Materials) confirm values between 1.9–2.0 among the different classes, with the only exception for the combination back–reef sand–algae and patch reef coral–sand with a JM value of 1.8, probably linked to the sedimentological features and texture of these classes due to the presence of biodetritic sand in the two classes. The image processing steps confirm the classes

spectral reflectance, since bottom-types are statistically separable and identifiable on the base of their reflectance spectra [35]. The scatter plot shows a high spectral separability among coral reef classes (Figure 8).

**Figure 8.** Scatter plot of the coral reef class coverage, showing a high classes separability.

Unsupervised isodata classification allowed the identification of 12 habitat classes/subclasses for the barrier reef (site no. 1) and 10 classes/subclasses for the coastal reef (site no. 2) (Figures 9 and 10). To test the goodness of the habitat classification performed on Sentinel-2 images, data were compared with NOAA habitat maps from Ikonos satellite imagery and sea-truth validations carried out in the period 2002–2004 [33]. Comparison of habitat classifications on the barrier reef (Figure 9) shows how the classification with Sentinel-2 image was not able to identify the class bank/shelf escarpment, i.e., the deeper area of the outward barrier found on NOAA images. In Sentinel-2 classification, the whole area was identified as fore reef. The reef crest and the back reef showed differences between NOAA findings [33] and our results. The area of reef crest was thinner and better defined in NOAA maps, wider in our estimation. On the contrary, the back reef and the patch reef areas were more defined in Sentinel-2 imagery where both scattered corals and patch reef were identified as well as the nature of the bottom (sandy or hard bottom).

**Figure 9.** Comparison between habitat classifications in the barrier reef zone (site no. 1, see Figure 1). On the left, NOAA's map and 2004 Ikonos imagery (informative layers downloaded from https: //products.coastalscience.noaa.gov/collections/benthic/e102palau/); on the right, present work (ENEA) habitat map and 2017 Sentinel-2 imagery.

**Figure 10.** Comparison between habitat classifications in the Coastal Reef zone (site no. 2, see Figure 1). On the left, NOAA's map and 2004 Ikonos imagery (informative layers downloaded from https://products.coastalscience.noaa.gov/collections/benthic/e102palau/); on the right, present work (ENEA) habitat map and 2017 Sentinel-2 imagery.

The extension in hectares (ha) of the three main barrier reef areas (Table 4) showed an overestimation of Sentinel-2 classification of 12.3% for bank/shelf–fore reef and 20.0% for reef crest and back reef areas, differently an underestimation of lagoon area of 20.3%.

**Table 4.** Palau: Habitat classes cover (ha) in the two zones of the barrier reef (site no. 1) and coastal reef (site no. 2) computed for the 2017 in the present work (ENEA) and by NOAA (Analytical Laboratories of Hawaii, 2007) in 2004. Differences in accuracy are in percent (%).


The comparison for the coastal reef area of Palau (Figure 10) showed the suitability of Sentinel-2 images for the recognition of the main habitat classes on the reef platform found by NOAA in 2004. Both the reef crest (or fringing reef) and reef flat were recognized in Sentinel-2 images and classification.

The extension in hectares (ha) of the three main coastal reef areas (Table 4) showed an underestimation of Sentinel-2 classification of 33.9% for the reef crest and an overestimation of reef flat of 34.1%. Differently, the seagrasses were clearly identified by the classification and the two sub-classes showed differences in density or species. In this case clear differences in the spatial extent of the seagrass beds were found in comparison with NOAA data back to nearly 15 years.

#### **4. Discussion**

#### *4.1. Shoreline, Mangroves, Coral Reefs, and Seagrasses*

The difference in shoreline length (Table 2 in results paragraph) is due to the diverse techniques applied for shoreline extraction: manual digitalization (NOAA) and the use of classification algorithms by satellite imagery (ENEA). Furthermore, the wide range of shoreline values could be produced from the different nominal scale of source maps applied in order to extract the coastline [36].

Regarding the mangrove habitat, notwithstanding the different methodologies or imagery applied in NOAA and ENEA studies, the mangrove cover comparison with previous data is interesting. The difference in the extension of the mangroves from 2003–2006 imagery to 2017 is evidenced, but the disparity of NOAA and ENEA habitat mapping methods imposes caution when interpreting this result. The difference equal to 1100 ha between NOAA and ENEA (Table 2 in results paragraph), meaning an increase of the mangrove cover equal to 25.0% in the last decade (2006–2017), is not convincing. The difference between NOAA and ENEA mangrove habitat mapping is likely due to gaps in previous maps. Diachronic maps allow the measurements of temporal changes through concordance and discordance maps. Figure 11 shows a focus of this evidence in the eastern side of the main island of Palau, where it is possible to observe a great concordance between ENEA and NOAA mangrove cover (light green, Figure 11) but also significative discordances like a wide mangrove cover identified only by the present work (dark green, Figure 11) and a very limited mangrove cover identified only by NOAA.

**Figure 11.** (**a**) A detail of the east sector of the Republic of Palau, with Sentinel-2 enhanced Band 11, where dark grey is the mangrove cover; (**b**) mangrove cover diachronic map showing a concordance wide area of mangroves identified by ENEA and NOAA (light green); two discordance areas with ENEA mangrove cover only (dark green) and NOAA mangrove cover only (yellow), respectively.

On the basis of these results, it is important to note that RS comparisons with different satellite imageries and methods could have some embedded drawbacks. Moreover, the present work's mangrove cover, equal to 5500 ha (Table 2), is in good agreement with 2011 estimations [37]. These latter were achieved with Landsat Enhanced Thematic Mapper (ETM+) data, with images collected during epoch centered on the year 2000, along with very high resolution images such as Ikonos and Quickbird. The overall estimation of 5666 ha for Palau Archipelago [37] is nearly coincident with our results (5500 ha) suggesting an overall equilibrium in Palau Archipelago in the last decades (only 3% difference in hectares in around 17 years). For these reasons an overall equilibrium in the mangrove forests is considered reliable.

With regard to coral reefs and seagrasses, the classification analysis demonstrates that the coral bottom-types are statistically separable and identifiable based on their reflectance spectra. We reason that features in reflectance arise primarily as a result of spectral absorption processes. Radiative transfer modeling shows that in typically clear coral reef waters, dark substrates such as corals have a depth-of-detection limit on the order of 10–20 m [35].

To measure the accuracy of the method proposed on barrier reef and coastal reef zones, we compared results of habitat mapping cover data of main classes from Sentinel-2 image versus data from NOAA maps assumed as reference data. The accuracy varied from 44.3 to 54.6% for the barrier reef estimations and from 39.8 to 56.8% for the coastal reef estimations (Table 4). These results are in the range of overall accuracy estimated for large area studies such as those performed for the Great Barrier Reef (GBR), where overall accuracy estimations of user/producer on Landsat 8 OLI classification for geomorphic zonation and benthic cover types were respectively 59.5 and 32.8% [24]. On the Southern Section Cairn GBR management region, the Sentinel-2 benthic mapping analysis performed on eight reefs with sea truth data validations produced an overall accuracy of 49% using six categories [32].

Reef flat subclasses showed differences between our and NOAA results (Figure 10): two areas were identified by unsupervised classification, but it was impossible to attribute them to turf or algae with different cover. The seagrasses were clearly identified and the cover showed a variation and two further sub-classes probably related to different density or different species were found.

#### *4.2. Limits and Challenges of the Present Approach vs. Similar Studies*

In the present work, Palau Island habitat classes/subclasses were chosen with the objective to achieve a RS habitat map centred on the most ecologically important habitats sensitive to climate changes in the Pacific region, like coral reefs, seagrasses, and mangroves. Habitat classes were combined in a unique geomorphological and ecological classification, a method considered as the most appropriate for RS of tropical coastal areas [29]. The unsupervised classification used for benthic habitat classification exploits the advantages of statistical segmentation to find natural boundaries in a dataset and provides consistent classification at multiple sites with little to no ground truth required [30]. However, without a ground validation, it was not possible to map turf, coralline algae, macroalgae, and seagrasses at species level. Differently, the mangrove cover along the island coastline has been identified with high detail, in relation to the high radiometry of the Sentinel-2 spectral signals (14 bit). The method showed its potentiality for medium (seasonal) and long-time (decadal) assessment of changes in sensible communities like seagrasses, also in conjunction with past Landsat 8 imagery collections, but it is limited to main habitat reconnaissance of large-scale habitat mapping [31].

Sensitivity of coral reefs, seagrasses, and mangroves to any effect of climate change varies for the three habitats. It is difficult to separate any effects of climate change from those produced by coastal development and land use practices. In the past six decades, the Palau Archipelago experienced in 1998 and 2010 two episodes of extreme heat associated with El Ni ˇno, considered conducive for coral bleaching [38]. Bleaching was most severe in the north-western lagoon, whereas in the bays where temperatures were higher than elsewhere, bleaching and mortality were low. Although the intensity and duration of elevated temperatures are strong predictors of a coral's fate, affecting survival and reproduction [4], the extent of bleaching is different not only by site but it can affect also deep reef. For this reason, the survey of shallow water coral habitats with RS assumes a greater importance as they are source for new propagules also for deeper outer slope coral communities [39]. Seagrasses are likely to be highly sensitive to increases in sea surface temperature, whether they occur as short-term spikes over periods of hours or as chronic exposures for weeks or months. Temperature extremes are known to reduce seagrass growth and lead to plant mortality with reduced carbon sequestration and habitat alteration [4,40]. Increases in seawater temperature also translate into increased disturbance via stronger extreme weather events such as heavy rainfall and tropical cyclones and storms, which put additional stress on seagrass habitats. Direct climate change impacts on mangrove ecosystems are likely to be less significant than the effects of associated sea level rise [2,41]. Sensitivity of mangroves

to increased sea surface temperature is likely to be moderate. Rise in temperature and the direct effects of increased CO2 levels are likely to increase mangrove productivity, change the timing of flowering and fruiting, and expand the ranges of mangrove species into higher latitudes.

#### **5. Conclusions**

The proposed approach showed how nowadays the availability of different resources such as free Sentinel-2 imagery and ancillary information allows a reliable classification of marine habitats and the quantification of high value tropical habitats colonized by coral, seagrasses, and mangroves. The accuracy of the method and a revisiting period of the Sentinel-2 satellites of 5–15 days offer the possibility to follow communities covers from season to season and to assess environmental changes over time. Although the value of ground validation is evident to disentangle some uncertainty of interpretation, we demonstrate that the proposed approach is appropriate for extensive large-scale habitat classifications in remote sites like Palau Republic and all Pacific islands. Furthermore, the present methodology can be a good base for future monitoring programmes to be conducted also by resident personnel trained in studies for Marine Spatial Planning and Marine Protected Areas.

Finally, the present work demonstrates the chance offered by free availability of imagery and information to optimize time and resources through worldwide collaboration of research teams to mitigate the effects of climate change in remote Pacific islands.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2077-1312/7/9/316/s1, Table S1: Jeffries-Matusita spectral separability coefficients: values are between 1.9–2.0 among the different classes, except the combination "back–reef sand–algae and patch reef coral–sand" with a JM value of 1.8.

**Author Contributions:** Conceptualization, F.I. and E.C.; Methodology, F.I., A.P., M.B., I.D., E.C., and S.C.; Formal analysis, F.I., E.C., A.P., M.B., I.D., and S.C.; Writing—original draft preparation, A.P., I.D., S.C., and M.B.; Writing—review and editing, A.P., I.D., S.C., and F.I.

**Funding:** This research received no external funding.

**Acknowledgments:** We thank N.M. Caminiti for supporting our work, and D. Dominici S. Petric, E. Wang and three anonymous referees whose suggestions greatly improved the manuscript.

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

#### **References**


© 2019 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/).

## *Article* **Statistical Deviations in Shoreline Detection Obtained with Direct and Remote Observations**

**Giovanni Pugliano 1, Umberto Robustelli 1,\*, Diana Di Luccio 2,\*, Luigi Mucerino 3, Guido Benassai <sup>1</sup> and Raffaele Montella <sup>2</sup>**


Received: 29 March 2019; Accepted: 6 May 2019; Published: 11 May 2019

**Abstract:** Remote video imagery is widely used for shoreline detection, which plays a fundamental role in geomorphological studies and in risk assessment, but, up to now, few measurements of accuracy have been undertaken. In this paper, the comparison of video-based and GPS-derived shoreline measurements was performed on a sandy micro-tidal beach located in Italy (central Tyrrhenian Sea). The GPS survey was performed using a single frequency, code, and carrier phase receiver as a rover. Raw measurements have been post-processed by using a carrier-based positioning algorithm. The comparison between video camera and DGPS coastline has been carried out on the whole beach, measuring the error as the deviation from the DGPS line computed along the normal to the DGPS itself. The deviations between the two dataset were examined in order to establish possible spatial dependence on video camera point of view and on beach slope in the intertidal zone. The results revealed that, generally, the error increased with the distance from the acquisition system and with the wash up length (inversely proportional to the beach slope).

**Keywords:** DGPS measurements; video camera observation; shoreline position; beach survey

#### **1. Introduction**

Coastal areas are highly dynamic environments that provide important benefits, but are also subject to a variety of natural hazards such as beach erosion, tsunamis [1,2], and floods [3]. For coastal zone monitoring and coastal risk assessment [4], shoreline detection is a fundamental work. The shoreline is the line where land meets the sea, and due to the dynamic nature of the sea, sometimes it is difficult to determine a precise line that can be called the "shoreline". As reported by Boak and Turner [5], a functional definition of the "shoreline", which has to consider the shoreline in both a temporal and spatial sense, is required.

Many authors highlighted the existence of different methodologies for coastal monitoring [6–9], not only limited to shoreline detection, based on direct and remote acquisition systems.

Direct shoreline surveys are normally carried out using the GPS technique by post-processing or real-time methods [10]. The main limitation of this method is the huge time required for covering large stretches of the coastline and the difficulties inherent in doing ad hoc timely post-storm measurements.

Remote shoreline observations can be distinguished in remote sensing [11–13], UAV (Unmanned Aerial Vehicles) [14], and video monitoring [15,16], which were first introduced by Aarninkhof [17] and Turner et al. [18]. These remote observations have been also extensively used to validate wind and wave numerical models [19–23] and to outline the environmental big picture in marine spatial planning applications [24–26].

In addition to using UAVs in beach survey [27,28], video monitoring can provide a remotely-sensed measurement, fixed at a secure location (e.g., a tower or high-point), with the capability of acquiring imagery at a frequency ranging from fractions of seconds to hours. The technology is relatively low-cost, but the main issue is the processing method, especially the image rectification process, considering that the imagery is strongly oblique and relies on a number of GCPs (Ground Control Points) for finding the best geometry solutions. This technique has been successfully applied for both shoreline monitoring [29,30], rip current measurement [14] and morphodynamics classification of sandy beach [31,32].

Several techniques are proposed to extract shoreline from images, based on discriminating sea from sand. Plant and Holman [33] used a method initially developed for gray-scale cameras, called Shoreline Intensity Maximum (SLIM), which defines the shoreline as the cross-shore position at which wave breaking is maximized and consequently corresponds to a maximum in pixel intensity. With the adoption of color cameras, spectral information could also be used to identify the shoreline, making use of the water property to absorb the Red signal (R) and the sand property to absorb the Green (G) and blue signal [34,35]. Following Almar et al. [36], we identify beach pixels as channel values with a high R/G ratio, whereas water pixels as high G channel values. The shoreline represents the transition zone between the two peaks. In order to smooth out high-frequency signals that are caused by disturbances like foam, we used time-exposure images [17,18].

In this paper, we performed shoreline detection by means of remote video camera observations and DGPS direct shoreline measurements on a sandy microtidal beach located in central Tyrrhenian Sea. Almar et al. [36], among other authors, have presented video-based methods to determine shoreline position in different tidal and wave energy contexts, showing an error dependence on local swash length, which is inversely proportional to the beach slope. Following this research item, we first performed a preliminary analysis of the GPS survey, which was followed by the evaluation of the statistical deviations of the video camera coastline measurements from the DGPS ones. The errors have been spatially processed, in order to examine if systematic deviations can be ascribed only to the intertidal beach slope or to additional factors, like the distance from the video camera. In particular, we examined the longitudinal and cross-shore errors associated with the video camera point of view.

The paper is structured as follows: after a brief description of the study area (Section 2), we illustrate the survey and validation methods used in Section 3. The results with the relative discussion are reported in Sections 5 and 6, respectively.

#### **2. Study Area**

The study area includes the beach of Serapo, located in the Central Thyrrenian Sea (Italy), which exhibits a rather regular morphology and rectilinear shoreline that is complexly long, approximately 1.45 km (Figure 1), aligned in the NW-SE direction (for the statistical processing, we considered only a 1.34 km-long shoreline, directly sampled with GPS). This beach is characterized by homogeneous grain size features with median diameter around 0.38 mm, indicating medium sands [37]. The topographic survey data, acquired in September 2017, refer to coastline and to five beach profiles T1–T5 as reported in Table 1. Coastline investigations were conducted in September 2018 (for the location, see Figure 1). Based on the difference in beach width and berm height, the studied beach can be subdivided into the following two stretches: a western beach stretch including profiles T1 and T2, which is 45–60 m wide and characterized by lower berm height (1.10–1.30 m) and a mean emerged beach slope of 5.5%; an eastern beach stretch including profiles T3, T4, and T5, which is characterized by higher beach width of 74–87 m and higher berm height of 1.29–1.46 m, with a lower emerged beach slope of 1.48–3.94%. This different beach width was already explained by Di Luccio et al. [37] in terms of partial clockwise beach rotation around a central pivotal point.

**Figure 1.** Serapo beach study area (**a**) with the location of the beach video camera system (red star) and the five investigated cross-shore transects (red lines), obtained by a beach survey conducted in September 2017; the beach profile along the transect with the location of the anthropic structures (black trapezium), Mean Sea Level (MSL), and the Mean Spring Tidal Range (MHSW), extracted by the official Italian tide archives (http://www.mareografico.it, last access: 30 April 2019) was reported for T1, T2, T3, T4, and T5. (basemap c 2018 Google product).

**Table 1.** Summary of Serapo beach's morphological characterization. The reported parameters are relative to the beach cross-shore profiles T1–T5.


#### **3. Materials and Methods**

#### *3.1. Sea Waves' Analysis*

The wave climate is connected to the methodology used for the definition of the coastline since the morphological changes affecting the beach after a sea storm can change the longitudinal beach profile (beach slope) in a more or less accentuated manner, influencing, as we shall see then, the remote acquisition processes. Moreover, if storm patterns and wave distributions (time scales) change, the coastline shape can evolve with erosion/accretion processes [38], inducing possible beach rotation, well evident with continuous video monitoring systems.

In order to define the sea condition characterizing the study area, we analyzed the historical time series extracted by the Ponza buoy 3-hourly data, provided by the Italian Sea Wave Measurement Network [39]. This directional pitch-roll buoy is moored in deep water [40], a few miles to the south of Ponza island (40◦52 00.10 N, 12◦56 60.00 E). The available parameters are significant wave height Hs, mean wave period Tm, and mean wave direction Dm, covering the years 1989–2014 (with some data gaps). Climatological wave analysis is presented, highlighting the wave storms that exposed a lower limit threshold of Hsequal to 2 m.

#### *3.2. Kinematic GPS Survey*

The Global Positioning System (GPS) surveying was employed to get the reference shoreline position. The collection of shore-parallel GPS positions was carried out in September 2017 using a single-frequency code and carrier phase receiver (Trimble Pathfinder ProXH) as the rover. The most common shoreline detection technique applied to visibly-discernible shoreline features is manual visual interpretation in the field, as reported by Boak and Turner [5]. In this paper, the shoreline was compiled by interpolating between described shore-normal beach profiles and a series of points collected along the beach face, which included maximum run-up limit and 0.5-m depth.

In order to compute the positioning solutions by various modes including single-point and relative positioning, post-processing of the raw data was performed using the RTKLIB software [41], an open-source software for standard and precise positioning [41]. The processing option for the relative positioning mode was set to carrier-based kinematic positioning; the ambiguity resolution was performed both recalculating the phase bias estimates every epoch (instantaneous) and using the Kalman filter to estimate the phase biases over many epochs (continuous). The last solution was much less noisy, as expected, and was adopted as the reference line. The analysis of the single point and carrier-based instantaneous solutions was included as well. In particular, the use of the single-point positioning was investigated as having the potential for being a low-cost method for shoreline mapping. In order to improve the results obtainable with this technique, the researchers have directed their studies towards the identification and reduction of multipath [42–44] as it is the highest source of error in the single-point positioning. In the last few months, smart phones equipped with a dual frequency multi-constellation GNSS receiver have appeared on the market. They could represent a valid and low-cost alternative to carry out the survey as shown by Robustelli et al. [45]. The position accuracy was evaluated on the basis of the calculation of the distance and azimuth of the estimated horizontal position error vector (Figure 2).

**Figure 2.** Scheme of the quantification of errors in the kinematic GPS horizontal positions referred to a portion of the trajectory.

Further analysis was carried out to discern also the effect of random GPS errors on the trajectory of the kinematic path. The magnitude of the effect due to random errors was identified through the computation of the deflection angle given by the angle between a V'W' line and the prolongation of the preceding U'V' line (Figure 3). Deflection angles were computed as positive or negative values depending on whether the line lied to the right (clockwise) or left (counterclockwise) of the prolongation of the preceding line.

**Figure 3.** Scheme of the deflection angles referred to a portion of the trajectory.

#### *3.3. Video Camera Observations*

A video monitoring system was installed in the western part of Serapo beach (41◦12 41.81 N–13◦33 29.66 E) as shown in Figure 4. Two cameras provided a total view of the bay with 1920 ×1080 pixel resolution from an elevation of about 11 m above Mean Sea Level (MSL) and 100 m from the coastline. The mean pixel resolution was evaluated using the methodology suggested in Holland et al. [46]. Images were collected every second from 08:00 until 16:00 local time on a local video-recorder with a 2-TB storage capability from May 2017–June 2018. Optical measurements are subject to image distortions due to inherent camera characteristics [47]. According to Stumpf et al. [48] and Mucerino et al. [29], the images were calibrated using Ground Control Points (GCPs), which were placed in the view area of the cameras and were acquired in UTM-WGS84 using DGPS. The calibration process required at least nine GCPs for each camera that covered the entire camera view; camera system calibration was performed by using 32 GCPs: 21 on the southeast side camera and 11 on west side camera. Finally, the image dataset was processed using Beachkeeper plus software [15].

Shoreline detection from the image was based on the physical consideration that the color contrast between beach and water is sufficient, lighting is strong enough, and the number of pixels in the water and beach groups is sufficient. Following Almar et al. [36] and other less recent authors [49–52], in order to reduce errors due to sea level variations, the shoreline was detected using the swash signature on timex average images, using images averaged over short periods (30 s), which significantly improved the accuracy of shoreline determination. The shoreline thus identified on the oblique image was converted to real-world coordinates using ortho-rectification techniques. Finally, shorelines obtained from corresponding images from the two cameras were managed to form a whole continuous shoreline.

**Figure 4.** Camera system located on Serapo beach: (**a**) camera records, southeast side of the bay; (**b**) camera detections, west side of Serapo shoreline; (**c**) a detail of the camera installation.

#### *3.4. Validation Method*

The comparison between video and DGPS coastline has been performed on the whole beach, measuring the error as the deviation from the continuous DGPS line computed, depicted in green, along the normal to the DGPS itself for a sample of about 500 points. These points were chosen keeping constant the mutual distance (about 2.7 m) along the track.

In this paper, an error model for video camera measurements is proposed. The displacement offset with respect to the normal direction to the coastline can be computed by an exact trigonometric relation. In Figure 5, let *ε<sup>N</sup>* represent the error in the normal direction to the coastline, *ε<sup>L</sup>* the longitudinal error along the line of sight from the camera to the object, and *ε<sup>T</sup>* the transverse error, perpendicular to the line of sight. The angle *α* is the angle between the coastline and the line of sight. If *α* has been observed, then:

$$
\varepsilon\_N = \varepsilon\_L \sin \alpha + \varepsilon\_T \cos \alpha \tag{1}
$$

**Figure 5.** Scheme of the longitudinal and transverse error projections along the normal direction to the coastline.

The error model was based on two error components, the longitudinal *ε<sup>L</sup>* and the transverse *ε<sup>T</sup>* error, which were projected in the direction of the coastline's normal. The longitudinal and transverse errors were computed from the following equation:

$$
\varepsilon\_L = \varepsilon\_T = aD - be^{\text{cs}} \tag{2}
$$

where *D* is the distance of the shoreline from the video camera, *s* is the intertidal beach slope, *e* is Napier's constant, and *a*, *b*, *c* are three constants to be determined experimentally. In Equation (2), the correction *becs* may be applied to the distance-dependent error *aD* as the intertidal beach slope increases, which is in agreement with the inverse proportionality of the normal error with the beach slope [36].

Substituting Equation (2) into Equation 1 and if the longitudinal and transverse errors are grouped, this yields:

$$
\varepsilon\_N = (aD - b\epsilon^{cs})\left(\sin\kappa + \cos\kappa\right) \tag{3}
$$

#### **4. Results**

#### *4.1. Wave Climate Conditions*

Significant wave height Hs, mean wave period Tm, and mean wave direction Dm, covering the years 1989–2014 (with some data gaps) are reported in Figure 6a–c, respectively. Maximum values of Hs for the ten highest wave storms recorded in the observation period are reported in Table 2, together with the mean period Tm, peak period Tp, and mean direction Dm associated with the storm peak. The highest Hs value was observed in December 1999, in agreement with Piscopia [53]. The selection of the ten highest storms showed that the directions of the storm waves were confined between 238◦ N and 272◦ N. This result is consistent with the one obtained with the selection of all the storms (Hs > 2 m) highlighted in red in Figure 6, showing that the highest waves were associated with southwestern, western, and northwestern directions [37]. With regard to astronomical sea level variations, the study area experiences a typical semi-diurnal tide, with a mean tidal range of 0.35 m, following the official Italian tide archives (http://www.mareografico.it, last access: 30 April 2019). However, main sea level variations due to meteorological surges can reach values up to1m[54].

**Figure 6.** Offshore wave parameters (Hs in (**a**), Tm in (**b**), and Dm in (**c**)) recorded by the Ponza buoy during the years 1989–2014 with three-hourly time steps. The red color highlights the wave parameters associated with the highest waves (Hs > 2 m).

**Table 2.** Summary of the ten highest wave storm events recorded by the Ponza buoy in the period 1989–2014.


#### *4.2. Comparison between Different Accuracies of GPS Solutions*

The results of the comparison of the single point and kinematic instantaneous GPS solutions with the continuous one are shown below.

As shown in Figure 7, the errors of the horizontal positions obtained in single-point mode reached a maximum value of 2.45 m, with a mean horizontal error of 0.93 m (Figure 7a). It should be noted that this average value was greater than the average error of the kinematic instantaneous solutions (Figure 7b) equal to 0.49 m. The standard deviation of the single point position errors was, however, comparable to that of the kinematic case, with a value of 0.28 m. Figure 7a also denotes a clear systematic error in terms of direction with prevailing azimuths towards northeast.

**Figure 7.** Single-point (**a**) and kinematic instantaneous (**b**) horizontal position error.

From the polar histogram of the azimuths shown in Figure 8, it is noted that about 85% of the single point error directions fell in the first quadrant from 0–90 degrees, with about 65% of the directions concentrated between the azimuth of 15 degrees and the azimuth of 75 degrees. This systematic error was reduced for the directions of the kinematic instantaneous errors; in fact, 85% of the errors were distributed in a range of greater amplitude, which extended from 0–240 degrees.

**Figure 8.** Histogram of the single-point and kinematic instantaneous error position azimuths.

With regard to the deflection angle analysis reported in the Methods, we computed its standard deviation, which was much lower for the kinematic continuous trajectory with respect to the single point and kinematic instantaneous ones, as reported in Figure 9. This circumstance confirms the capability of the kinematic continuous option for smoothing GPS data due to the use of Kalman filter.

**Figure 9.** Histogram of the deflection angles.

#### *4.3. Comparison between DGPS and Video Camera Coastline*

In Figure 10, the comparison between the video camera coastline (in green) and the DGPS one (in red) is reported for four adjacent west-east beach sectors. The error distribution is not only dependent on the camera distance, as explained in the following.

**Figure 10.** Video camera coastline (in green) compared with the DGPS coastline (in red) along Serapo beach. (basemap c 2018 Google product). In order to guarantee a better representation of the two coastlines, the beach has been divided into four sections, shown in the panels **a**–**d**, starting from west (panel **a**) to east (panel **d**).

In Figure 11a, we report in the chromatic scale the deviation between the two coastlines along the beach. Moreover, in Figure 11b, we report the error magnitude, neglecting the sign, as a function of the shoreline length. The following results can be synthesized.

First, the error exhibited a different amount between points located at a distance lower or higher than approximately one third of the beach extent, as shown from Figure 11b. In particular, for the eastern camera images, from the shoreline length of 200 m to a distance of 500 m, the error was around 1 m, while for the more distant points, the error was higher, with a maximum value of 5 m.

**Figure 11.** Observed deviations of the video coastline in the normal direction to the DGPS line for about 500 sample points. Panel a shows the deviations in the chromatic scale, in panel b the deviations are reported as function of the shoreline length.

Secondly, the error showed a lower standard deviation for the closer points; this is evident from the quite uniform color in Figure 11a if compared with the color variation along the second part of the coastline. The summary of the error statistics is reported in Table 3.


**Table 3.** Summary of the error statistics.

The inspection of the deviations between the remotely-measured coastline and the direct DGPS ones showed a lower error amount for points located at a distance closer than one third of the beach extent. This result has been obtained both in terms of maximum error (less than 2.0 m) and for the error standard deviation (less than 0.3 m). On the contrary, the more distant points exhibited a maximum error up to 5.0 m and an error standard deviation up to 1.3 m.

#### **5. Discussion**

In this paper, a shoreline detection technique from a low-cost video monitoring station was used, applying an automatic extraction technique that undertakes a water/beach distinction from color band ratios and a shoreline slope recognition module.

Prior to the validation of the video-derived shoreline with the direct survey, the differences between the single point and the unfiltered kinematic GPS have been examined. The results showed that the unfiltered kinematic GPS positions exhibited a better performance of the single-point GPS with respect to the systematic error. Nevertheless, the latter presented a comparable standard deviation and a reduced cost both in terms of instruments and operational speed.

Compared to a direct DGPS survey, the video-derived coastline acquisition was less time consuming and more cost effective. The possibility to acquire a beach topography with a high temporal frequency can potentially highlight coastal processes during the winter season, when a direct survey is difficult due to harsh weather conditions.

To perform the error analysis, as stated previously (Equation (3)), it is useful to break the normal error into a distance-dependent part *ε <sup>N</sup>* and a slope-dependent one *ε N*.

The distance-dependent part can be estimated as:

$$
\varepsilon\_N' = aD\left(\sin\alpha + \cos\alpha\right) \tag{4}
$$

where the constant *a* is computed experimentally from the field data; the values were 0.0125 and 0.0079 for the west side camera (from 0–200 m) and the southeast side camera (from 200–1340 m), respectively. Thus, the distance-dependent error is represented in Figure 12, with positive (negative) values indicating seaward (landward) offsets of the shoreline. The light blue curve represents normal errors that increase with the distance compared to the observed values (red curve). Figure 12 shows results obtained when the model is not corrected for the beach slope.

**Figure 12.** Distance-dependent part of the normal error.

In performing normal error estimation, it should be assumed that there are possible corrections due to the increasing intertidal beach slope. The model adopted to estimate these corrections is:

$$
\epsilon\_N^{\prime\prime} = b \epsilon^{\rm cs} \left( \sin \mathfrak{a} + \cos \mathfrak{a} \right) \tag{5}
$$

where the constants *b* and *c* are computed experimentally from the field data and the slope *s* is expressed in percent units. The constants *b*, *c* and the slope correction *ε <sup>N</sup>* have been calculated only for the stretch of coast between Transect 1 and Transect 5 for which the slope values have been obtained by interpolation. The values were 0.005 and 1/3 for the constants *b* and *c*, respectively. Figure 13 shows the relationship between the slope and the correction due to the beach slope.

Substituting the appropriate values into Equation (3), yields:

$$
\varepsilon\_N = (0.0079D - 0.005e^{\frac{\theta}{2}}) \left(\sin \kappa + \cos \kappa\right) \tag{6}
$$

The result depends on the combination of three components including camera distance *D*, beach slope *s*, and the angle *α* between the coastline and the line of sight. The error model for video camera measurements has been tested to the extent between Transect 1 and Transect 5 due to the field data not enabling calculation of the beach slope and determining the error offset for the entire study area. Figure 14 shows the estimated normal errors determined by the final computed Equation (6). The particular combination produced results corresponding to those observed with the exception of the value at the shoreline length of 600 m.

**Figure 13.** Slope correction due to the increasing intertidal beach slope between Transects 1 and 5.

**Figure 14.** Estimated normal errors.

#### **6. Conclusions**

Currently, risk assessment plays a fundamental role in preventing irreversible erosion processes, as well as flooding damages. In this context, video monitoring represents a low-cost instrument to store a large amount of data with a high temporal frequency irrespective of weather conditions.

In this paper, we analyzed the limitations of the video-based coastline in terms of distance from the point of view and influence of beach topography. A shoreline detection technique from a low-cost video monitoring station was validated against DGPS-derived shoreline.

We introduced an error model, based on the transverse and longitudinal error along the line of sight, and computed a correction factor that can be applied to the distance-dependent error.

Wave storms alter the amounts of wave energy approaching a shore and can change the beach cross-shore profile, influencing the applicability of video monitoring acquisition systems.

Several types of coastal information can be provided from images. In Di Luccio et al. [16], we analyzed wave run-up; in the near future, we will use data for rip-current detection, which we already analyzed with direct and remote observation tools [14].

Day-light limitations of the video-derived data can be overcome by thermal cameras, which can operate also during night hours in order to monitor the beach during severe events. This goal will be the continuation of the present research, together with the challenge of modeling new algorithms to lower the deviation from direct measurements.

**Author Contributions:** Conceptualization G.P., G.B., U.R., D.D.L., L.M. Supervision G.B. Methodology G.P., G.B., U.R., D.D.L., L.M. Formal Analysis G.P., U.R., D.D.L. Validation G.P., G.B., U.R., D.D.L. Visualization G.P. and D.D.L. Software R.M. Data Curation R.M. Investigation, L.M. Writing—Original Draft Preparation G.P., U.R. and D.D.L. Writing—Review and Editing G.P., U.R. and D.D.L.

*J. Mar. Sci. Eng.* **2019**, *7*, 137

**Funding:** The research activity has been supported by the Parthenope University of Naples with a grant within the call "Support for Individual Research for the 2015–17 Period". The above support is gratefully acknowledged.

**Acknowledgments:** The authors are grateful to the forecast service of the University of Napoli "*Parthenope*" (http: //meteo.uniparthenope.it) for the HPC facilities and to the beach club "Nave di Serapo" (www.navediserapo.it) for having hosted the video monitoring equipment.

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

#### **References**


c 2019 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/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Journal of Marine Science and Engineering* Editorial Office E-mail: jmse@mdpi.com www.mdpi.com/journal/jmse

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18