Quantifying Estuarine Hydrometeorological Coastal Hazards Using a Combined Field Observation and Modeling Approach

: Coastal development and its associated site management have rapidly expanded to estuarine environments while continuing to increase worldwide. With the growth of coastal management projects, ﬁeld observations are required to understand how anthropogenic construction, coastal defense, environmental restoration, and conservation areas will react to the typical, extreme, and long-term conditions at the proposed sites. To address these unknowns, we present a multi-faceted coastal risk assessment of a unique, recently nourished estuarine beach near the mouth of the Delaware Bay Estuary by merging rapid-response remote sensing platforms, hydrodynamic models, and publically available monitoring datasets. Speciﬁcally, hydrometeorological events from 2015 to 2019 were the focus of peak-over-threshold statistics, event type deﬁnition, and clustered event interval determination. The 95th percentile thresholds were determined to be the following: 0.84 m for the signiﬁcant wave height, 13.5 m/s for the 10-m elevation wind speed, and 0.4 m for the total water level residuals. Tropical and extra-tropical cyclones, light gales, or cold and stationary fronts proved to be the meteorological causes of the sediment mobility, inducing the hydrodynamics at the site. Using these event types and exceedance instances, clustered meteorological events were deﬁned as having an interval greater than twelve hours but less than ﬁve days to be considered clustered. Clustered events were observed to cause greater volumetric change than individual events, and are currently underrepresented in coastal risk planning and response in the region. Coastal monitoring ﬁeld measurements should consider clustered events when conducting post-hazardous or erosional event response surveys. This work highlights the importance of clustered hydrometeorological events causing estuarine coastal risk, and how to quantify these effects through combined ﬁeld observations and modeling approaches.


Introduction
The coastal management of estuarine environments faces a multitude of risk types beyond the traditional open-coast events.While both environments on the East Coast of the United States are subject to tropical (TC) and extra-tropical (XC) cyclones, frontal passages, and relative sea-level rise, the shallow nature of estuaries like the Delaware Bay Estuary and the Chesapeake Bay means that they regularly experience high freshwater river discharge events, seiching, swell events, and localized increased total water levels during sustained wind events [1].These event types are capable of producing hazardous coastal conditions due to inundation and erosion, particularly in estuarine environments where the dune height, berm width, and beach slope are mild and provide less protection to the surrounding infrastructure [2].Low-energy beaches such as sandy or marsh-fronted coastlines within estuaries are influenced by large sedimentation events such as TC, XC, or river discharges, rather than the seasonal wave climates which are typically responsible for sandy, open-coast morphologies [3].Because many estuarine coastlines are microtidal, a minor increase or decrease in the total water level (TWL) exposes additional areas to wave and wind energy, resulting in potentially hazardous sediment transport and morphological change [4].
Due to the continued rapid development in all coastal areas, including estuaries, increased amounts of infrastructure and human populations are more frequently exposed to these hazards, potentially causing economic and life-threatening damages [5,6].Vulnerability assessments, early-warning systems (EWS), and in situ monitoring for threshold exceedances are historical and current methods for the protection of human life and infrastructure [7].These systems combine in situ, modeled, and probabilistic data for sites of interest, such that managers and decision-makers can assess site-specific vulnerabilities to inundation or erosion, or create a post-event response if it is needed from support services [6].The advantage of these systems over historical approaches such as storm indices is the possibility of forecasted and real-time hazard identifications when up-to-date in situ data such as topobathymetric surveys are coupled with real-time water level observations and forecasted wave heights [8].The disadvantage is the need for up-to-date topobathymetric datasets and in situ wave height observations in estuaries where shallow depths prevent vessel-based bathymetric surveys, and the high turbidity reduces topobathymetric lidar depth penetration [9].Numerical and probabilistic models are able to fill gaps in the hydrodynamic and morphological data; however, these methods do not account for rapid elevation changes due to coastal defense projects or the response of estuarine coastlines to multiple hydrometeorological events in rapid succession, which are often called clusters [10].
Clustered events in coastal hazard management are traditionally not considered and are rarely measured in situ, given the unpredictability of multiple events impacting the same area and the logistical resources needed to observe such events [11].Compounding the difficulty of measuring clustered event effects are the multiple definitions in the existing literature on how to define the interval length and what events are significant to be considered [12,13].The more obvious events are TCs, XCs, and mesoscale tropical-like cyclones such as 'Medicanes' [14,15].However, regions such as the Mid-Atlantic of the United States and estuarine environments are inundated or eroded under less severe or clustered lowerenergy events.Examples of these hazards include stationary fronts coupled with an XC or several days of a slow-moving cold front bringing strong winds, fetch-generated waves, and localized surge [16].Looking at the effects of clustered meteorological events on coastal risk, the peak-over-threshold approach (POT) using the 95th or higher percentile value for TWL or wave height is standard practice for the identification of hazardous events [17].Cluster intervals may be selected from 3 h up to several days or weeks depending on the site and the purpose of the interval.From a morphological perspective, events could be clustered if the beach profile has not recovered to pre-event conditions, which may take years depending on the location.From an estuarine coastal hazard perspective, the profile recovery approach is unrealistic, as sites may never recover to a pre-storm morphology [18].Therefore, a combined approach of event identification using POT methods coupled with in situ morphological and hydrodynamic measurements is presented here as a realistic method for risk management decision-making and inclusion into existing EWS.
In order to achieve this goal, a case study of Broadkill Beach, Delaware, located within the DBE on the East Coast of the United States, is presented as an example of the combination of pre-and post-storm topobathymetric surveys, in situ high-resolution hydrodynamic time series, forecasted Wave Watch III (WW3) wave heights, and publically available monitoring datasets of the predicted and observed water level, wind velocity, and cyclone presence to aid estuarine risk assessment and post-event hazard response.The main hypotheses addressed are:

•
Wind speeds and TWL residuals created by cyclonic and non-cyclonic meteorological events cause potentially hazardous inundation and erosion events; • Temporally grouped or clustered meteorological events cause increased inundation and erosion hazards compared to clustered events that occur individually.
Estuarine coastal hazard monitoring and risk assessments would benefit from multiple parameter thresholds rather than a singular threshold due to the codependency of sediment transport processes in these areas.This work aims to generate utilizable estuarine coastal hazard thresholds and event type characterization for rapid incorporation into existing early warning systems (EWS), current management practices, and monitoring approaches in the DBE.

Materials and Methods
Multiple in situ acoustic instruments and remote-sensing survey techniques were employed to quantify the hydrometeorological extreme event parameters within the bay portion of the DBE from 2016 to 2019.These observations were combined with publically available topobathymetric surveys, along with meteorological, water level, and climate datasets, to determine the statistical extreme-event thresholds and the prevalence of clustered meteorological events within the DBE region.Appendix A contains information relating to data collection, processing, and analyses.The following sections detail the field campaigns, publically available dataset metadata, and analytical methods.

Study Site Characteristics
The Mid-Atlantic region of the United States lies along the East Coast of the country and is comprised of the following states: New York, New Jersey, Delaware, Maryland, and Virginia.The Chesapeake Bay and Delaware Bay estuaries reside within this region and are two of the largest estuaries-in terms of surface area-within the continental U.S.These estuaries host dense infrastructure, tidal saltmarshes, freshwater wetlands, agriculture, and forested areas [19].Meteorologically, the Mid-Atlantic experiences intermittent tropical cyclones, more-frequent extra-tropical cyclones, and strong frontal passages with gale force or greater winds [20,21].
Broadkill Beach, on the western shore of the DBE, was chosen as the study site for this work.The estuary itself is critical to the U.S. East Coast shipping industry and supports over 50% of the oil and goods coming into the coast [22].It is a trumpet-shaped, drowned river valley consisting of the Delaware River at the northern extent (also referred to as the upper estuary); the lower estuary, with its fringing saltmarsh and muddy shores; and the southernmost portion, denoted as the bay, with low-energy, sandy beaches [23].Broadkill Beach is a 5 km stretch of the bay portion of the estuary within Sussex County, Delaware (Figure 1) [24].
Broadkill Beach is approximately 11 km north of the estuary mouth, and it hosts a small community of full-time and part-time human residents, migratory bird groups (Calidaris canutus), spawning horseshoe crabs (Limulus Polyphemus), reef-building polychaete worms (Sabellaria vulgaris), and the largest singular volume of nourishment sediment placed on an estuarine beach on the U.S. East Coast [9,[25][26][27][28].The beach faces 55 • cross-shore and 315 • -145 • alongshore, with winds typically occurring from the northwest to the south, with northern to southeastern winds generating the longest fetches over the DBE.This area is microtidal, with a mean range of 1.24 m and a spring range of 1.42 m.While the sediments are typically medium-to-fine quartz sand, due to extreme events and consistent small nourishments over the decades, the grain size distribution is heavily influenced by anthropogenic coastal defense projects [27,29].The most recent project was a placement of nearly one million cubic meters of dredged navigation channel material, forming a 120 m wide beach platform completed by a trapezoidal dune of 2-m height and 7-m width spanning the 5 km of Broadkill Beach [26,30].
The DBE water levels have been continuously monitored by NOAA for decades, and meteorological sensors have been added in recent years.Modeling efforts have increased in recent years for the area, and now include the Physical Oceanographic Real-Time System (PORTS), also known as DBOFS, and Wave Watch III wave height forecasts.It is important to note that these outputs are intended for EWS and navigational safety, rather than scientific research, although several studies have leveraged these long-term datasets for climatic observations such as sea-level rise and flood events [1].Biannual topographic profiles are performed for all of Delaware by the Department of Natural Resources and Environmental Control (DNREC) and are made available upon request.
The DBE water levels have been continuously monitored by NOAA for decades, and meteorological sensors have been added in recent years.Modeling efforts have increased in recent years for the area, and now include the Physical Oceanographic Real-Time System (PORTS), also known as DBOFS, and Wave Watch III wave height forecasts.It is important to note that these outputs are intended for EWS and navigational safety, rather than scientific research, although several studies have leveraged these long-term datasets for climatic observations such as sea-level rise and flood events [1].Biannual topographic profiles are performed for all of Delaware by the Department of Natural Resources and Environmental Control (DNREC) and are made available upon request.

Remote-Sensing Surveys
The surveys consisted of subaqueous sonar mapping using vessel-mounted Phase-Measuring Echosounders (PMES), subaerial mapping using Unoccupied Aerial Systems (UASs), and Real-Time Kinematics Global Position Systems (RTK GPS).These surveys generated the digital surface models, georectified imagery, and GPS point measurements used in the analyses described below.The standardized volumetric change values for cyclonic events from 2012 to 2018 were leveraged from Dohner et al. [9] and were used here with hydroacoustic data to assess meteorologically induced coastal hazards.In-depth details of the survey methods can be found in Dohner et al. [9], including field collection parameters, data post-processing, and volumetric change calculations; as such, a brief synopsis is presented here with additional details in Table A1.
The topographic surveys consisted of two approaches: RTK GPS surveys from 2012 to 2015, and UAS aerial surveys from 2016 to 2018.Both methods aim to measure the preand post-cyclone coastal response at Broadkill Beach.RTK GPS surveys auto-collected points every two meters in the cross-shore direction, with profiles spaced every 20 to 30 m alongshore, depending on the survey.The UAS surveys consisted of 80% along-and cross-track overlap of the images with ground control points (GCPs) placed approximately every 100 m.Natural control points such as pilings, beach chairs, and decks were used to determine the DSM error for each survey.The GCPs and natural control points were surveyed using RTK GPS.The survey images were stitched using Structure-from-Motion techniques using Pix4D and Metashape (formally Agisoft Photoscan) software to create the surfaces.
The bathymetric surveys were completed using two PMES systems: Ping DSP 3DSS (450 kHz) for depths of less than 3 m, and Edgetech 6205 (230 and 550 kHz) for depths greater than 3 m.The Ping 3DSS system used a D-GPS (Hemisphere GPS) and internal IMU system for navigation and heave correction while the Edgetech 6205 was coupled to two Novatel RTK Antennas with a Coda F190 IMU.Periodic CTD casts were made during the surveying using a Castaway CTD system to collect the sound speed velocity profiles.Tidal records were taken from NOAA station 8,557,380 for post-processing purposes.The raw sonar data was processed in SonarWiz (Cheseapeake Technology Inc., Fall River, MA, USA) software in order to create gridded bathymetric surfaces.The topo-and bathymetric surfaces were imported into ESRI ArcGIS Pro, corrected to MSL, and combined into single surfaces for analysis.

Hydroacoustic Time Series
In situ hydroacoustic data were collected intermittently from 2016 to 2019 using a suite of acoustic sensors mounted on an aluminum frame which was deployed on the seafloor for approximately 1-2 months per deployment at depths ranging from 3.0 to 6.4 m (MLLW).The acoustic frame configuration consisted of a 600 kHz Teledyne RDI Acoustic Doppler Current Profiler (TRDI ADCP), a 2 MHz Nortek Aquadopp HR Profiler (Pulse-Coherent), a 2.25 MHz Imagenex 881 fan-beam rotary sonar (RS), and an AML Plus-X CTD with a turbidity sensor.This frame created a 2 m × 2 m footprint on the seafloor, with the TRDI ADCP-mounted 0.9 m off the seafloor-collecting burst waves and currents every 30 min in 25 cm bins.The downward-facing Aquadopp was mounted 0.5 m from the seafloor, and measured burst currents every 30 min in 3 cm bins.The Aquadopp was oriented offshore to reduce frame effects as much as possible.The raw data were copied from each instrument following recovery and were processed in propriety software from Teledyne and Nortek to output usable formats of wave (height, direction, period), current velocity, and backscatter amplitude bursts.The TRDI and Aquadopp values were concatenated into a single column for a near-complete record of the full water column.Blanking distances were included based on the instrument settings, mounting distance, and frame location.The cumulative time in the water resulted in 12 months of acoustic data spanning February-November during the 3 years of the deployments.No deployments occurred during December and January due to instrument maintenance, dangerous weather, and vessel availability.Surface sediment grabs were performed for the frame deployment and survey operations using a Young's grab sampler.The samples were processed using the pipette method owing to the estuarine samples containing organics, clays, and silts [31].Any sand fractions were dried and sieved using a Ro-Tap machine and sieves sizes of −2 to 3 Φ in half-phi steps.The sediment characteristics and statistics were determined using GRADISTAT [32].

Publically Available Data
Hourly and 6-min meteorological data were gathered from NOAA stations 8,557,380 and LSWD1.These stations are both located at the Lewes-Cape May Ferry terminal.Station 8,557,380 also provided predicted and verified water levels in the same sampling frequencies as the Met data.The data from station 8,557,380 are managed by the CO-OPs program, and are available from their websites, while data from station LSWD1 are managed by the NDBC program.
TWL above the MHHW exceedance values were found at the CO-OPs station 8,557,380 Extreme Water Levels webpage.Additional significant wave height, wave periods, and TWL exceedance values for station 8,557,380 were collected from the USACE ERDC North Atlantic Coast Comprehensive Study (NACCS) [33,34].Supplemental cyclonic characteristics from historical events which were not measured in the field campaigns were included in the analyses.
Significant wave heights from 2015 to 2019 were collected from the Wave Watch III (WW3) system [35].The grid cells were selected by matching the acoustic frame deployment locations.The WW3 and TRDI burst-average significant wave heights were time-synchronized and correlated using linear regression.This regression was used to correct the WW3 wave heights before further analysis [36,37].
Cyclonic events (TC, XC) were identified from the NWS Daily Weather Maps (DWM) from 1945 to 2019 and existing data from Dr. Daniel Leathers of the University of Delaware [38].The cyclonic events were characterized as closed isobars surrounding a low-pressure system if it fell between 70 • 00 -80 • 00 W and 30 • 00 -40 • 00 N (as previously set by Leathers et al. 2011).The cyclones were divided into types, such that tropical depression, tropical storm, hurricane, and extra-tropical were identified.Additional details and figures for this topic are found in Supplementary S1.
Biannual topographic profiles from DNREC spanning 1995 to 2019 were acquired directly from the agency for the study site.These profiles were performed with various RTK GPS systems throughout the years.These data provided long-term morphological context and feature tracking for comparison with the remotely-sensed data from 2015 to 2018.Pre-and post-nourishment bathymetric surveys (cross-shore profiles) from USACE were incorporated into the nearshore slope and normalized volumetric calculations.

Analyses 2.5.1. Spatial
Pre-and post-cyclone topobathymetric surfaces and profiles were used to determine the normalized volumetric change (∆V), dune toe and crest elevations, and beach slope (nearshore and foreshore).∆V was calculated by surface differencing sequential DSMs, determining the volume change of the surface and normalizing for the overlapping surface area of the surveys used.Dune crests were identified as the maximum elevation within the sandy beach area, while dune toes were identified as the maximum curvature before the dune crest [39].Post-storm UAS orthomosaics provided the visual identification of wrack lines and wet-to-dry sediment lines for the digitization of the HWL and shoreline.Previously published values for these parameters from 2012 to 2017 are incorporated here for further analysis alongside new values obtained from multiple cyclones in 2018 and 2019 at the study site [9].

Coastal Risk
Using the available in situ times series, the 95th percentile values (95%) were determined for the H s , T p , wind speed (U 10 ), total water level greater than MHHW (TWL MHHW ), and total water level residuals (TWL res ).TWL MHHW was determined as the magnitude of the difference between the observed water levels and the station MHHW vertical benchmark of 1.418 m above MLLW.TWL residuals are the absolute values of the difference between the predicted and observed water levels.The distributions for H s , U 10 , and TWL res were fitted with Generalized Extreme Value functions, while the wind direction (θ wind ) utilized a kernel function, and TWL MHHW used a logistic function.These distributions were sampled and paired to create the Joint Probability heat plots (JP) found in the Results section.The distribution fit values and JP look-up tables are presented in Appendix B. Extreme event values for measured coastal risk events were incorporated into the return period equations created by the NOAA station extreme event analysis and the NAACS modeling report at station 8557380.
The 95% threshold values of H s and U 10 coupled with barometric pressure (P atm ) were then used to identify instances in time when these thresholds were exceeded, with the inclusion of WW3 H s to extend and fill gaps in the in situ data.These exceedance instances were classified into three meteorological phenomena based on the following: Any instances meeting the above classifications were cross-referenced with the DWMs to ground-truth the classification scheme.The 'Storms' were TC or XC based on the methods described in Section 2.4.'Wind' events were confirmed as instances of large gradients between cyclonic and anticyclonic systems near the East Coast or frontal passages.'Sunny' events comprised front passages or distant strong cyclonic events outside the geographic zone identified for 'Cyclones'.The categories and their justifications are discussed further in the following sections.
Aeolian and bedload sediment mobility thresholds were calculated based on in situ sediment grabs, U 10 six-minute time series, and Aquadopp currents in the bin nearest to the seafloor.The Aeolian critical threshold for sediment mobility was determined using Shields' formula for quartz sand (d = 0.2 mm, ρ s = 2650 kg/m 3 ): where A = 0.1, ρ = 1.225 kg/m 3 (uncorrected for temperature and humidity), and g = 9.81 m/s 2 [40].The wind velocities (U * ) at the surface (z = 0 m) were calculated using the following equations from Zingg [41]: where z = the sensor height in meters, U z = the wind velocity at the sensor height in m/s, z 0 = the roughness length, κ = 0.4, and d = the sediment size in meters.The corresponding U 10 for U *c was found by solving for U 10 in Equation ( 2).The bedload thresholds were determined using the instantaneous Shields dimensionless parameter with a conditional wave and current shear stress equations: where ρ s is the density of the sediment, ρ is the density of the water, g is the gravitational acceleration constant, d is the mean grain size diameter, τ w is the shear stress under wave conditions, U w is the near-bed velocities due to waves (Equation ( 9)), τ c is the shear stress under current conditions, C d is the drag coefficient determined from the grain diameter and associated roughness length, and U c is the near-bed velocities due to currents measured by the Aquadopp.The values from Equation ( 4) were compared to the sediment mobility critical threshold from Whitehouse and Soulsby [42] when the sediments were noncohesive (i.e., they had a mud content <5%): Equation ( 5) includes the wave friction factor of [43] and the wave orbital velocity calculated as follows from the TRDI wave data: where h is the depth, k is the wavenumber (2π/L), H s is the significant wave height, and T p is the peak period.Run-up exceedance and depth of closure (h c ) parameters are applied methods for management planning, coastal risk assessment, and project design lifetime.These values are calculated for the available TRDI H s and depth time series, and then later compared with the dune toe elevation as a management criterion for major morphological change following meteorological events or clusters of events.The run-up heights were determined using the Stockdon et al. [44] general formulation, and the beach slopes were determined from the UAS DSMs: where R 2% is the run-up height exceeded by 2% of the swash events, β f is the beach slope over the portion of the beach where the run-up occurs (foreshore slope), H s,0 is the deepwater significant wave height, and L p,0 is the deepwater peak wavelength.This analysis utilized the available frame H s time series from the site to calculate h c using the simplified version of Birkemeier [45]: where h c is the depth of closure and H e is the significant wave height of the top 0.137% waves in a year [46].Select storm h c values were determined using the Shields parameter form of the equation by Hallermeier [47]: where U b is the near-bed orbital velocity from Equation (9).

Clustered Morphological Events
As presented in the introduction, the concept and implementation of clustered events as a coastal hazard is typically decided site by site, when it is considered at all.Because this work aims to aid coastal managers, monitors, engineers, and policy-makers, several definitions for the cluster interval were applied to the available data, and the applicability of each was weighed in the discussion.Intervals of 2, 3, 4, 5, 7, 9, and 14 days were chosen based on the existing interval definitions [48].An alternative method and interval definition for the DBE is presented as well.This interval was determined based on the 95% threshold of the six-minute TWL res from 2015 to 2019.Individual events were identified when the TWL res was greater than the 95% threshold (the start of an event) and fell below the threshold (the end of an event).The intervals between the end of one event and the start of the next event were fitted to determine the mean of the distribution.Intervals of less than 12 h were excluded, as these were considered to be the same event.
To ground-truth the event and cluster identification, 95% threshold dates were compared to the non-tidal water column backscatter (~0.3 m from the seafloor), in which the backscatter amplitude is considered to be a qualitative turbidity estimate.This was carried out by demeaning the backscatter of each bin over individual deployments using a 25-h moving window average due to the semidiurnal tides; this is referred to as the non-tidal backscatter.Then, identified events greater than the 95% threshold for each bin were time-synchronized with the backscatter residuals.The ground-truthing for the hazard event type resulted in 95% threshold exceedance events temporally corresponding with 95% backscatter residual exceedance events.Finally, a site-specific cluster interval was determined from the acoustic backscatter residual time series.Clustered events were defined as events that exceeded the 95% backscatter residual threshold and did not return to the mean (zero) before another event.This was an additional method for the determination of the cluster interval length based on subaqueous sediment mobility in the form of acoustic backscatter, rather than inferred mobility from hydrodynamic or meteorological parameters.It is important to note that this is one method to confirm whether the identified threshold events were capable of potentially causing morphological change.There are limitations to this method, which are discussed in later sections.

Morphological Change
Critical subaerial morphological features such as the elevations of the dune crest, dune toe, and berm were mapped and extracted from the northern, central, and southern portions of the study site from 1995 to 2019 (Figure 2).An overall loss of feature elevations occurred from 1995 to 2008 in the northern and central regions, while the southern portion was not mapped until 2008.Following 2008, the elevations appear erratic and inconsistently switch been high erosion and accretion.These changes were due to TC Sandy (2012) and the immediate XC, which caused an overtopping of the existing dune crests, the failure of the dune, and widespread flooding at the site.
Following this sequence of events, an emergency sediment placement occurred, creating the accretional event in 2012.This cycle continued in 2013 with 27 XCs, TD Dorian, and TS Andrea, followed by sediment nourishment.The year 2014 was similar, with 43 XCs, TS Bertha, and TCs Arthur and Cristobal, followed by sediment nourishment at the site.The year 2015 gave way to 31 XCs, TS Ana, and TC Joaquin, with the largest site nourishment to date starting in the fall of this year and completing in the spring of 2016.After 2016, the berm elevations at all of the locations eroded consistently through 2019, while the dune crest elevations remained stable.The dune toe elevations oscillated between erosion and accretion due to the installment of cross-shore and alongshore sediment fencings in 2017.While the minor accretion from 2017 to 2019 may suggest that the coastal defense structures are stable, the overall trend of Figure 2 for the dune toe and berm features is erosional, which is supported by the volumetric measurements [9] and obvious field observations of dune toe scarping and walk-over damage.Whilst they are unlikely to experience overtop-ping due to the nourishment dune crest elevation, the failure of the dune during collision events of mild wave heights and water elevation is possible.The project design suggested maintenance nourishments every 3-5 years; at the time of this work, no nourishments have taken place at the site [26,49].
J. Mar.Sci.Eng.2022, 10, x FOR PEER REVIEW 10 of 26 obvious field observations of dune toe scarping and walk-over damage.Whilst they are unlikely to experience overtopping due to the nourishment dune crest elevation, the failure of the dune during collision events of mild wave heights and water elevation is possible.The project design suggested maintenance nourishments every 3-5 years; at the time of this work, no nourishments have taken place at the site [26,49].Considering the whole profile, rather than discrete features, the normalized volumetric change due to individual and clustered storm events was calculated and compared to the peak TWLMHHW for each event or cluster (Figure 3).Note that TC Joaquin occurred in 2015 before much of the large nourishment was placed, and XC Jonas occurred before the project completion in 2016.All of the other events in Figure 3 occurred post-nourishment.TC Sandy (2012) is considered a benchmark morphological event for the site due to the major erosion and flooding, making it a reasonable comparison for the following events.Several interesting observations come from Figure 3: the large values of volumetric change, XC Jonas being comparable to TC Sandy, clustered XCs causing measurable erosion, and TCs Jose and Hermine being potentially accretional events.Considering the whole profile, rather than discrete features, the normalized volumetric change due to individual and clustered storm events was calculated and compared to the peak TWL MHHW for each event or cluster (Figure 3).Note that TC Joaquin occurred in 2015 before much of the large nourishment was placed, and XC Jonas occurred before the project completion in 2016.All of the other events in Figure 3 occurred post-nourishment.TC Sandy (2012) is considered a benchmark morphological event for the site due to the major erosion and flooding, making it a reasonable comparison for the following events.Several interesting observations come from Figure 3: the large values of volumetric change, XC Jonas being comparable to TC Sandy, clustered XCs causing measurable erosion, and TCs Jose and Hermine being potentially accretional events.
Regarding the magnitude of the values, several factors may contribute to the scale of change observed.Several events occurred following a large, unconsolidated placement of sediment.This placement was exposed to record-setting water levels and wave heights from multiple events, as well as strong, clustered XCs.While each pre-and post-storm survey was temporally near to each cyclonic event, other hydrometeorological events occurred around these storms and potentially contributed to the volumetric change.Therefore, the magnitude of the change is reasonable, and the relative comparison between events in Figure 3 holds.XC Jonas created record-setting wave heights, as calculated by the correct WW3 of 1.99 m on top of the 1.38 m of the TWL MHHW .This combination contributed to the large volume loss, which was comparable to TC Sandy.As will be presented later in this section, XC events occurring in rapid succession-such as the four March 2018 XCs-are capable of causing a great deal of erosion.This is most likely due to the estuarine hydrodynamic conditions not recovering to a pre-storm state and allowing sediment to settle out of the water column.This is acutely important in estuarine settings due to the mixture of organic particles, fine sediments, cohesive sediments, and the nonlinear hydraulics of a semi-enclosed basin [50].An unexpected result of Figure 3 was the accretional changes caused by TCs Hermine (2016) and Jose (2017).Both occurred after the nourishment, and the pre-storm surveys were taken post-nourishment and pre-storm, thereby avoiding bias due to the sediment placement.Both TCs occurred geographically near the site to be captured in the DWM datasets but did not meet the criteria of the storm events category based on the wave, wind, and barometric pressure.Neither event made landfall; however, both events moved slowly and lingered offshore for several days.Ultimately, these two events were classified as sunny-day swell (both) and wind events (TC Hermine only), rather than storm events.TC Jose showed strong non-tidal backscatter residuals in the near-bed acoustic time series, suggesting high suspended and bedload particle concentrations.The acoustic frame was not deployed for TC Hermine and therefore was not present in the non-tidal backscatter dataset.This shows a curious statistic for coastal defense projects; simply treating all of the storm events as erosional in design life calculations is unrealistic, as these events are potentially accretional for estuarine environments.Regarding the magnitude of the values, several factors may contribute to the scale of change observed.Several events occurred following a large, unconsolidated placement of sediment.This placement was exposed to record-setting water levels and wave heights from multiple events, as well as strong, clustered XCs.While each pre-and post-storm survey was temporally near to each cyclonic event, other hydrometeorological events occurred around these storms and potentially contributed to the volumetric change.Therefore, the magnitude of the change is reasonable, and the relative comparison between events in Figure 3 holds.XC Jonas created record-setting wave heights, as calculated by the correct WW3 of 1.99 m on top of the 1.38m of the TWLMHHW.This combination contributed to the large volume loss, which was comparable to TC Sandy.As will be presented later in this section, XC events occurring in rapid succession-such as the four March 2018 XCs-are capable of causing a great deal of erosion.This is most likely due to the estuarine hydrodynamic conditions not recovering to a pre-storm state and allowing sediment to settle out of the water column.This is acutely important in estuarine settings due to the mixture of organic particles, fine sediments, cohesive sediments, and the nonlinear hydraulics of a semi-enclosed basin [50].An unexpected result of Figure 3 was the accretional changes caused by TCs Hermine (2016) and Jose (2017).Both occurred after the nourishment, and the pre-storm surveys were taken post-nourishment and pre-storm, thereby avoiding bias due to the sediment placement.Both TCs occurred geographically near the site to be captured in the DWM datasets but did not meet the criteria of the storm events category based on the wave, wind, and barometric pressure.Neither event made landfall; however, both events moved slowly and lingered offshore for several days.Ulti- Shifting from a site-wide, morphological viewpoint of coastal events, Table 1 presents the hydrometeorological parameters based on the sediment grain size, depth of closure, wave run-up heights, and sediment mobility thresholds for the study site.The in situ d 50 enabled the calculation of the critical mobility thresholds based on the wind speed and near-bed current speed (presented as shear stress).These values allowed for the additional confirmation of sediment mobility during the identified event types, showing that the transport occurred subaerially and subaqueously during non-storm events.Using the depth-of-closure and run-up elevation values, a combined conceptual model of the coastal hazards was created based on the 2019 elevation profiles (Figure 4).The instantaneous values for the four March 2018 XCs were the following: Riley (3.1 m), Quinn (2.2 m), Skylar (2.5 m), and Toby (3.3 m).The point here is that it impresses the importance of coupling up-to-date topobathymetric values with values derived from in situ observations to assess the risk within the DBE due to hydrometeorological events.The berm and dune toe across all of the profiles were subject to continuous mobility events-either storm, sunny-day swell, or wind-supporting the results of Figure 4.
for the additional confirmation of sediment mobility during the identified event types, showing that the transport occurred subaerially and subaqueously during non-storm events.Using the depth-of-closure and run-up elevation values, a combined conceptual model of the coastal hazards was created based on the 2019 elevation profiles (Figure 4).The instantaneous values for the four March 2018 XCs were the following: Riley (3.1 m), Quinn (2.2 m), Skylar (2.5 m), and Toby (3.3 m).The point here is that it impresses the importance of coupling up-to-date topobathymetric values with values derived from in situ observations to assess the risk within the DBE due to hydrometeorological events.The berm and dune toe across all of the profiles were subject to continuous mobility events-either storm, sunny-day swell, or wind-supporting the results of Figure 4.

Exceedance Thresholds and Event Identification
Having established the effects of individual and clustered hydrometeorological events on the study site morphology, exceedance thresholds for the identification of these events were created from the in situ datasets, and are presented in Table 2.

Hazard Identification Exceedance Thresholds and Event Identification
Having established the effects of individual and clustered hydrometeorological events on the study site morphology, exceedance thresholds for the identification of these events were created from the in situ datasets, and are presented in Table 2.These values were either derived from in situ monitoring systems in the DBE or forecast using numerical models, and are thus readily available for ingestion and utilization by existing EWS.Using these values, management groups are able to decide when and where to assess potential damage to coastal areas based on the topobathymetric values at sites of concern, as conditions exceeding these thresholds are likely to have caused considerable morphological change in a short period.
Combining the exceedance thresholds of Table 2 and the event type classifications, the relationships between the individual hazards are considered in Figure 5. Immediately, the large gap in events from 200 to 270 • is evident.However, this is unsurprising, as this sector contains land rather than open water.Also expected is the fact that storm and wind events generate higher wave heights, and occur from the NW to E directions.Sunny-day swell events are the only type generated from the E to S directions.The directions are typical for offshore swells in the open ocean of this region, particularly from offshore TCs and XCs.Interestingly, all of the event types coincide with the exceedance of the TWL res threshold, which is more closely associated with non-tidal backscatter events than any other parameter.The backscatter exceedance events were found to be statistically different from typical backscatter values (α = 0.05, p < 0.001), and were therefore considered a strong indicator of sediment flux in the area.
the relationships between the individual hazards are considered in Figure 5. Immediately, the large gap in events from 200 to 270° is evident.However, this is unsurprising, as this sector contains land rather than open water.Also expected is the fact that storm and wind events generate higher wave heights, and occur from the NW to E directions.Sunny-day swell events are the only type generated from the E to S directions.The directions are typical for offshore swells in the open ocean of this region, particularly from offshore TCs and XCs.Interestingly, all of the event types coincide with the exceedance of the TWLres threshold, which is more closely associated with non-tidal backscatter events than any other parameter.The backscatter exceedance events were found to be statistically different from typical backscatter values (α = 0.05, p < 0.001), and were therefore considered a strong indicator of sediment flux in the area.  2 for the values).
Being mostly meteorological in origin, the event types were expected to be seasonally dependent.Figure 6 shows instantaneous occurrences matching the event type criteria  2 for the values).
Being mostly meteorological in origin, the event types were expected to be seasonally dependent.Figure 6 shows instantaneous occurrences matching the event type criteria and the associated date.These event occurrences were compared to instances of positive non-tidal backscatter residuals to confirm that turbidity was present in the water column during an event.Of the identified events that overlapped with frame deployments, 26 out of 26 storm, 54 out of 55 sunny, and 7 out of 8 wind events were associated with positive backscatter residuals (100%, 98%, and 87.5%, respectively).When comparing the events to the backscatter 95% thresholds, 20 out of 26 storm, 27 out of 55 sunny, and 5 out of 8 wind events exceeded the 95% threshold (76.9%, 49.1%, and 62.5%, respectively).This equated to 97.8% of the events being corroborated with positive backscatter residuals, and 58.4% of the events being confirmed to exceed the 95% threshold.While not all of the identified events corresponded with higher turbidity residuals, the aforementioned percentages confirm that the event type criteria captured sedimentation events related to hydrometeorological forces.It is interesting to note that the threshold was exceeded more often by sunny events than the combined total of storm and wind events.Sunny events were also more frequently responsible for positive backscatter residuals, but fewer sunny events rose above the 95% threshold.This suggests that sunny events more frequently create turbidity events, while storm and wind events are more likely to be extreme turbidity events.
gion.From a coastal hazard perspective, sunny-day events would correspond to wavy conditions in the bay with no gale or small craft advisory and no obvious approaching storms.While they are unlikely to cause catastrophic profile change or extreme hazards to boating individually, these events last several days and may lead to unexpected flooding and profile change, which may exacerbate the effects from previous or approaching storm and wind events.The three event types corresponded with TWL res events over the study period, and support the importance of these types as coastal hazards in an estuarine environment.While they are not a hazardous event category, the TWL res exceedance occurrences were far more frequent than any event category.The occurrences centered on the winter and early spring months of December-April, the peak time for XCs and strong wind events in the region.Sunny-day swell events occurred as the most frequent event type and were associated with the following synoptic phenomena based on the DWMs: stationary fronts, cold fronts, and gradients between high-to low-pressure systems in the surrounding region.From a coastal hazard perspective, sunny-day events would correspond to wavy conditions in the bay with no gale or small craft advisory and no obvious approaching storms.While they are unlikely to cause catastrophic profile change or extreme hazards to boating individually, these events last several days and may lead to unexpected flooding and profile change, which may exacerbate the effects from previous or approaching storm and wind events.

Return Periods
Pivoting from defining coastal hazard types to risk assessment, return period values for the DBE are presented in the context of the event types determined in the previous subsection (Figure 7).The values for a selection of historical storm events and hydrometeorological events were compared using probabilistic modeling of H s return periods and extreme-value rank methods for TWL MHHW .As expected, XC Jonas (2016) was a rare event for this area, as was TC Sandy and the March 1962 XC.As presented earlier in the Results section, TWL MHHW is a reliable indicator of sediment flux, but it is not the only hazard that is important to major coastal risk.The joint return periods of TWL MHHW and H s highlight an unnamed event in January 2017 that was greater than the existing 100-year return period.Unfortunately, there were no field or modeling data available for the wave heights during TC Sandy to compare with these XCs.Interestingly, multiple wind events corresponded with high return periods similar to extreme storm events.Conversely, sunny-day swell events were consistently lower than wind and storm events but did cluster with typical storm events.As a reminder, wind events must have barometric pressures greater than 1015 mb, and sunny-day swell events have wind below the 95% threshold.riod.Unfortunately, there were no field or modeling data available for the wave heights during TC Sandy to compare with these XCs.Interestingly, multiple wind events corresponded with high return periods similar to extreme storm events.Conversely, sunnyday swell events were consistently lower than wind and storm events but did cluster with typical storm events.As a reminder, wind events must have barometric pressures greater than 1015mb, and sunny-day swell events have wind below the 95% threshold.

JPA
Seeing the connectivity of coastal hazards in the estuarine environment and the potential increased risk when occurring in tandem, the individual hazard probabilities were fitted and compared to determine the JP of such coupled hazards (Figure 8).These heat maps visualize the most likely conditions experienced by the coastlines within the DBE and are useful for project design, maintenance intervals, and design life expectancies.The most probable wave heights are above 0.5m when TWLMHHW is below the MHHW line.Translating this to the real-world beach profile, larger waves were most probable from the dune toe to the berm, suggesting that the collision regimes for this area occur frequently (Figure 4).This result points to the foreshore being frequently exposed to higher wave

JPA
Seeing the connectivity of coastal hazards in the estuarine environment and the potential increased risk when occurring in tandem, the individual hazard probabilities were fitted and compared to determine the JP of such coupled hazards (Figure 8).These heat maps visualize the most likely conditions experienced by the coastlines within the DBE and are useful for project design, maintenance intervals, and design life expectancies.The most probable wave heights are above 0.5 m when TWL MHHW is below the MHHW line.Translating this to the real-world beach profile, larger waves were most probable from the dune toe to the berm, suggesting that the collision regimes for this area occur frequently (Figure 4).This result points to the foreshore being frequently exposed to higher wave energies, such that scarping and dune collapse are possible.Interestingly, positive TWL MHHW are most probable from the E-S directions, while the largest TWL res are associated with the N-E directions.This is a combination of the DBE orientation relative to the predominant wind patterns of the region, as well as the average frontal and storm passage tracks.The wind direction for U 10 , TWL res , and H s showed higher probabilities from N-E quadrants than other directions, while TWL MHHW was most probable from the E-SE directions.This disparity represents the climatic or more-common wind directions associated with the longest fetch in the DBE (N-E) and the more extreme events generating winds and higher water levels from the E-SE directions.
the N-E directions.This is a combination of the DBE orientation relative to the predomi-nant wind patterns of the region, as well as the average frontal and storm passage tracks.The wind direction for U10, TWLres, and Hs showed higher probabilities from N-E quadrants than other directions, while TWLMHHW was most probable from the E-SE directions.This disparity represents the climatic or more-common wind directions associated with the longest fetch in the DBE (N-E) and the more extreme events generating winds and higher water levels from the E-SE directions.

Event Clusters
Leveraging the DWM storm identifications with the in situ data and exceedance thresholds, the intervals between each event type, as well as all of the exceedance events, were determined (Figure 9).The intervals were chosen based on the existing literature and the in situ data with a minimum interval of 12 h to prevent the identification of a lingering event as a cluster.

Event Clusters
Leveraging the DWM storm identifications with the in situ data and exceedance thresholds, the intervals between each event type, as well as all of the exceedance events, were determined (Figure 9).The intervals were chosen based on the existing literature and the in situ data with a minimum interval of 12 h to prevent the identification of a lingering event as a cluster.
From the varying intervals and event types, there are several observations.First, the long-term DWM dataset shows the clustered storms of TC or XC origin that occurred in this area since 1851.Thus, at a minimum, clustered cyclonic events are an additional category of coastal hazards for the DBE, and need to be incorporated in management plans for improved safety.Second, events purely above the U 10 threshold are a yearly occurrence for this area and are clustered as well.This is an important consideration for aeolian transport because the U 10 is nearly double the critical threshold for the initiation of sediment motion.Particularly for the sandy beaches of the DBE, aeolian transport most likely constitutes consistent profile change until a larger event passes through the region.Third, storms and sunny-day events cluster in larger numbers than wind events.This seems unexpected at first, but considering the sequence of synoptic-scale meteorological events of the area, fewer stationary fronts, offshore storms, and high-to low-pressure gradients occur than simply strong wind events.Finally, Figure 9 demonstrates that multiple hydrometeorological events cluster, and that the interval can vary based on the analysis.
Using the in situ backscatter amplitude residuals and the TWL res , the sediment mobility-based cluster interval in the DBE region was determined to be 5 days.This interval length arose from an event initially exceeding the 95% backscatter residual and remaining above the mean before a successive event (Figure 10).An example for one bin is shown for 0.3 m off of the seafloor (Figure 10).The three example clusters are a combination of Sunny to Wind to Sunny (blue), Storm to Wind to Sunny to Storm (cyan), and Sunny to Wind to Sunny (orange).These examples support Figure 9's results in that Wind and Sunny-day events occurred frequently in this interval length or shorter alongside fewer, but still evident, storm events.These are expected pairings because frontal passages, pressure gradients, and offshore swells are more likely to occur in succession due to mesoscale meteorological forcings with less-frequent cyclonic events occurring back to back.This mixture of hazardous events highlights the impacts of hydrometeorological hazards before, after, and unconnected with cyclonic systems.From the varying intervals and event types, there are several observations.First, the long-term DWM dataset shows the clustered storms of TC or XC origin that occurred in this area since 1851.Thus, at a minimum, clustered cyclonic events are an additional category of coastal hazards for the DBE, and need to be incorporated in management plans for improved safety.Second, events purely above the U10 threshold are a yearly occurrence for this area and are clustered as well.This is an important consideration for aeolian transport because the U10 is nearly double the critical threshold for the initiation of sediment motion.Particularly for the sandy beaches of the DBE, aeolian transport most likely constitutes consistent profile change until a larger event passes through the region.Third, storms and sunny-day events cluster in larger numbers than wind events.This seems unexpected at first, but considering the sequence of synoptic-scale meteorological events of the area, fewer stationary fronts, offshore storms, and high-to low-pressure gradients occur than simply strong wind events.Finally, Figure 9 demonstrates that multiple hydrometeorological events cluster, and that the interval can vary based on the analysis.
Using the in situ backscatter amplitude residuals and the TWLres, the sediment mobility-based cluster interval in the DBE region was determined to be 5 days.This interval length arose from an event initially exceeding the 95% backscatter residual and remaining above the mean before a successive event (Figure 10).An example for one bin is shown for 0.3 m off of the seafloor (Figure 10).The three example clusters are a combination of Sunny to Wind to Sunny (blue), Storm to Wind to Sunny to Storm (cyan), and Sunny to Wind to Sunny (orange).These examples support Figure 9's results in that Wind and Sunny-day events occurred frequently in this interval length or shorter alongside fewer, but still evident, storm events.These are expected pairings because frontal passages, pressure gradients, and offshore swells are more likely to occur in succession due to mesoscale meteorological forcings with less-frequent cyclonic events occurring back to back.This

Discussion
The in situ extreme event observations, pre-and post-event mapping, and regional event climatological data were combined to determine the extreme-event thresholds and storm cluster effects within the DBE.This low-energy, sandy estuarine beach is located 11

Discussion
The in situ extreme event observations, pre-and post-event mapping, and regional event climatological data were combined to determine the extreme-event thresholds and storm cluster effects within the DBE.This low-energy, sandy estuarine beach is located 11 km north of the mouth and is nourished typically following extreme erosion.Field observations of estuarine beach dynamics and extreme event response are poorly represented in the existing literature relative to those of open, sandy coasts, due in part to the difficulty in observing these systems at high spatiotemporal resolutions [3].XC morphological impacts, as well as clustered meteorological events, are rarely captured in field observations, particularly along the entire subaerial-to-subaqueous profile [3].Nor are these clustered events considered in risk assessment or project planning [6].Compounding these issues, the existing literature provides multiple definitions for cluster intervals depending upon the discipline presenting the interval [48,[51][52][53][54]. Therefore, the following discussion leverages unique morphological and hydrodynamic datasets alongside the publically available data for the DBE to present a method grounded in field data to establish extreme-event thresholds, event types, and event cluster intervals.These results are immediately applicable to existing EWS and future coastal management plans.

Extreme Events Types and Thresholds
The extreme event statistical analysis of decades of in situ data covering hazardous events is a useful measure for managers, project designers, and risk assessments as inputs to vulnerability and probability-of-occurrence calculations [6].The applicability of these statistics depends upon the assumption that point data from a single gauge or instrument station are valid for the surrounding area.This assumption is commonly applied to opencoast sites so long as the shoreline orientation, sediment type, and regional hydrodynamics are consistent.These assumptions for open-coast extreme events do not apply to large estuarine sites such as the DBE, as the previously stated morphological and hydrodynamic parameters change drastically within a small geographic region [16].This is seen in other sites around the world such as Venice, Italy; Rethymno, Greece; San Francisco Bay, U.S.A; Catalonia, Spain; and New South Wales, Australia [48,[55][56][57][58][59].These works chose the event parameters which are most applicable to their sites, such as TWL in low-lying Venice or H s related to massive erosion on sandy beaches.No such analysis has been conducted for the interior of the DBE owing to a lack of long-term hydrodynamic monitoring within the estuary.However, the studies stated previously (e.g., NACCS and NOAA) determined the return period relationships for Delaware's Atlantic coast, or modeled probabilistic TC parameters to determine such extreme event return periods.
This work presented multiple in situ observations from 2012 to 2019 to determine extreme-event thresholds using the 'peaks-over-threshold' method, in which the threshold was established at 95% of the available parameter distribution (Table 2) [48].Using these thresholds and publically available meteorological and water level data, additional hydrometeorological events became apparent as potential estuarine coastal hazards, rather than the typical thresholds focused on H s [52][53][54].Figure 2 showed the rapid morphological change associated with synoptic meteorological events, but these singular events alone did not account for the volumetric change measured in Figure 3. Instead, additional hydrometeorological events contributed to the elevated TWL res , H s , and ∆V, such as stationary fronts, regional gradients between high-and low-pressure systems, troughs, and offshore swell events.These event types were ground-truthed as sediment mobility events through the non-tidal acoustic backscatter residuals functioning as a qualitative turbidity estimate (Figure 10, Supplementary S2).
The incorporation of these additional event types into morphologic and flood risk assessments adds another layer of warning for the existing EWS in the DBE.While the return period of many non-storm events is lower those of than historic TC or XC events, Figure 7 demonstrated the nontrivial impact of XC and windy events with return periods greater than 20 years.Certainly, wind events leading up to or following TC and XC may have contributed to the large volumetric changes observed in Figure 3, as the pre-and poststorm surveys were weather-condition dependent due to UAS and small-craft restrictions.Particularly with the wind events, the critical threshold of motion for the subaerial sediment was half of the threshold speed used as an event identifier.Therefore, the wind velocities consistently exceeded the critical threshold of motion before the extreme-event threshold, and suggest that aeolian transport is a major process for the sandy beaches near the estuary mouth.Directionally, these higher wind speeds are from the NW-NE and flow alongshore uninterrupted.This transport was seen in situ as the net accretion of the nourishment dune and raised dune crest elevation (Figure 2).These changes were supported by overlaying UAS orthomosaics with the accretional areas that geographically matched the cross-shore and dune toe sand fencing, as well as the planted dune crest vegetation [9].
The subaqueous sandy sediment supply to this site historically arrives from flood tidal currents around the Cape Henlopen spit complex [60,61].Due to the littoral transport nodal point at the center of this site, the net littoral transport for the southern portion is southeast towards the mouth, while the net littoral transport for the northern portion flows northwest [25].Because of this, the typical sediment conditions for the nearshore region of the site are erosional.Exceptions occur when offshore extreme events create the potential to force additional coastal sandy sediment into the bay, as was seen during TCs Hermine and Jose (Figure 3).
Recognizing the logistical difficulties of monitoring a large estuary such as the DBE, the statistical relationships between physical hazard parameters were determined using the JP method between the dominant environmental parameters (Figure 8).Rather than highlighting extreme events as the threshold analysis did, these values serve as resources for the planning, risk assessment, and management of the region.The usefulness of the JP method lies in the values incorporating quiescent periods, extreme events, and nonlinear relationships between the parameters which are not readily captured by deterministic models such as WW3 and DBOFS [55].Statistical modeling methods such as those used in the NAACS report make assumptions regarding the storm shapes (i.e., the rate of wave height change, duration, etc.), paths, and associated surge, and typically do not consider non-cyclonic events in morphological risk assessments.

Clustered Events
While debate remains regarding the definition of storm clusters and little field data exist to evaluate the coastal response to storm clusters, the effects of coupled or grouped events differ from singular events [51,56].TCs, XCs, and frontal passages are capable of generating a high degree of wave chop, and of transporting suspended and bedload sediments [56].Here, several cluster intervals were considered based on existing literature which used meteorological data, morphological data, or a combination of the two as guidelines to decide a site-specific interval.We recommend using a 5-day interval for the DBE based upon a meteorological and sediment mobility approach, as shown in Figure 10.The Mid-Atlantic region is frequently exposed to rapidly moving XCs with lingering high-wind events due to the resulting regional pressure gradient [62,63].Occurrences of strong winds and frontal passages mobilize and suspend the mud and mixed sediments, particularly during neap tidal cycles when the water depths are reduced and wave energy travels deeper into the water column.This leaves DBE conditions well above non-storm levels and lends credence to the inclusion of events such as gales, squalls, and frontal passages in clusters from a geomorphological perspective.The Mid-Atlantic, the Mediterranean, the Iberian Peninsula, and Southeastern Australia experience several types of energetic meteorological events which contribute to nearshore and aeolian fluxes [9,17].These events would likely be lower in magnitude than the tidal or storm fluxes in the overall system but may significantly contribute to small-scale biological, chemical, and geological processes [64,65].The importance of, at minimum, clustered TC and XC was shown in Figure 3 for the four March 2018 XCs and 2012, with TC Sandy followed swiftly by an XC.This work demonstrates the importance of including clustered hazardous events in EWS within estuarine environments, as these shallow, low-energy systems are more morphologically responsive to multiple hydrometeorological processes than open-coast environments.

Limitations and Future Research
It is important to note the limitations of these findings and to present measures to address the limitations in future work.The normalized volumetric changes at the site are biased due to the particularly large amount of sediment placed from 2015 to 2016.This is why the fitted line is a qualitative representation, rather than a useful quantitative predictive equation.The morphology of this particular beach is vastly different from the remainder of the estuarine beaches; therefore, the applicability of the volumetric change relationship to other DBE beaches is questionable.In the same respect, the dune toe elevation mapped in July 2019 serves as a threshold for extreme events only at Broadkill Beach (Figure 4).The remaining have little or no existing dune system, and the lower estuary transitions to fringing marsh followed by developed waterfronts in the upper estuary near the Delaware River [23].A consideration in the use of TWL res as an in situ threshold is the lack of run-up and swash zone processes incorporated into that parameter.Calculated run-ups of 0.5 m or more for each March 2018 ETCs further eroded beach features like the nourishment dune toe, as observed from the orthomosaics and volumetric change maps [9].Given these limitations, the 95% exceedance values for H s (Table 2) have greater applications for the region, but require a permanent monitoring buoy within the DBE-where one currently does not exist-or more ground-truthing of WW3 forecasts, as they consistently under-predict extreme events [37].The relationship between volumetric change and H s is indirect and nonlinear, as other parameters such as the tidal phase, storm duration, river discharge, and wind velocities are significant contributors within estuarine systems [65,66].Therefore, the TWL relative to the dune toe elevation is a practical and useful method for the assessment of extreme events across the estuarine beaches if the semiannual RTK GPS profiles taken along the Delaware shoreline by DNREC are continually incorporated into EWS.Coupling this threshold with the measured beach slopes (e.g., RTK GPS profiles) provides inundation vulnerabilities at each profile location [6].Utilizing the continuous RTK GPS monitoring profiles, the inundation vulnerability could be determined for each profile along the western DBE.Thus, the hazard findings presented here serve as useful elements for coastal risk assessment for the lower DBE using high-resolution field observations.
The subaqueous sediments within the DBE vary spatial, and are typically mixed and cohesive rather than noncohesive [67].This points to potential cohesive effects on the bedload and suspended transport near the bed, while the upper water column showed high turbidity, as is expected for a shallow water body during an energetic event [68].There is a potential for the suspended sediment flux to be larger than the bedload due to the amount of mud present at the site as well as the high levels of backscatter seen throughout the water column [23,42].Cohesive sediments are expected to be agitated and suspended during energetic events due to the small grain size, and they remain in suspension until flocculation and lower energies persist [69].This supports the use of the non-tidal backscatter as confirmation of sedimentation events, rather than simply H s or critical thresholds of motion for sediments.
Given these limitations, the statistical methods coupled with the modeled and in situ data address the major concerns presented above.Future research aimed at collecting hydrodynamic data from December to January, as well as inside the depth of closure, would close the data gaps of the in situ data used here.Ideally, field collection at the different coastal types in the DBE of alongshore and cross-shore subaerial-to-subaqueous sediment flux would provide high-resolution spatiotemporal data for the validation of the thresholds, event types, and cluster relationships presented in this case study.These types of field experiments would aid in the modeling of sea-level rise impacts on eroding coastlines, coastal defense projects such as the nourishment dune at this site, and how sediment supply is affected by changing sea levels and changing coastal hazards [70][71][72][73].

Figure 1 .Figure 1 .
Figure 1.Study site characteristics: (a) outline of the United States with a red box outlining the DBE; (b) zoomed-in map of the DBE with the study site marked with a black star.NOAA wave buoys 44,009 and 44,001 are marked with black circles, while NOAA station 8557380/LSWD1 is marked with a black triangle.(c) UAS imagery (2018) from the center of the study site looking southwardFigure 1. Study site characteristics: (a) outline of the United States with a red box outlining the DBE; (b) zoomed-in map of the DBE with the study site marked with a black star.NOAA wave buoys 44,009 and 44,001 are marked with black circles, while NOAA station 8557380/LSWD1 is marked with a black triangle.(c) UAS imagery (2018) from the center of the study site looking southward towards the bay entrance; (d) frame deployments (red circles) with dates and the survey area for the study site (green hashed box).Inlaid image of the acoustic frame with an upward-facing 600 kHz Teledyne RDI ADCP (yellow box), downward-facing 2 MHz Aquadopp HR Profiler (blue box), AML CTD with a turbidity sensor (red box), and an Imagenex 881 Rotary Sonar (green box).

Figure 2 .
Figure 2. Elevations relative to MHHW for the northern, central, and southern regions of the study site.The features tracked include the dune crest, dune toe, and berm.The missing points indicate that the feature was not detectable in the data and most likely not present along the profile.Each line represents a profile feature, in which the colors are the same location and the line types are the same feature.The profile surveys were performed by DNREC and are made publically available upon request.

Figure 2 .
Figure 2. Elevations relative to MHHW for the northern, central, and southern regions of the study site.The features tracked include the dune crest, dune toe, and berm.The missing points indicate that the feature was not detectable in the data and most likely not present along the profile.Each line represents a profile feature, in which the colors are the same location and the line types are the same feature.The profile surveys were performed by DNREC and are made publically available upon request.

Figure 3 .
Figure 3. Normalized volumetric change (ΔV) at the site from 2012 to 2018 using various topobathymetric sensors and platforms.Individual TC and XC events are shown, as well as the combined values for the four XC occurring during March 2018.The fitted line shows the negative log relationship between the ΔV and the TWLMHHW for the events.

Figure 3 .
Figure 3. Normalized volumetric change (∆V) at the site from 2012 to 2018 using various topobathymetric sensors and platforms.Individual TC and XC events are shown, as well as the combined values for the four XC occurring during March 2018.The fitted line shows the negative log relationship between the ∆V and the TWL MHHW for the events.

Figure 4 .
Figure 4. Topobathymetric profiles were completed in July of 2019 for the northern, central, and southern portions of the study site.The yellow boxes indicate areas of submergence during normal and hydrometeorological conditions, relative to MSL.

Figure 4 .
Figure 4. Topobathymetric profiles were completed in July of 2019 for the northern, central, and southern portions of the study site.The yellow boxes indicate areas of submergence during normal and hydrometeorological conditions, relative to MSL.

Figure 5 .
Figure 5.The above scatter plots show the results of the conditional event identification process based on 95% thresholds (Hs, TWLres, and U10).These types are compared to the barometric pressure (top row) and wind direction (right column) for 2015 to 2019.The grey shaded region denotes compass directions overland.The black dashed lines are the threshold for the respective parameters (see Table2for the values).

Figure 5 .
Figure 5.The above scatter plots show the results of the conditional event identification process based on 95% thresholds (H s , TWL res , and U 10 ).These types are compared to the barometric pressure (top row) and wind direction (right column) for 2015 to 2019.The grey shaded region denotes compass directions overland.The black dashed lines are the threshold for the respective parameters (see Table2for the values).

Figure 6 .
Figure 6.Occurrences of threshold exceedances based on the storm, wind, sunny, and TWL res 95% criteria.The bar widths are one month.

Figure 7 .
Figure 7. (Left) Return period plot of the total water level above MHHW for the Lewes gauge station for a selection of historical storms compared to the statistical return period values created by USACE ERDC during the NACCS project.See Appendix B for a list of storm names, dates, and values used in this figure.(Right) ERDC return period values for Hs and TWLMHHW plotted with corrected WW3 wave heights, NOAA tide station water levels, and U10 wind speeds from 2015 to 2019.

Figure 7 .
Figure 7. (Left) Return period plot of the total water level above MHHW for the Lewes gauge station for a selection of historical storms compared to the statistical return period values created by USACE ERDC during the NACCS project.See Appendix B for a list of storm names, dates, and values used in this figure.(Right) ERDC return period values for H s and TWL MHHW plotted with corrected WW3 wave heights, NOAA tide station water levels, and U 10 wind speeds from 2015 to 2019.

Figure 8 .
Figure 8. JP heat maps for the combinations of the following parameters: U10 velocity, Hs, TWLres, and TWLMHHW.The distributions consisted of field data from 2015 to 2019 for the water level and wind velocity, while the wave heights originated from in situ acoustic deployments.The heat maps surrounded in red boxes correspond to the JP tables in Appendix B.

Figure 8 .
Figure 8. JP heat maps for the combinations of the following parameters: U 10 velocity, H s , TWL res , and TWL MHHW .The distributions consisted of field data from 2015 to 2019 for the water level and wind velocity, while the wave heights originated from in situ acoustic deployments.The heat maps surrounded in red boxes correspond to the JP tables in Appendix B.

Figure 9 .
Figure 9. Bar plot showing the intervals between various datasets such as the DWM (1851-2019), the categorical event types, the TWLres 95% threshold, and the event dates from 2015 to 2019.The bar widths are one day, and the intervals are cumulative occurrences greater than 12 h.Not all of the event types had clusters for every interval amount.The interval amounts were chosen based on the existing literature and the in situ datasets.The values for this figure are found in Appendix B.

Figure 9 .
Figure 9. Bar plot showing the intervals between various datasets such as the DWM (1851-2019), the categorical event types, the TWL res 95% threshold, and the event dates from 2015 to 2019.The bar widths are one day, and the intervals are cumulative occurrences greater than 12 h.Not all of the event types had clusters for every interval amount.The interval amounts were chosen based on the existing literature and the in situ datasets.The values for this figure are found in Appendix B. J. Mar.Sci.Eng.2022, 10, x FOR PEER REVIEW 18 of 26

Figure 10 .
Figure 10.Residual backscatter values from all of the available deployments (top).The values here are an example from 0.3 m off of the seabed, while the analysis included all of the bins of the water column.The cluster examples include Sunny-Wind-Sunny (blue), Storm-Wind-Sunny-Storm (cyan), and Sunny-Wind-Sunny (orange).

Figure 10 .
Figure 10.Residual backscatter values from all of the available deployments (top).The values here are an example from 0.3 m off of the seabed, while the analysis included all of the bins of the water column.The cluster examples include Sunny-Wind-Sunny (blue), Storm-Wind-Sunny-Storm (cyan), and Sunny-Wind-Sunny (orange).

Table 1 .
Parameters determined from the acoustic and meteorological time series.

Table 2 .
The 95th percentile values used as the exceedance thresholds to determine coastal hazard event occurrences from 2015 to 2019 and the morphological exceedance elevations.

Table A1 .
[31]data for the field campaigns and publically available data used in the analyses.NDBC Station LWSD1 is the same location as CO-OPs station 8557380; however, the data were only available from the NDBC website during June 2016 to March 2017. 2.3https://tidesandcurrents.noaa.gov/est/est_station.shtml?stnid=8557380 (29 November 2021).