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

: 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. Speciﬁcally, (1) the entire domain of the National Parks Service o ﬀ shore 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 identiﬁed, and (3) newly obtained data were compared to that from the pre-storm survey to determined changes that could be attributed to speciﬁc storms such as Sandy and Joaquin. Capturing event speciﬁc 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. to document storm-related changes in and marine habitats on multiple scales: From large scale, side-scan sonar maps and interpretation of acoustic bottom types, to characterize as as


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 2 of 30 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.

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][5][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][8][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][9][10][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][16][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.

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.

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.

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

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  Island National Seashore mapping project: from EdgeTech raw data to SonarWiz to shapefiles and Google Earth .kml data products.

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

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. 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 .75 meters, 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. 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 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.

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.

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 2 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. 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 2 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

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.

Acoustic Mapping Products and Findings?
In total, over 65 km 2 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 m 2 /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).     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. 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. 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. 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 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. 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. 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.   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.

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.

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 14 and 15) and up to 80 m in the southern portion ( Figure 16, Figure 17, Figure 18 and Figure  19). The shift in mud-coarse sediment class boundary however, was measured to be up to 180 m in the southern part ( Figure 17).  (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).

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

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.    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 meters) were recognized in bathymetry cross section profiles ( Figure 22 and Figure 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 centimeters are recognized (Figure 24). The largest vertical change was recognized in October 2015 profile J2 (Figure 24), which suggests axial incision of about 65 centimeters. 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.        (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).

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. 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 ( Figure 22 and Figure 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 (Figure 20 and Figure 21).

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. 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 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. inherent uncertainty in the study, this benthic morphodynamic study is a valuable resource to inform management of publicly-valued and historically under-studied habitats. 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 ( Figure 11 and Figure 12) and in the seabed morphology revealed in the bathymetry (Figure 21 and Figure 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 shortterm program for immediate storm-response surveys during the hurricane season at Assateaguewhile 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. 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.