MODELING THE MORPHOLOGICAL EVOLUTION OF STABILIZED DUNE SYSTEM DURING EXTREME STORM EVENTS

The key role of beach dunes in protecting coastal developments from damages caused by seasonal storm and hurricanes (e.g. Hurricane Sandy 2012), prompted series of studies devoted to assess their stability during these conditions. A 2-D numerical study for a “100-year storm” event (1% probability of annual exceedance) in Rhode Island indicate extreme dune erosion can occur when these structural barriers become submerged during coastal flooding events, making them ephemeral solutions only (Schambach et al., 2018). In contrast, areas with dense vegetation experienced less morphological changes due to reduced wave energy. Model predictions however remain uncertain due to the lack of comparative studies with field data. The study provides a performance evaluation assessment of four modeling scenarios used to assess the stability of a vegetated dune system in barrier beaches during storm events. These are combination of two wave models (phase averaging vs. phase resolving) with two approaches to modeling the affect of vegetation on sediment transport (bed friction formulations vs. wave damping formulations). Best estimates of the post dune profiles were obtained when using the phase averaged wave model with vegetation described by bed friction coefficients. This scenario was able to predict accurate changes of the along shore crest height elevations, making it a good assessment tool for determining the vulnerability of coastal communities in Rhode Island. In a second chapter, we address the potential use of Geotextile Sand-filled Containers (GSCs) to reinforce and stabilize a dune system during storm events. A new classification of the damage states associated with these “soft-structures” used for coastal protection are identified and described. Modeling approaches to identify each damage level during a storm were validated using post survey measurements taken from a storm which caused damage to a new project site where reinforced dune system with GSCs was constructed. Two hydraulic stability formulas for GSCs were used to determine the critical conditions for the instability of the GSC structure. Results indicate lower confidence when using Hudson’s (1956) based stability equations when compared with Recio and Oumeraci (2008) semi-empirical formulations. Overall good prediction of the damage levels was obtained when compared with measurements, suggesting a potential tool for predicting the stability of the structure.


Abstract
Sea level rise and storm intensification lead to re-evaluate inundation assessments along the North Atlantic US shoreline. A particular effort is devoted to assess the coastal community risk to the "100-year storm" event (1% probability of annual exceedance) in Rhode Island (RI), US, using a chain of state of the art storm surge, wave propagation, and coastal erosion 2-D models. As part of this effort, a recent study indicated that RI beach barriers covered with dense vegetation would experience less erosion and morphological changes than similar beach barriers barely vegetated, as expected from theoretical considerations and past studies  The present study is cast in this context and provides a comparative analysis of selected hydromorphodynamic approaches for predicting dune erosion and morphological changes of the dune system, in dissipative beaches and barrier islands, for various storm regimes as described by  using the 2-D model XBeach. In particular, we investigate the applicability, performance and limitations of the selected hydro-morphodynamic models regarding their ability to assess the vegetation's effect on sediment transport during extreme storm events. We compare results of simulation for an historical storm Surfbeat ) models are investigated. Various performance capabilities and limitations are observed when each model was assessed during collision, overwash and inundation regime. Unresolved depth profile in XBeach lead to significant reduction in the depth averaged velocities and resulted in underestimation of the sediment transport rate when using the vegetation module in overwash and inundation regime. Less energy dissipation occurred when using the bed friction parameterization approach and therefore better predictions of the morphological shape of the dune profile. Comparison between the two wave models indicated incompatibility when using the Van Thiel -Van Rijin 2009 sediment transport formulation with the Non-hydrostatic approach to estimate dune erosion rates in the overwash and inundation regime due to the dominancy of the short waves orbital velocities. Reasonable predictions of the crest elevations and morphological profile of the dune were obtained when using the Surbeat wave model for all regimes.

Introduction
Extreme stages of dune erosion have been common occurrences along the US North Atlantic coast exposing coastal communities to higher inundation risk than expected based on historical storms (e.g. , 2006. The anticipated increase frequency of extreme events due to climate change has incited coastal management agencies to take protective actions. Planting vegetation along dunes crest combined with protective measures restricting their usage is a traditional approach to stabilize dunes (e.g. Woodhouse et al. 1978;Knutson, 1977).
The vegetation indeed helps to maintain the dune. From an aeolian perspective, dunes increase roughness length (Charnock, 1955), increasing the loss of energy by friction and turbulence, resulting in a lower flow velocity, ultimately trapping the sand into the dune system. When the vegetation is submerged, similar hydrodynamics processes to the aerodynamics one affect the flow promoting local sedimentation. In addition, trapped-sediments are protected from re-suspension due to the local deficit of energy in the flow, reinforcing the protective effect of the vegetation on the dune system (e.g. Gacia and Duarte, 2001;Manca et al., 2012).
The functional role of dune vegetation as a management tool to preserve the dune morphology and protect inland area from inundation has been abundantly documented, beginning with Cowles's ecological lesson (1899) relating dune morphology and "plant society. The idea of using natural and nature-based features (NNBF) between ocean and land such as dunes in particular acting as buffers to protect coastal communities however became recently increasingly popular in both the scientific and the coastal management community (e.g. Smith et al., 2016) leading to collaborative works between both stake holders (e.g., BEACH SAMP, RI). Some authors have attempted to assess the value of specific coastal ecosystems such as vegetated dunes, salt marshes and mangrove habitats in an Ecosystem-Based Management (EBM) approach (Aburto et al., 2012). Barbier et al.
(2008) demonstrate an increase value in "coastal protection service", proportional to the vegetated habitat area.
In parallel, a need to provide numerical simulations for extreme storm scenarios with adequate accuracy to address the coastal communities risk, has emerged. Although modeling storm surge and wave propagation can be done with relative accuracy, modeling dune erosion accurately is more complex.
In this work we investigate the adequacy of four selected hydro-morphodynamic approaches for predicting dune erosion and morphological changes of the dune system in dissipative beaches and barrier islands, for various storm regimes as described by  using the 2-D model XBeach . In particular, we investigate the applicability, performance and limitations of the selected hydro-morphodynamic models regarding their ability to assess the vegetation's effect on sediment transport during extreme storm events. We compare results of simulation for the historical storm Sandy using either a vegetation module based on Mendez and Losada's formulation (2004), or a friction parameterization based on a Manning formulation. In addition, the sensitivity to the selected wave model, the Non-hydrostatic (Smit et al., 2010 and Surf-beat   A background summary is provided in Section 2. The methodology is presented in Section 3 with a presentation of the selected site and of the simulated scenarios, including a brief overview of the XBeach model. Data and model set up are presented in Section 4.
Results are provided in Section 5. Discussion and conclusion are presented in Section 6 and 7 respectively.

Key Conceptual Wor k
The use of vegetation to damp wave energy has been investigated for many decades theoretically and empirically. Price et al. (1969) proposed the deployment of artificial seaweeds in the surf-zone to limit coastal erosion in England, based on theoretical development as well as tank experiments. Dalrymple et al. (1984) modified Radder's wave propagation parabolic model (1979) to include a term, which allows for the dissipation of wave energy associated to wave propagation across vegetation (Booij, 1981). The theory was restricted to damping induced by vegetation for non-breaking regular waves and horizontal bottom. Mendez and Losada (2004) expanded Dalrymple et al.'s (1984) theory to take into account the bottom variations, the randomness of the waves and the dissipation due to breaking in the surf zone. Although various formulations of the dissipation terms have been proposed, the standard formulation is based on a hypothetical array of rigid cylinders, to mimic stems, creating an additional drag force on the flow. Typically the plant-induced force acting on the fluid is expressed in terms of a Morison-type equation (Morison et al., 1950) neglecting swaying motion and inertia forces, lumping the uncertainty on the flexibility of and the relative configuration of the stems in an empirical drag coefficient. Consequently, the drag coefficient must be calibrated for specific vegetation stem types and patterns in either field studies or laboratory experiment. In later studies, this formulation has been modified to include the plant motion, coupling flow and vegetation motions (e.g. Stratigaki et al., 2011;Mendez et al., 1999;Mazda et al., 2013).

Implementation of Vegetation Modules in Numerical Modules
Semi-empirical formulations of energy dissipation associated with vegetation cover based on Mendez and Losada (2004) have been implemented into hydrodynamics models such as SWAN (Suzuki et al., 2012), STWAVE (Smith et al., 2016), NHWAVE (Ma et al., 2013) and into the hydro-morphodynamic model XBeach. This formulation provides an alternative method to assess the energy dissipation associated with vegetation to the standard friction approach in which the energy is dissipated through a friction coefficient, such as a Manning coefficient, reflecting the bed roughness associate to the specific land cover (Arcement and Schneider, 1989).
This later standard parametric method, using a spatially variable Manning coefficient with value associated to the specific land-cover, has been applied and discussed in many applications using several hydrodynamics models (e.g, in STWAVE, Wamsley et al., 2009Wamsley et al., , 2010, and more recently using the hydro-morphodynamic model XBeach . However, while the use of the Manning coefficient has been extensively validated for river discharges and

Methodology
In this work we investigate the adequacy of both (1) the friction parameterization The topography before and after Sandy as provided by USACE is used to assess the validity of the two vegetation modules.
The choice of the test site is based on two criteria, (1) the availability of pre-and post storm data of any historical storm over-washing the dune system and (2), the similarity of the vegetation to the RI shoreline, with the objective of adopting a reliable protocol and eventually apply the method in RI (where no data for storms in overwash regime are available). The model is first briefly described followed by a presentation of the selected scenarios. Data are described in the next section.

Numerical Model
XBeach ("eXtreme Beach behavior") is a 2-D coupled hydrodynamic and morphodynamic model that dynamically simulates the coastal response during timevarying storm conditions, and consequently the storm-induced changes in bed level . The model includes a choice of hydrodynamic models and a sediment transport model combined to a morphology change model. The model is designed to simulate processes occurring during the four erosion regimes defined by , swash, collision, overwash and inundation. The modeled coastal response includes dune erosion, breaching and accretion. In this work we use and compare results using two hydrodynamic modules, the Surfbeat and the Non-hydrostatic models. For both hydrodynamics approaches we use the two vegetation parameterizations. Each module of the XBeach model is briefly described hereafter.

Surfbeat Mode (phase averaging)
The "Surfbeat" model simulates mean currents in combination with a wave action conservation equation. The characteristic feature of this model is the modeling of the long wave component of swash motions, the surf-beats resulting from the forcing by wave groups of infragravity waves and causing slow oscillations of the mean super-elevated water level (period of 20 to 250 sec) directly related to foreshore erosion (Longuet-Higgins and Stewart, 1964;Schaffer et al., 1993). Surf-beats are indeed particularly dominant in the surf zone in storm conditions, since the incident-wave frequency band (0:05 < f < 0:18 Hz) is saturated, while the infra-gravity frequency band is not, in particular on dissipative beaches . While the full directional distribution of wave action is maintained, the frequency spectrum is limited to a single frequency, the peak frequency. The propagation of wind waves is based on the conservation of wave action equation, coupled to a roller energy balance equation  in the surf zone through the wave dissipation term. Radiation stress tensors acting as forcing terms for depth-and period-averaged mean current equations are  . Indeed on dissipative beaches the short waves are mostly dissipated by the time they reach the shoreline, and only infragravity waves reach the shoreline.

Non-hydrostatic Mode (phase resolving)
In Surfbeat mode, parametric functions are used to determine the nonlinear evolution of wave fields, whereas in non-hydrostatic mode these are fully resolved. Instantaneous wave velocities and water level variations are captured within the nonlinear shallow water equations by including the non-hydrostatic pressure force. In this way the dependency on the wave action equation becomes redundant. The non-hydrostatic model is particularly applicable for reflective beaches, where short waves have more significant impact on beach erosion (Smit et al., 2010(Smit et al., , 2013(Smit et al., , 2014. The numerical scheme used is as introduced by Stelling and Zejlima (2003), and further developed by Zejlima and Stelling (2008) which gives competitive predictions in shallow waters compared to higher order Boussinesq models (e.g. Chen et al. 2000), while only maintaining a single layer (Smit et al., 2010(Smit et al., , 2013(Smit et al., , 2014. By using an edge based finite method, and a momentum conservation scheme, the possibility to model wave breaking without using a separate model was shown. In XBeach, this was further developed by Smit et al. (2010) to include a limited version of a scheme introduced by MacCormack (1969) to allow for more accurate predictions of shock waves and breaking.

Vegetation Module in XBeach
The energy dissipation through drag and turbulence when waves are propagating across vegetation is modeled as a dissipation term, ! , in the wave action equation = ! ] is solved according to Mendez and Losada's (2004), with E, the energy , ! the group velocity , the onshore coordinate and , the gravitational acceleration.
with ! the time-averaged rate of energy dissipation per unit horizontal area induced by the vegetation, ! a depth-averaged drag coefficient, is the water density, the local wave number, !"# the root mean square wave height, the number of vegetation stem per unit of horizontal area, ! the plant area per unit height of each vegetation stem normal to the velocity, the wave angular frequency, the relative vegetation height (%) relative to the local water depth ℎ.
The model can include layered vegetation (e.g., to represent heterogeneous vegetation, such as mangroves; Suzuki et al. 2011). In that case the dissipation term is simply the sum of several dissipation terms, each associated with specific layer. The model can also include the swaying motion of the vegetation, which requires a recalibration of the drag coefficient (e.g. Maza et al., 2013).

Bed Friction Parameterization
The standard approach to account for the effect of vegetation on the flow velocity is using a spatially variable friction coefficient function of the bed roughness associated with the specific land cover (Kothyari et al., 1997;Wamsley et al., 2009Wamsley et al., , 2010. The Manning coefficient (n) can be expressed as (e.g., Wamsley et al., 2009): where !" is a bed friction coefficient, and h is the water depth. The value of the Manning coefficient is assigned according to land cover following Wamsley's e al.
(2009), as used in  and shown in

Sediment Transport Module
The sediment transport model used in this study is the default model in XBeach depth-averaged advection-diffusion equation . The sediment transport is controlled by the sediment concentration in the water column relative to an equilibrium concentration (Soulsby et al., 1997). While wave propagation equations are based on linear theory, the effect of non-linear waves on sediment transport is included by adding an arbitrarily parameterized advection velocity u a (Stokes drift) to the Eulerian velocity, based on a wave skewness S k and asymmetry A s parameters both weighted by arbitrary empirical coefficient, f s and f a (so called facua parameter) .

Model Parameters
Default model parameter values were used, except for cases where extensive sets of validation and recommendations were provided through comparative field analysis and laboratory experiments Devet et al. 2015;. In XBeach, these parameters are bed friction coefficients, wave skewness and assymetry "facua" parameter, morphodynamic accelerator factor "morfac", critical wet and dry slope avalanching values, and scaling factor "eps". Further sensitivity analyses of these parameters were also provided. XBeach is presented in detail in .
[Note that we used XBeach version XBeachX 1.23.5446M), released with the fortran and

Simulation Scenarios
At the test site we model the coastal impact of Superstorm Sandy (October 2012) using four different approaches based on the choice between the wave hydrodynamic module and the vegetation module as summarized in Table 1

Scenario 1 and 2
In both of these scenarios the Manning Roughness approach is used to model the effect of vegetation on sediment transport. This allows comparison between the Surfbeat and Non-hydrostatic approach. The models are used in 2-D and the friction is spatially variable based on the vegetation cover found at the study site. Vegetation data is

Scenario 3 and 4
In both of these scenarios the vegetation module is used to describe the effect of formulation in XBeach. The objective of these numerical experiments is to obtain a general understanding on how the sediment transport rate is affected by the vegetation module as compared to the Manning Approach. Results from these experiments and the field study are used to evaluate the performance and applicability of each scenario.

Assessment and Evaluation Methods
The ability of the model to simulate accurately the eroded volume is assessed by computing several skill parameters. A summary of the skill parameters used to assess the accuracy of the model for each scenario is presented in Table 1 that compares the relative value of the eroded volume. Table 1.3 -Summary of skill parameters used to assess the simulation results.

Statistical Skill Parameter
Conceptual assessment

SK
Relative value of the eroded volume

Site Description and Model Set-up
Hurricane Sandy (2012) was one of the most destructive storms to hit the US North Atlantic Coast to date. Estimated cost of damage was at least $50 billion as reported by the National Hurricane Center . Approaching its landfall location with a forward speed of 10 m/s and wind speed 38.5 m/s, it caused severe overwash, flooding and inundation along the coastlines of New York and New Jersey in particular . Many of the undeveloped sections of the barrier islands close to the landfall location left pampered after the incident.
In response to this event, the USGS in collaboration with National Oceanic and Atmospheric Agency (NOAA) and other state and federal agencies developed a database of the areas affected by Sandy, consisting of large scale aerial photographs, bathymetric and topographic LiDAR survey measurements, and land cover surveys. This data captured many of the morphological changes that occurred during the storm and is utilized for this study.   .

Pre and Post Survey Measurements
Topography for the pre and post event was extracted from the NOAA DEM database

Vegetation and Land Cover
The spatial variability of land cover and vegetation type was determined using the USGS land cover surveys (30m resolution) interpolated on our high resolution computational grid.  of the dunes. These were manually adjusted when converting the land cover type into Manning roughness as described in Section 3.1.4. Figure 1.3(b) shows the spatial distribution of Manning roughness (n) values as specified in XBeach simulations. Figure   1.3(c) shows a typical pre storm transect across the dune at the site ( Fig. 1), showing a Manning coefficient varying progressively from 0.025 in the sandy beach foreface (green), to 0.03 when the herbaceous dunes appears (orange), to 0.05 in the shrub area in the back dune (blue), and finally to 0.035 when Spartina Alterniflora appears in the wetlands at the backside of the barrier (yellow).

Computational Grid
To construct a computational grid for XBeach, bathymetric data in the surfzone was required. These were obtained from the NOAA DEM database. The publication date of the data is December 2015, and therefore some uncertainty of the pre-storm bathymetry remains. These were combined with the topographic data described in section 4.1. The final grid has a resolution of 2x2m near the dune area, and 2x5m offshore and is shown in

Boundary Conditions
The start date of the simulation was October 28, 2012 00:00:00 UTC, and the end date was October 31, 2012 00:00:00 UTC. The hydrodynamic boundary conditions i.e.
waves and water levels were obtained from available simulation results of a coupled model (SWAN and ADCIRC) at the offshore boundary location. A peak significant wave height of H s = 6.5 m, spectral peak period T p = 16, mean direction θ = 170 degrees in nautical convention and a peak water level of 2.7 m were specified at the offshore boundary location of our computational grid during landfall. In addition water level was

Results of Field Study for Scenario 1 and 2 Using Manning Approach
Scenario 1 and 2 (Table 1.2) are compared to field survey measurements. In both cases the vegetation effect is modeled by a spatially variable Manning coefficient. In Scenario 1, the model is run in SurfBeat (SB) mode while in Scenario 2 the model is run which did not prevent overwash and inundation. This is illustrated by Figure 1.6, which shows the maximum water depth throughout the simulation for each location, identifying the regions that experienced overwash and consequently became inundated. Both scenarios gave similar prediction of these locations, however SB mode predicted higher maximum landward distance than NH mode. This however could not be validated since no water level measurements were taken. Therefore model performance will be based entirely on morphological predictions.

Comparison of simulation results and measurements indicates various morphological
behaviors and sediment transport sequences occurring at different regions of the site.  Although Scenario 1 (SB) provides a better estimation of dune crest elevation and identification of sensitive locations than Scenario 2 (NH), the NH approach seems to predict more accurately than the SB approach the maximum water level and therefore the maximum landward distance of sediment deposition. This is illustrated in Figure 1.8, which shows overestimation of the maximum landward distance for sediment deposition in SB mode as compared to measurements. This better representation of the maximum water level in NH mode is expected since, by definition in the NH mode, the hydrodynamic module is phase resolving and simulates the storm in real time while the SB mode is phase averaged. However, while the NH mode estimates well the sediment distribution pattern, it underestimates the eroded volume from the dune. At the opposite, the SB mode predicts a shorter deposition distance of the eroded sediment volume resulting in different sediment distribution patterns. Survey measurements show that erosion at sites experiencing a collision regime only, results in relatively less deposition in front of the dune foreface and beach foreshore than further offshore where most of the accretion occurs. By contrast, in areas experiencing overwash and inundation some sediment deposition can be observed on the beach foreshore. Though the uncertainty of the exact bathymetry after the storm however still remains questionable due to the time lag between the storm and post measurements.  The model performance was evaluated for three of Sallenger's regime, collision, overwash, and inundation. The site was further subcategorized into 3 areas, which were assessed separately. Since overwash and inundation regimes were prevented in the NH mode, its performance in the collision regime was evaluated and compared to the SB mode. The northern region having larger of dune elevations was subcategorized as Area A. The middle region, which experienced crest lowering but not landward sediment deposition was categorized as Area B. The southern part and the most sensitive part was categorized as Area C, is where overwash and inundation regime have occurred. The subcategorization of the site is illustrated in Figure 1.10. A summary of the skill parameters of the model used to assess accuracy for each scenario and each subzone A, B, C is presented in Table 1.4, with the definition of these summarized in Table 1 SB was able to predict the overwash regime in these areas, predictions were rated poorly.
In zone C where inundation and landward sediment deposition occurred, SB gave better morphological predictions with a BSS of 0.19 as compared to -2.30 for NH mode. This region consists of lower dune elevations, which were wiped out during hurricane Sandy 2012. SB predicted elevation of dune crest lowering reasonably well as shown in Figure   1.7, however exact morphological profiles were not obtained. In NH dune crest elevation remained high during the simulation and prevented actual changes of the profile. The average results of each skill parameter for all zones are listed in Table 1.5. A 1-D model was therefore constructed similar to a wave flume tank, but with a constant depth chosen to be relatively to the incident wave height. A vegetation field was applied at the bottom surface starting at 250m after the seaward boundary. The length of the vegetation field is 700m, and was chosen to be at least twice the wavelength of the maximum incident wave. A slope of 1:20 was specified at the landward side of the domain to dissipate incoming waves. Morphology and sediment transport modules in XBeach were turned off to simulate wave propagation and dissipation by vegetation only.
A cross-sectional plot of the final model set up is shown in Figure 1.11. The physical characteristics of the plant are given in Table 1   Results are presented in Figure 1.12 in a non-dimensional formulation of the distance required to dissipate 50% of the wave height based on a relative water depth k p d versus X50/L p , with k p and L p the peak wave number and wavelength respectively, d the   Smith et al. 2016, Surfbeat model, Non-hydrostatic. Colors refer to Incident wave periods (black=3s, red=4s, green=5s, blue=6s, magenta=7s, cyan=8s, yellow=10s, Dark green=12s) A similar trend to Smith's experiment, with however still some disparities: the SB mode systematically estimates more vegetated distance to limit the wave height while the NH predict that a shorter distance of vegetation can limit the wave height, and is therefore more sensitive to vegetation.

Numerical Experiment 2
In order to validate these approaches for modeling dune erosion during overwash and inundation regime, the results of a conceptual study is presented. In this study Scenarios 3 and 4, which uses the vegetation module in XBeach are compared with Scenario 1 and 2, which uses a parameterized friction coefficient. For this experiment Scenario 1 is considered as a baseline to assess results of the vegetation module, since it was found to give reasonable morphological predictions and estimation of dune crest elevation and erosion in Area C, which experienced inundation regime during Hurricane Sandy simulation (Table 1.

4). A typical cross section of a barrier island covered with
Spartina Alterniflora will be submerged under water for different depths and subjected to an incident wave with Hmo = 5m, T p = 10s. The correlation of the drag coefficient as presented in "Numerical Experiment 1" for this plant will also be used in this study.
To compare the approaches a Manning roughness coefficient was required to be selected. The "n" value selected for Spartina Alterniflora was 0.04, which was interpolated from  This can be illustrated by comparing the Maximum GLM and orbital velocities that occurred at each relative depth. Figure 1.14 shows that in a relative depth close to zero (i.e. at the dune crest) whether using SB or NH, similar estimates of the maximum GLM velocities are obtained. When the dune crest is submerged (relative depth < 0),

Application to Hurricane Sandy Simulation
Numerical experiment 2 indicated that minimal dune erosion occurs when using the vegetation module in XBeach. However prediction from Scenario 3 and 4 could not be validated directly using the field data set on hurricane Sandy (2012), since it was not possible to apply it to the field study. This is due to lack of information on the physical characteristics of the vegetation, and more importantly, the drag coefficients for each Results for bed level change after hurricane Sandy simulations for Scenario 3 and 4 representing vegetation module predictions are presented in Figure 1.15. As stated, the vegetation module predicts less erosion at the dune crest where vegetation is presented.
Although Scenario 3 was able to predict location of sensitive areas, Scenario 4 using the NH did not. For Area A which experienced a collision regime only during sandy, and no vegetation interaction occurred, similar estimation of bed level changes was obtained for both SB and NH.  Accurate prediction of dune crest lowering but with error greater than 0.5m.

3
Green Accurate prediction of dune crest lowering with error less than 0.5m. Excellent

Discussion
The results of this study allowed us to explore the performance of 4 different approaches to modeling the effect of vegetation on sediment transport and dune erosion during the 3 storm regimes identified by . Significant differences were found when these approaches were compared for each regime. In this section the sensitivity to both wave model and vegetation model are discussed for each storm regimes. Further elaboration on the models capabilities and limitations are provided.

Collision Regime
The ability of XBeach to model the morphological changes of beaches, barrier for Area A with the lowest BSS rated as "reasonable". It must be noted also that discrepancies between the pre and post survey measurements were noticed in the vegetated areas due to filtering errors. It is therefore expected that these errors lead to much lower BSS and SK values.
Sensitivity to the wave model showed that better morphological predictions of the dune profile is obtained when using SurfBeat. Although this cannot be interpreted from the evaluation assessment since the NH gave better BSS and SK, it is clearly illustrated in hence sediment deposition resulting in better morphological skill, while the profile shape and crest elevations of the dunes were poorly predicted.

Overwash Regime
An important skill required for a model to predict flooding events is the capability to accurately estimate the dune crest elevation during a storm regime. This feature is most dominant when applying Scenario 1, which showed highest performance index for predicting dune crest elevation. Again since higher eroded volume occurred in Area B due to overwash, more deposition occurred in front of the dune leading to lower BSS.
Nevertheless, the post-storm dune profile shape and crest elevations from Scenario 1 were the most accurate.
Despite more accuracy in water level, NH mode underestimates the eroded volume in collision and overwash modes. While it simulates well the landward transport distance of sediments in the over-wash regime, the volume is underestimated leading to a crest elevation remaining too high, which further restricts overtopping, and overwash. This was further observed when applying Scenario 3, and 4 that predicted less dune crest elevation, resulting in a lower index score when compared to Scenario 1 and 2 respectively (Table 1.8), and consequently less flooding and landward sediment transport.

Inundation Regime
Area C showed to be the most sensitive region with low dune crest elevation. Many of these areas became inundated and extreme morphological changes occurred. Best morphological predictions were obtained when using SB with a Manning roughness value as suggested by Wamsley et al. (2009). This has resulted in better BSS and estimation of landward sediment deposition volume compared to other scenarios. In addition the dune crest elevation was predicted with very good accuracy and an overall performance index of 2.63. All other scenarios failed to accurately predict the rate of dune erosion, and therefore landward sediment transport leading to inaccurate morphological profile predictions and "poor" BSS.

Vegetation Module in XBeach
The sensitivity of sediment transport to the vegetation formulation introduced by (2004) were explored through Numerical Experiment 1 and 2. In

Numerical Experiment 1 confidence of the hydrodynamic conditions and wave
propagation across a vegetation field were obtained when XBeach results showed similarities to Smith et al. (2016) study. In which, the suggested distance required to dissipate 50% of the incident wave heights were predicted by XBeach.
While the kinematic and dynamic free surface conditions in the presence of integrating the velocity profile to calculate the depth-averaged velocities that are used in the sediment transport module, a significant reduction is introduced leading to much lower estimation of the sediment transport rate.
The vegetation formulations implemented in XBeach and described in Section 3.1.3

assumes rigid cylindrical structures and neglects the sway motion of the plant. For
XBeach this is seen to be ideal, since the computational grid consists of a single vertical layer, and hence it is not possible to resolve the oscillatory deformation of the plant. More profound equations were introduced which accounts for the displacement of the plant (Dupont et al., 2010). These are used to calculate the relative velocity between the water particle and the plant stem and determine an appropriate drag coefficient. Thus, plants with flexible stems will have lower relative velocity, and hence smaller drag (e.g. Stratigaki et al., 2011;Mendez et al., 1999;Maza et al., 2013). Mazda et al. (2013) showed that when the formulation is extended to include the sway motion of plant, the correlation of "bulk drag coefficient" becomes different than when no swaying. In the case of swaying, the drag increases compared to no swaying.
The sensitivity of these equations to the changing flow regimes during a storm event is not considered in which, only a constant drag coefficient value is chosen based on the average flow regime during a storm. In such case the evolution of the storm regime becomes restricted. As a result better predictions were obtained during collision in where the vegetation module restricted the overestimation of dune eroded volume in SB mode, but prevented overwash and inundation in the same time.

Bed Friction Parameterization
Scenario 1 and 2 gave best performance in overwash and inundation when using the Manning roughness coefficients suggested by Wamsley et al. (2009). In addition the estimated dune eroded volume and the landward sediment deposition volume were accurately predicted by Scenario 1 in the inundation regime. The capability of the Manning roughness to describe the sediment concentration has also been well validated in other studies (e.g. . When using Manning roughness less energy dissipation is introduced to the profile if compared to the vegetation module and therefore higher depth average velocities leading to better sediment concentration estimation. However this doesn't necessarily mean that accurate wave propagation was obtained when using the Manning approach, since higher land ward distance of sediment deposition occurred in areas experiencing inundation.

Non-hydrostatic mode
Low performance was observed when using the Non-hydrostatic mode, mainly

Surfbeat mode
The main capability of Surfbeat is to simulate bed level changes and morphological changes of dunes, and barrier islands during Sallenger's Storm regimes in a feasible fashion. It remains the only model to accurately predict the dune crest elevations in the sensitive areas of our site. Using the recommended parameters suggested by literature and described in Section 3.1.6 proved to show reasonable results. Although the seaward deposition of the dune eroded volume is not accurately simulated, due to unresolved phase velocities. Comparison of simulation results and measurements indicate that further offshore deposition occurred of the eroded volume. Inaccurate predictions of sedimentation pattern in front of the dune, has therefore leading to lower performance scores.

Conclusion
The functionality of a dune system to protect coastal communities against hurricanes and seasonal storms, lead US federal and state agencies to assess the sustainability of these structures. Particularly in Rhode Island, dune crest elevation remains the major factor to reduce damage levels associated with flooding and overtopping during overwash and inundation storm regimes (e.g. . Further studies on climate change and sea level rise scenarios indicated severe dune erosion may occur during a 100 year storm event, leaving many coastal communities vulnerable to austere environments . Except for areas of dense vegetation, these showed more stability of the dune system.

Recent formulations describing the energy dissipation rate due to vegetation created
more profound ways to model sediment transport in vegetated areas. The lack of comparative analysis with field data sets however, remains uncertain and questions the reliability of these models. This study evaluated the performance of 4 scenarios to modeling the effect of vegetation on sediment transport and dune erosion during  storm regimes. Combination of wave models and vegetation models were used to simulate the morphological changes occurred at the study site during Hurricane Sandy (2012). The capability and limitation of each scenario to model collision, overwash and inundation regime were identified when compared to pre and post storm measurements of the bed level.
Sensitivity to wave models indicated better performance when using the long wave resolving model in XBeach. Scenarios using "Surfbeat" mode led to overall better prediction of the post dune profile shape and crest elevation allowing the development of the storm regimes. Offshore and landward deposition distances however where inaccurately estimated leading to low performance scores. These were better predicted when using the wave-resolving model. Although, incompatibility of the default sediment transport model implemented in XBeach did not allow proper assessment of the Nonhydrostatic mode.
Depth averaging introduced large modulation errors to the GLM velocities when using the vegetation module in XBeach. While vegetation formulations simulates well the wave energy dissipation, and hence wave propagation, using them explicitly for modeling sediment transport resulted in poor performance scores. Manning roughness underestimates the wave dissipation, and predicts higher depth average velocities eventually leading to better morphological predictions.
In conclusion, Scenario 1 showed to be the most applicable approach for evaluating the stability of a dune system during the storm regimes. Dune crest elevations can be accurately predicted by proper calibration of the model. Such feature makes it a feasible assessment tool for determining dune stability and identifying vulnerable areas.

Introduction
Hurricanes and seasonal storms that have increasingly reached and damaged the U.S Atlantic coastline have driven the US federal and state agencies to impose innovative solutions for coastal protection and beach erosion against future events. In a parallel struggle for environmental preservations and habitat restoration, coastal protection measures have evolved to more soft natural based solutions (e.g. Smith et al., 2016). In particular in Rhode Island, the Coastal Resource Management Council (RI CRMC) has rebuked structural shoreline protection in beaches, dunes, and barrier islands and relies yet on hefty beach and dune nourishment programs.
Over the past decade an increasing number of studies have been addressing coastal hazard and coastal resilience at a global scale (e.g. Baquerizo and Losada, 2008;Vitousek et al., 2017) or more specifically, at regional and local scales. In RI, in particular a series of studies are devoted to assess the risk of coastal communities and barrier islands against extreme storms, sea level rise and climate change (Grilli et al. , 2017cSpaudling et al., 2017;. Results from wave and damage risk models indicated that severe flooding and large damage indexes for residential and commercial structures are expected for the annual 1% of exceedance event (the "100- year" storm). Results of these studies suggest that dune elevation and residential First Furnished Floor elevation (FFFE) are the main factors contributing to reduce the current damage risk of the local coastal communities. Further investigations to evaluate the stability of the dune system during extreme events using a 2-D morphodynamic numerical model has demonstrated the critical importance of a healthy dune and in particular the importance of the vegetation to limit the coastal barrier beaches' erosion and to protect the coastal communities .

Tekmarine
Indeed, while there is a large number of detailed studies investigating GSC's stability based on numerical modeling with validation from flume/tank experiment, there is a lack of validation of the GSC's deployed in the field, simply because there are very little data available. In addition, the complexity of the structures shaped as breakwaters but formed of deformable elements, each behaving differently according to its position in the structure and the specific of the environmental forcing (e.g. up-rush, down-rush) led recent studies to question the applicability of the fundamental equations used for the design , 2009a, 2009b, Dassanayake and Oumeraci 2012.
The potential use in Rhode Island of GSCs structures as dune reinforcement to protect vulnerable structures against austere environment is currently being considered.
Post-construction surveys have been undertaken to assess and validate available design codes and to assess the structure performance under future storms. In order to optimize the future design of the GSC's structure in RI as well as to predict its behavior and resilience in the current wave climate, including extreme storms, as well as in the predicted close future climate (25 years), we propose to study the GSCs' behavior recently deployed southerly from RI on the New England shoreline, in Montauk, NY. We propose to model the environmental forcing and verify the state-of -the art hydraulic stability's formulations originally derived from in-situ's experiments (Dassanayake and Oumeraci, 2012). We use the 2-D hydro-morphodynamic model XBeach  to simulate recent storms and their associated erosion at the Montauk GSC's site. Field data were indeed collected at the site for the recent Tropical storm Hermine (Septemner 2016) 6 months after completion of the GSC's structure. The main steps of the study are: [1] Validation of XBeach using Tropical storm Hermine to simulate the erosion at the site.
[3] Development of Numerical model that predicts the damage level of a dune structure reinforced with GSCs.
While the present approach is deterministic, in a phase 2, we address the damage assessment in a stochastic approach providing fragility curves for the structure as deployed at the site. A description of the site is presented in Section 2; in Section 3, we describe the methodology. Input data required for model and equations are described in Section 4. Results are provided in Section 5, followed by discussion and conclusion in Section 6 and 7 respectively.

Site Description
The increase in storm frequency combined with the urbanization of the dune prevents the natural restoration of beach systems and limit their role of protective buffer zones along the sandy coastlines of Long Island, resulting in an increasingly higher risk imposed on coastal developments and insurance agencies . To improve coastal resiliency of the downtown area of Montauk, NY and to minimize the risk sustained by coastal front developments, the USACE has built a dune system reinforced with GSCs. The objective of the new synthetic dune system is to limit the beach erosion and protect the dense community settled at the top of the dune barrier.  The 800m stretch of the GSC structure expands at the top of the beach parallel to the line of coastal front development. The project consumed 11,000 non-woven GSCs of typical dimensions (1.67x0.9x0.3m) using 109,000 cubic yards of sand with similar characteristics to that found in site for 80% filling ratio. GSCs were placed longitudinally using two armor layers at a slope of 1:2. The structure is 3 meter high and has a crest elevation of 3.8m NAVD88. The design of the structure is shown in Figure 2 Post construction surveys recorded several consecutive storms that hit the project.
Tropical Storm (TS) Hermine (September 2016) hit the site 6 months after completion followed by a Nor'easter 4 months after. These storms generated high waves with offshore significant wave heights of 4 and 6 m respectively (NOAA), causing severe erosion of the beach's berm and exposure of the GSC layer. Post analysis of the site revealed that the cyclonic nature of storms causes erosion in the east side of the site and deposition on the west side. The middle region is considered to be the most sensitive for the local economy with 3 major hotels located in this area. Post survey photographs are shown in Figure 2.3 for this region. (USACE-NAN, 2014). c) Aerial photography taken during construction of the project (Photo credit by TenCate).

Methodology
In this study we investigate the adequacy of available hydraulic stability equations for predicting damage levels for reinforced dunes with GSCs in barrier beaches during a storm regime. Field survey performed after TS Hermine provide data to validate the models used in the study. A brief description of formulations and models used is presented in this section. While this study is part of a larger research project seeking to develop fragility curves for the GSC based on a stochastic approach of the hazard, the objectives of this component of the project are: [1] Validate XBeach to predict the dune erosion at the GSC site using field survey performed after TS Hermine damaged the site.
[2] Investigate the adequacy of state of the art's stability equations for predicting damage level on reinforced dunes with GSCs (Oumeraci, 2012) by comparing with field observations.
[3] Assess if these models and approaches are adequate to be used for phase 2, which is developing a fragility curve to estimate the damage levels associate with reinforced dunes with GSCs.

Damage levels of Reinforced Dune with GSCs
Damage levels refer to Dassanayake and Oumeraci's classifications (2013) in which damages are classified into two sub-classification, I classification 1, damages are identified at the scale of a single GSC, while in classification 2 damages are identified at the scale of the entire GSC's structure. Some authors (Shirlal and Mallidi, 2015) summarized the results in a level 3 classification, where the total fraction of the damage structure is estimated as the percentage of GSCs displaced with respect to total GSCs in the structure. We introduce a level 4 classification where the full beach is considered and the erosion on the entire ecosystem beach and GSC's structure is considered (

Hydro-morphodynamic Model XBeach
The XBeach ("eXtreme Beach behavior") software package was selected to quantify damage levels 0,1,2,3,4 and 6, as well as required hydrodynamic conditions to calculate hydraulic instabilities of GSCs for damage level 5. XBeach is a 2-D coupled hydrodynamic and morphodynamic model that dynamically simulates the coastal response during time-varying storm conditions, and consequently the storm-induced changes in bed level . The model is designed to simulate processes occurring during the four erosion regimes defined by , swash, collision, overwash and inundation, on coarse-grained beaches, dunes and barrier islands.
The key advantage of using XBeach for this study is its algorithm for calculating dune avalanching and ability to resolve long waves and ultimately the swash motions that dominate the sediment transport at the foreface of the dune. Furthermore, it offers a choice of hydrodynamic models and sediment transport models combined to a morphology change model. Wave and sediment transport modules selected for this study are briefly described.

Wave Propagation and Transformation Module
The "Surfbeat" model was selected which simulates mean currents in combination with a wave action conservation equation. The characteristic feature of this model is the modeling of the long wave component of swash motions, the surf-beats resulting from the forcing by wave groups of infragravity waves and causing slow oscillations of the mean water level directly related to foreshore erosion (Longuet-Higgins and Stewart, 1964;. Surf-beats are indeed particularly dominant in the surf zone in storm conditions, since the high frequency band is saturated, while the low infra-gravity frequency band is not, particularly on dissipative beaches . While the full directional distribution of wave action is maintained, the frequency spectrum is limited to a single frequency, the peak frequency. The propagation of wind waves is based on the conservation of wave action equation, coupled to a roller energy balance equation  in the surf zone through the wave dissipation term.
Radiation stress tensors acting as forcing terms for depth-and period-averaged mean Although non-linear properties of short waves which are not represented by linear wave theories, such as asymmetry, skewness and turbulence and are only parameterized; the surf-beat approach has been validated on dissipative beaches with a surf similarity parameter ranging from 0 -0.55 (mild conditions) (e.g. . Indeed on dissipative beaches the short waves are mostly dissipated by the time they reach the shoreline, and only infragravity reach the shoreline.

Sediment Transport Module
The sediment transport model used in this study is the default model in XBeach  . The sediment transport is controlled by the sediment concentration in the water column relative to an equilibrium concentration (Soulsby, 1997). While wave propagation equations are based on linear theory, the effect of non-linear waves on sediment transport is included by adding an arbitrary parameterized advection velocity u a (Stokes drift) to the Eulerian velocity, based on a wave skewness S k and asymmetry A s , both weighted by arbitrary coefficient, f s and f a (facua parameter) .
Higher facua value will result in more nonlinearity of the wave shape and hence more landward sediment transport.

Model Parameters
Default model parameters were used, except for those where extensive set of validation and recommendations were provided through comparative field analysis and laboratory experiments Devet et al., 2015;. Such parameters are bed friction coefficients, wave skewness and Assymetry "facua" parameter, morphodynamic accelerator factor "morfac", critical wet and dry slope avalanching value, and scaling factor "eps".

Hydraulic stability equations
While early GSC-structures were designed using the hydraulic stability formula for rubble mound armor layers such as Hudson's formula (1956), this formulations were not appropriate for deformable structures in marine environment and were progressively replaced with new semi-empirical formulations based on tank/flume experiments using specific GSC arranged in selected structural shape within a controlled wave environment.
Wouters (1998) using experimental data (Bouyze and Schram, 1990), developed Hudson's equation to account for wave period and steepness using the surf similarity parameter, as well as using different characteristic length to account for geometric variation of the structure's elements. Using Wouters' formula (1998) Oumeraci (2003, proposed new stability criteria for GSC's placed either on the slope or on the crest of the structure, expressed as: where N s is the stability number of the element, H s is the significant wave height at the toe of the structure, ! = 1 − • ! + ! • is density of the GSC element , ! and ! is the density of water and sand respectively, = • is the thickness of the first armor layer, is the length of the GSC, α is the structure slope, ! = locations of the GSC (i.e. slope or crest).  provided the first "process-based" stability equations for GSCs. Such equations are presented hereafter for two failure modes, sliding and overturning: where l ! is the length of the container, u is the horizontal flow velocity. KS !" , KS !" , KS ! , KS !" are deformation factors for drag, lift, resisting, and advection forces respectively.
Similarly C ! , C ! , C ! are coefficients for drag, lift, and advection respectively. µ is the viscosity of water, and g is gravitational acceleration.

Assessment and Evaluation Methods
The ability of the model to simulate accurately the eroded volume is assessed applying several skill parameters. A summary of the skill parameters of the model used to assess the accuracy of the model for each scenario is presented in

SK
Relative value of the eroded volume

Data and Model Setup
Bathymetric, topographic and hydrodynamic data were needed to simulate TS Hermine, the first storm to hit the project 6 month after completion. These were obtained from various sources; and processed as input data for our 2-D numerical model XBeach.
Further description of required input data and model set up can be found in their manual (XBeach Manual). This section only gives a brief description of our model set up and data used.

Topographic Data
Pre topographic conditions for the TS Hermine simulation used several data sources.
As-Built drawings and survey measurements provided by the USACE were used to construct the backshore of the beach and the cover layer of the GSCs. These Post survey measurements of TS Hermine were used to validate XBeach results.
These were taken on September 9, 2016 one day after the storm has dissipated. Using RTK-GPS, 22 transects along the project were measured. These are also provided in similar coordinate systems to pre survey measurements.

Geotextile Sand-filled Containers (GSCs)
The GSCs' locations were required to specify them as a non-erodible layer in XBeach. They were also specified as save points to output velocities and significant wave heights and to ultimately determine the hydraulic stability of the GSCs (Equations 2.1 -2.4).

Computational grid and Model Setup
To have a full XBeach model domain, bathymetric data was combined with the topographic data described in section 4.1. These were obtained from National Ocean and Atmospheric Agency (NOAA) Digital Coast database available on their website; and covers the full area domain. The resolution of the data is 1x1m. This was further interpolated to 5x1.2m offshore and 2x1.2 near-shore. The bathymetric data was published in 2015, and therefor some differences in the bed level at the interface of each data set (topo vs. bathy) was observed when combining the data. In addition several offshore bars were observed in the data and were removed to avoid un necessary dissipation of the incident waves.

Hydrodynamic Boundary Conditions
Hydrodynamic Boundary conditions at the offshore location of our computational Energy density spectrums were obtained from NOAA's wave buoy in 1D. A directional spreading function based on Mitsuyasu-type (Mitsuyasu et al. 1975) was applied to generate a 2-D wave spectrum used in input to simulate the wave propagation    Due to the uncertainty in the data for both (1) STWAVE Results due to the large discretization of the computational grid which is expected to result in lower peak wave heights, as well as the (2) bathymetric conditions used for the XBeach computational grid; certain adjustments to the hydrodynamic boundary conditions were made to account for these uncertainties and to obtain an upper confident level of the results. Especially when XBeach is shown to be very sensitive to the significant wave heights  and water depths. Any under specification of the offshore wave conditions might lead to lower estimation of erosion. Therefore to ensure the development of Damage level 3 (exposure of the bag) and validate the hydraulic stability equations (Section 3.3) for predicting Damage level 4, the significant wave heights at the offshore boundary were increased by 25%. In addition, the water levels obtained from the NOAA offshore Buoy were increased by 10% to account for uncertainty in the near shore bathymetry which may lead to additional dissipation of the incident waves.

Results
First, results of predicted bed level changes are compared with post survey measurements using standard morphological skill parameters. Then, predicted damage levels (Section 3.1) are assessed using both Oumeraci (2003) and Recio and Oumeraci 's stability formula (2008).

TS Hermine Bed Level Changes
Bed level change analysis of the XBeach results indicated intense erosion of the berm as well as the cover layer of the GSCs. Since TS Hermine hit the coastline in collision regime, avalanching was the principal process responsible for the beach morphological changes resulting in cover layer and berm's erosion and offshore sediment transport, with rapid deposition in the beach backshore (Figure 2.10). Erosion and deposition shows a relatively homogeneous behavior along the entire beach section. A detail comparison of measured and simulated erosion shows a spatial variability in observed eroded volume along the beach, reproduced with variable accuracy by simulations (Figure 2.11). Indeed measurements show an asymmetry in the sediment erosion/acretion pattern associated to a slight asymmetry in the coastal topography (east side is relatively higher) and to a likely longshore current, resulting in larger sediment accretion at the west side of the site. The middle area (in front the 3 major hotels) shows the most intense erosion of the berm, causing exposure of the GSC layer. Although numerical simulations do not reproduce the details of this observed asymmetric pattern, they capture relatively closely the erosion observed in the middle region where exposure of the GSC occurred (red line in Fig. 2.11).
The beach is schematized in 3 sub-zones characterized by their general pattern in erosion/accretion as well as model's performance (Figure 2.12). Zone A, the west side of the site where accretion was observed but however not predicted; indeed the model predicts exposure of the bag. Zone B, the middle region where GSC's exposure is observed, as well as predicted. Zone C, eastern side of the beach where only the erosion of the berm was observed and predicted. The model's performance is assessed for each zone based on standard model's skill parameters averaged for each morphological sub-zone (Section 3.4). Results are summarized in Table 2.3. Zone A show a "poor" performance. Indeed the predicted morphological response is inaccurate with large errors and low skill score; An "Excellent" performance is obtained in Zone B, with an average of the Brier Skill Score (BSS= 0.83), indicating good morphological predictions. Eroded volume was accurately calculated resulting in a high Gallegher Skill Score (SK=0.63) and low Bias Score (BI=0.18 m). A Reasonable performance is achieved in Zone C, predicting correctly the non-exposure of the GSCs, although the eroded volume is over estimated, leading to low scores. Typical cross section profiles of the bed level change for each zone is provided in Figure 2.13, to further illustrate the performance of the model.

TS Hermine Damage Levels
Results of the simulations are used to predict damage levels ( Table 2.1 at each location of the GSCs' structure. Bed level changes are used to determine damage level 0-3. For GSCs' exposed locations hydraulic stability equations (Section 3.2) are applied to determine if damage level 4 is likely to occur. Damage level 5, is only calculated for GSCs located along the structure's toe. Maximum water levels are used to determine if damage level 6 is likely to occur. Figure 2.14 provides the estimated damage levels based on Oumeraci (2003) and 's hydraulic stability equations.
Hudson's based stability formula proposed by Oumeraci 2003, did not give any predictions of GSC instability along the entire length of the structure. Recio and Omeraci (2008)'s stability equation, predicts 1 instable GSC occurring during peak conditions.
Although 1 displaced GSC was also observed during survey measurement (Figure 2.3), the location of the predicted displaced bag is not accurately predicted. Furthermore Damage 5 and 6 did not occur during TS Hermine, and was therefore not predicted by the model.

Hydraulic Stability of GSCs
Stability equations (both formulations, Oumeraci (2003) and ) are used at each exposed location to assess the GSCs' stability and the performance of the structure, providing a quantitative method to identify vulnerable locations. The Hudson's based stability equation (Oumeraci, 2003) suggests that all of exposed locations were still in stable conditions during the storm. Figure 2.15a plots the maximum calculated stability number at each location during the storm. It can be seen that these remained well under critical conditions. This however was not suggested by vulnerable during the storm. Two of these locations are located in Zone B, which is the most sensitive region. Figure 2.15b shows the suggested minimum required GSC's length to prevent displacement of each of the exposed GSCs. Except for 1 location, all the GSCs showed to have sufficient length to prevent displacement.

Discussion
The results of this study allowed us to investigate the adequacy of a selected hydromorphodynamic model in combination with two GSC hydraulic stability formulas to determine the damage level associated with a reinforced dune with GSC during a storm regime. Damage states associated with the structure were defined and calculated. Further elaboration and discussion of the results is provided in this section.

Performance of Hydro-morphodynamic Model
The ability of XBeach to model the morphological changes of beaches, barrier islands, and dunes have been well validated for collision regime through laboratory and field experiments Elsayed and Omeraci, 2017). Good performance was achieved particularly for the region of interest where the GSCs were exposed (Zone B), where accurate morphological predictions were obtained. The simulated eroded volume in Zone A and C is however overestimated by the model, and the longshore sediment transport and deposition on the west was not predicted.
This inaccuracy is expected since XBeach's "Surfbeat" mode has a phase averaged wave module and hence instantaneous velocities are not resolved. Let's also note that the temporal variability of the bathy-topo induces a large uncertainty in the model initial conditions, which can lead to significant variations in the extend of the swash zone, and consequently in the morphological changes of the beach.

Performance of Hydraulic Stability Equations
Hudson's based stability equations seem to underestimate the critical stability of the GSC. Several locations became vulnerable during TS Hermine, which was not detected by Hudson's based stability equations. These equations estimate the stability as a function of the significant wave height at the toe of the structure. Throughout the storm regime, the water depth at that location remained very shallow, and thus most short waves have The semi-empirical stability formulas such as introduced by Recio and Omeraci (2008) however suggested possible instabilities during TS Hermine, which was verified during the post storm surveys. Although the applicability of these equations to this study remains questionable, since Recio and Omeraci (2008) hydraulic stability equations have been validated for structures with 45 degrees slope (i.e. 1:1). The structure at the project site was constructed with a slope of 1:2, thus giving less contact area around the perimeter of the GSC. Different force coefficients calibration might be required for a 1:2 slope since this increases the effective area and hence the resultant forces and moments.

Conclusion
Assessment tools are required to determine the fragility of a reinforced dune system with GSC during a storm regime. This study presented a numerical model that identifies the damage states associated with these structures. The model combines an existing hydro-morphodynamic software package "XBeach", with two available hydraulic stability equations for GSC. A Hudson's (1956) based stability equation proposed by Oumeraci (2003) was compared with a more semi-empirical equation proposed by .
XBeach model was used to determine the initial damage states of erosion, while the hydraulic stability equations were used determine the stability of exposed GSCs. The model was validated against post storm measurements taken after TS Hermine, which caused exposure and minor damage to the GSC at the USACE project site. Results of bed level changes and eroded volume showed excellent performance in the region of interest, while some discrepancies occurred in other regions. Appropriate damage states as occurred during the storm was identified by the model, and consisted mainly of berm and cover layer erosion up to exposure of the GSC layer. The study indicated better confidence was obtained when using Recio and Omeraci (2008) stability equation for the exposed GSC, however the validity of this equation to the configuration of the structure used in this study remains questionable.