A Comparison of Wave and Erosion Modeling Methods for the 100-Year Storm in Southern Rhode Island

This body of work consists of two manuscripts related to modeling ocean waves and erosion associated with the 100-year storm in southern Rhode Island. Predicting the 100-year storm inundation along the Narragansett Bay shoreline The lack of confidence in FEMA maps after Hurricane Sandy (2012) led to question the accuracy of the methodology used in FEMA Flood Insurance Rate Maps (FIRM). The present analysis presents a case study, re-computing the 100year inundation maps in Washington County, RI, using a chain of two-dimensional models. Presented results focus on the West Passage of the Narragansett Bay, Rhode Island, USA, an area characterized by a complex shoreline. The selection of the initial and boundary conditions for ocean wave simulations are based on results of fully coupled 2-D surge and wave models, the ADvanced CIRCulation model (ADCIRC) and the phase-averaged STeady state spectral WAVE model (STWAVE), respectively, from the North Atlantic Coast Comprehensive Study (NACCS) performed by the U.S. Army Corps of Engineers. Waves are simulated using STWAVE over a series of high-resolution grids in near-shore and overland areas. The method is referred to as the NAST method. Results are mapped and compared with FEMA’s results along transects. Cross-sections across the Base Flood Elevation (BFE) at the site of the FEMA 1-D transects are in relatively good agreement using both methodologies, with larger discrepancies shown in the northern section of the bay. Discrepancies are primarily due to a difference in the 100-year wind assumption and secondarily to the 100-year storm surge assumption, which, besides defining the depth and extension of the inundation, controls the wave impact over land. Mapping the results reinforces the importance of adopting a 2-D approach to fully represent the inundation, showing that the selected transects often miss the sites of the most extreme wave incursion and most destructive impacts. The analysis highlights the differences in methods and results but also suggests adopting, besides a 2-D approach, a scenario-based approach to the 100-year storm rather than a deterministic single map to assess the uncertainty associated with this event. Modeling the erosion associated with the 100-year storm on the southern Rhode Island coast The erosion due to a synthetic 100-year storm is modeled on the southern shore of Rhode Island using the process-based morphodynamic model XBeach (Roelvink et al., 2009). The study area includes barrier beaches, coastal ponds, and residential areas. The model is forced using a synthetic storm time series extracted from the North Atlantic Coast Comprehensive Study (NACCS; Cialone et al., 2015, Nadal-Caraballo et al., 2015) database, selected as a proxy 100-year storm. An empirical method (Stockdon et al., 2006; Stockdon et al., 2012) is used to parameterize wave setup, swash, and runup to estimate the erosion impact regime (Sallenger, 2000), using NACCS water levels and offshore wave heights as inputs at seven cross-shore transects within the study area. As the storm is expected to evolve from the collision regime to the overwash and inundation regimes over time, the two-step XBeach model calibration process suggested by Nederhoff (2014) is followed. The first step consists of adjusting the facua wave skewness/asymmetry parameter to calibrate the collision regime. The second step involves testing the sensitivity to friction across the dune to calibrate the overwash regime. The sensitivity of the model to the facua parameter is tested by modeling a historical storm that stayed within the collision regime. Hurricane Irene (August, 2011) is selected and simulated erosion is compared to measured erosion at five cross-shore profiles monitored by the University of Rhode Island Graduate School of Oceanography (GSO). GSO performed measurements and data collection before and after Irene’s impact at the site. For the 100-year storm, results of four simulation scenarios are presented: the facua parameter is tested for the default and maximum suggested values of 0.1 and 0.3, respectively, and the bottom friction is either set to a constant Manning’s n of 0.02 or a variable Manning’s n bed friction coefficient based on land cover. The resulting eroded dunes for the four simulations are presented in 2-D maps as well as 1-D cross-sections at the seven cross-shore locations. Two additional transects to the GSO transects are included, located at the Federal Emergency Management Agency (FEMA) coastal transect sites that were used to create the 100-year Flood Insurance Rate Map (FIRM) for the region. Results are compared at these two sites with FEMAs dune erosion protocol. Additionally, the simulated dune is compared to the generalized barrier profile for a 100-year storm developed for the region by Oakley (2015).

but also suggests adopting, besides a 2-D approach, a scenario-based approach to the 100-year storm rather than a deterministic single map to assess the uncertainty associated with this event.

Introduction
The National Research Council (NRC) of the National Academies published a report in 2009 discussing theoretical and practical issues related to the Federal Emergency Management Agency's (FEMA) map accuracy (NRC, 2009). One of NRC's key findings was that the accuracy of the coastal flood maps could be improved significantly through the use of coupled two-dimensional storm surge and wave models and improved process models (NRC, 2009). In this analysis, we present a case study that addresses this particular concern.
Indeed, the severe local impact of Hurricane (2012) extended the method to a dynamical-statistical approach, in which they simulated synthetic storms associated with climate model simulations to account for future storm climatology changes. The latter method was applied to assess the storm surge risk on specific coastal communities (Lin and Emanuel, 2016).These state-of the art research approaches only focused on the storm surge risk.
In this study we use a similar statistical-deterministic approach, using the results of the North Atlantic Coast Comprehensive Study (NACCS; , performed by the U.S. Army Corps of Engineers (USACE), to determine the 100-year storm surge, which we use as initial and boundary conditions for wave simulations performed in high resolution nested coastal grids. In the NACCS, the USACE generated synthetic tropical and historical extra-tropical storms to simulate surge and waves along the U.S. East Coast, using the fully-coupled ADvanced CIRCulation hydrodynamic (ADCIRC) and the phase-averaged STeady state spectral wave (STWAVE) models (Smith et al., 1991;Massey et al., 2011;Anderson and Smith, 2015). Resulting storm surge and wave spectral parameters were saved at hundreds of virtual stations in the study area, and provided in a probabilistic form (return period) .
In RI, FEMA has historically simulated wave elevations over land using the 1-D Wave Height Analysis for Flood Insurance Studies model (WHAFIS), which was developed in 1977 and revised in 1988 to include wave growth and decay through vegetation (FEMA, 1988(FEMA, , 2007b (FEMA, 1991;Stoa, 1978). In FEMA's method, calculations are done along each transect using one of the above formulations and continuous values are estimated along the RI coastline by interpolation in between transects. FEMA's method uses erosion and runup modules when criteria for those processes to occur are met. In addition, when a potentially fragile structure is present on a transect, a new land profile is built, assuming that the structure has failed, thus limiting the protection of areas further inland.
In this study, we built on the results of the NACCS by integrating their regional scale surge estimations with local scale wave modeling. For the latter, we elected to use STWAVE in order to comply with FEMA's authorized model list.
Our chain of 2-D models is referred to in the following as the NAST (NACCS-STWAVE) approach.
Besides the 2-D extension, the main focus of the work is on the impact of waves in the inundation zone. Our method provides the 2-D approach that FEMA's method is currently lacking, providing a 2-D representation of the inundation zone, including wave elevation at a very fine scale (waves propagated in a 10 m grid and mapped on a 5 m interpolated grid), as well as a reliable method to capture specific 2-D processes associated with wave propagation, such as refraction, which is particularly relevant in an area defined by complex bathymetric and topographic features such as the Narragansett Bay (Figure 1.1).
The methodology as applied in this case study carries some simplified assumptions: (1) No protocol to assess the inundation zone in case of failure of a structure is applied; (2) Dunes are assumed kept intact. Although we have developed a protocol that simulates wave propagation across an eroded dune following a 100-year dune profile along the southern shoreline of RI Grilli et al., 2016), we have neglected this process in this more sheltered section of the coastline; (3) since STWAVE is not a phase resolving model, it does not simulate the runup, therefore this process is not included in the simulations. In the following, we first describe the NAST methodology as applied in the case study area; then, we present the resulting new inundation maps as well as selected cross-sections across the inundation zone to illustrate the comparison FEMA/NAST at the site of FEMA's transects.

Methodology
Washington County was divided into several Cartesian computational grids with resolution of ∆x = 10 m, designed to maximize topographic/bathymetric resolution and reflect the local wave climate. Waves are simulated across these grids using NACCS results as initial and boundary conditions.
Characteristics of the Narragansett Bay grid G4 are summarized in Table 1.1 ( Figure 1.1). The grid is oriented along the dominant wave direction of propagation (cross-shore, x) with the origin of the grid (x 0 ,y 0 ) at the south-east corner. The angle of rotation of the grid, α, is estimated counterclockwise from East (Table   1.1). Bathymetry and topography were derived from the merged 2011 LiDAR and NOAA bathymetry (RIGIS, 2013) (Figure 1.1). STWAVE was applied in half-plane mode (energy propagation offshore is neglected) using 2-D incident wave spectra reconstructed from the NACCS parameters as offshore boundary conditions. STWAVE simulates wave propagation in the horizontal plane, including refraction based on geometric optic theory and shoaling based on the conservation of wave action along wave rays (Smith et al., 1991;Massey et al., 2011). The sea state growth through the transfer of momentum from the wind field to the wave field is modeled using Resio's (1981Resio's ( , 1988 formulation. Wind energy fed into the waves is redistributed through non-linear interactions from the peak of the spectrum to both lower frequencies, increasing the peak period, and to higher frequencies, where it is dissipated through whitecapping, depth-induced wave breaking, and turbulent effects (Resio, 1987(Resio, , 1988. This nonlinear energy transfer due to wave-wave interactions, first described by the Boltzmann integrals of Hasselmann (1961), is implemented using Resio's methodology (Resio and Perrie, 1991). The loss of energy by wave breaking is modeled using Miche's (1951) breaking criterion, which includes the effects of both water depth and wave steepness-limited breaking (Smith et al., 1997;Battjes 1982;Battjes and Janssen, 1978). The energy loss is simulated by reducing the spectral energy, in each frequency and directional band proportionally to the amount of pre-breaking energy contained in each band. The bottom friction loss is implemented in this version of STWAVE (V6) using a Manning coefficient formulation (Holthuijsen, 2007), which can be specified as a spatially variable parameter.
The directional wave spectrum for the 100-year storm is specified as input offshore boundary conditions for each STWAVE computational grid, based on results and parameters computed in the NACCS study at save points, including the significant wave height, H s , peak spectral period, T p , and dominant wave direction, θ 0 .
To express the spectrum based on these parameters, we use a TMA energy spectrum since it was shown to optimally represent swell and wind-generated gravity waves in shallow water (Bouws et al., 1985).
The NACCS study used a state-of-the-art stochastic approach to generate a large number ( The 100-year return period values of these parameters (mean and 95% upper confidence interval) were extracted at each save point located in RI waters (about 1,000 points) and interpolated onto our local high-resolution 10 m grids ( Figure   1.1). The mean value of each relevant parameter at the offshore boundary was selected as the representative offshore value used to reconstruct the local TMA spectra. We selected the standard value for the spectral peak enhancement factor γ of 3.3, as well as for the directional spreading factor cos n θ with n = 4. Parameters for the 100-year storm spectrum and grid boundary conditions are defined in Table   1.2, and briefly described in the following.
The term still water level (SWL) is commonly used to describe the water level in the absence of wind waves and their effects, and thus includes the astronomical tide and the storm surge (due to wind effects); the static water level (STWL) is defined as the sum of the SWL and the static wave setup/setdown, the mean additional elevation of the water (positive or negative) associated with the presence of waves. NACCS provided the STWL for each 100-year synthetic storm at each save point; this level was interpolated and used as the reference level in our 10 m grids to compute near-shore and overland wave propagation with STWAVE. It should be noted that, in RI, the CRMC elected to use the 100-year upper 95% confidence interval values, rather than the mean value as the representative 100- year STWL. This selection was made in order to be conservative and account for uncertainties in the analyses.   (Wamsley et al., 2009(Wamsley et al., , 2010Arcement and Schneider, 1989). The BFE is the maximum water elevation of the controlling wave crest riding over the STWL, referenced to NAVD88 (BFE = STWL + η c ). FEMA's VE and AE zones are defined by a wave crest threshold of 0.9 m (3 ft). The VE zone is an inundated area with a controlling wave crest larger than 0.9 m (H s > 0.8 m); the AE zone is an inundated area with a controlling wave crest smaller than 0.9 m.

Results
The limit of moderate wave action (LiMWA) in the AE zone is the limit of wave crests larger than 0.45 m (1.5 ft; H s > 0.4 m). Note that neither dynamic setup nor runup is included in the current estimation of NAST BFE.
Wave simulations were performed for the 100-year spectral parameters, as defined in the previous section, specified as boundary and initial conditions. The model was first tested for a range of grid sizes, friction values, and wind speed and direction to assess the associated uncertainty. Significant wave height results were plotted along a transect SW-NE crossing the Narragansett Bay West Passage from the mouth of the bay up to Prudence Island in the north of the bay (Figure 1.2a).    Figure   1.5 shows results of NAST simulations for scenarios 2 and 3 (wind and extreme wind). We can see that in these extremely windy conditions, the local wind can regenerate significant waves in the bay. However, from a statistical perspective, winds associated with scenario 2 and 3 do not reflect 100-year winds, but much larger return period storms. However, in view of the current observed change in storm frequency and intensity it might be wise to use a design storm associated with more extreme winds than the current 100-year winds.    The wave propagation at the mouth of the Petaquamscutt river at the north of Narragansett beach is finely modeled, showing waves (on the order of 1 m of controlling wave crest) propagating upstream, at the top of the surge across the salt marsh surrounding the river, locally increasing the BFE from 4 m up to 5.5 m (scenario 1; base case); for the high wind case, waves are penetrating slightly further upstream into the river. Wave regeneration due to the extreme southerly winds in the river is also apparent, increasing the BFE on the far shore.

Sensitivity analysis
In Narragansett Beach, the dunes, assumed intact for the current simulations, are creating a strong barrier, forcing the waves to break and preventing them from overflowing inland as they do in the no-dune/parking lot areas. In these unprotected areas, waves with an order of 1 m wave crest are propagating across the main road. While FEMA uses a relatively high density of transects in this area, none of these are capturing the maximum amplitude of the inundation, which crosses the main road and therefore isolates the northern section of the town from the southern one.
Similarly in Bonnet Shores and Wickford, none of the FEMA transects are able to capture the largest waves propagating across the beach in the inundation zone.
In Bonnet Shores, Transect 48 intersects the main road, however, further north of the beach, where overland waves are already reduced. In Wickford, since the harbor is in a protected area, only relatively small waves are reaching the harbor.
However, since NAST's SWEL is higher than FEMA's, waves are likely propagating further overland as shown on Transect 60. Runup might be a significant issue on the wall protecting the marina and a phase resolving model would be recommended to more accurately assess the vulnerability of this area.  It is interesting to note that the discrepancies between FEMA's BFE, as reported in the 2012 report (FEMA, 2012), and NAST's BFE (Appendix 1) seem larger than observed when we compare the transects. The VE zone elevation is indeed defined as the average BFE observed in the VE zone; FEMA and NAST however use slightly different methodologies to define the VE zone offshore limit.
NAST restricts a priori the offshore limit of the VE zone to the contour line 0  2. NAST defines the offshore limit of the VE zone at the limit of the shoreline, the contour line 0 NAVD88. FEMA's VE zone offshore limit is, although not clearly defined, around the wave breaking zone, independently of happening overland or not. The VE zones might therefore be inflated in well protected areas, such as in transect 60, where no waves are propagating overland and where the VE elevation is based on an offshore (in the sense of not overland) limit.
Mapping the results reinforces the importance of adopting a 2-D approach to fully represent the inundation, showing that the selected transects often miss the sites of the most extreme wave incursion and most destructive impacts. The analysis highlights the differences in methods and results but also suggests adopting, besides a 2-D approach, a scenario-based approach to the 100-year storm rather than a deterministic single map to assess the uncertainty associated with the event.
It should be pointed out that the NAST method does not model wave runup nor reflection. We are currently testing a similar approach, but using a phase resolving model in coastal areas to model these processes, although such a model is not currently validated.   An empirical method (Stockdon et al., 2006;Stockdon et al., 2012) is used to parameterize wave setup, swash, and runup to estimate the erosion impact regime (Sallenger, 2000), using NACCS water levels and offshore wave heights as inputs at seven cross-shore transects within the study area. As the storm is expected to evolve from the collision regime to the overwash and inundation regimes over time, the two-step XBeach model calibration process suggested by Nederhoff (2014)

Introduction
Modeling the 100-year event is an integral part of coastal floodplain mapping and risk analysis. The Federal Emergency Management Agency (FEMA) defines special flood hazard areas as regions that will be inundated by a flood that has a 1percent chance of being equaled or exceeded in any given year (100-year flood). In coastal regions, the additional risk due to storm surge and waves is independently assessed for each of those hazards. FEMA defines two main coastal zones based on the expected wave crest elevation of the 1-percent of the highest waves (η c ): the VE zone with η c greater than or equal to 3 feet, and the AE zone with η c less than 3 feet.
Dune erosion occurs over the evolution of the storm, resulting in sediment transport landward, cross-shore, and offshore. Correlatively, the resulting changes in beach profile and coastline shape modify wave propagation and their impact on the coastline, as well as on the expected total inundation elevation. Consequently, the dynamics of the dune erosion can potentially have a large impact on the limit and characteristics of the 100-year flood zones.
Long-term time scale predictions of shoreline erosion due to sea-level rise (SLR) are part of a dynamic field of research. Modeling teams are attempting to assess the impact of climatic change on the potential hazards using either a scenario-based risk analysis (Oumeraci et al., 2015) or dynamical-statistical modeling . Few approaches include the complex process of shoreline changes (Bilksie et al., 2014;Grilli et al., 2016), but these long term erosion processes must be included if one expects the 100-year storm to occur later in the century. FEMA, however, does not currently consider the effects of SLR in the FIRMs. To be consistent with FEMA's current protocol, we assumed that the 100-year storm is occurring in the near future, when the SLR is negligible, and the erosion processes are at the time scale of the storm event only.
The study area is located on the southern coast of Rhode Island (RI) and includes East and Charlestown barrier beaches, Ninigret Pond, as well as the surrounding residential and wildlife preserve areas (Figure 2.1). The following paragraphs describe the two methods used to determine 100-year flood maps for the region, both taking dune erosion into consideration but using different approaches.
The first method is the official FEMA methodology (2013) used to create the RI FIRMs, based on simulations along 1-D transects interpolated along the shoreline.
The second method is based on using a suite of 2-D models and was developed as an alternative to FEMA's 1-D methodology for RI (Grilli et al., 2015;Spaulding et al., 2016;Schambach et al, 2016). PBL; Thompson and Cardone, 1996), along with the regional wave model WAM  Both methodologies, FEMA and NAST, make use of predetermined eroded dune profiles to modify the bathymetry/topography and then run stationary wave models assuming a static modified shoreline. An alternative approach to determining the beach erosion is a dynamic approach, which uses process-based erosion models to simulate the dune erosion over the time scale of a storm.
In order to properly model erosion, one would expect a wave propagation module containing the physics to simulate the swash motion. Holland and Puleo (2001) indeed showed that swash collisions are directly related to foreshore erosion. Swash motions, or surf-beats, result from wave group forcing of infragravity waves (0.004 < f < 0.05 Hz) (Tucker, 1954;Longuett-Higgins and Stewart, 1964). These are slow moving (minutes) oscillations of the mean super-elevated water level (static setup) Svendsen, 1988, Dean andBender, 2005).
Surf-beats are particularly dominant in storm conditions since the swell frequency band swash (nominally 0.05 < f < 0.18 Hz) is saturated, while the infragravity swash frequency band is not. Infragravity wave dissipation is relatively weak in the surf zone, while reflection from the beach face is strong (Raubenheimer and Guza, 1996). Swash models range from empirical formulations (e.g., Stockdon et al., 2006), analytical approaches (e.g., Schaeffer, 1993;Erikson et al., 2005), 1-D numerical models (e.g., , and 2-D models (e.g., van Dongeren et al., 2003;Reniers et al., 2004Reniers et al., , 2006. In this study we propose to use the 2-D wave propagation and sediment transport model XBeach (Roelvink et al, 2009) to simulate the near-shore wave propagation and dune erosion during a selected historical storm, Hurricane Irene (August, 2011), and during a NACCS 100-year storm. Irene was selected for model calibration because it is the only significant storm for which there are simultaneous wave observations at nearshore buoys in very shallow water and beach cross-section measurements. The XBeach model requires time series of waves and water elevation as boundary conditions. While this data is available for Irene as either direct data at buoys, hindcast data at Wave Information Studies (WIS, 2010) stations, or modeled data at the computational grid boundary, such time series are not available for the 100-year storm. In this case, one can select a historical storm representative of the 100-year storm, scale a historical storm to 100-year values, or create a synthetic design storm (e.g., Carley and Cox, 2003). We have chosen to select a time series from the NACCS database as a proxy 100-year storm. The storm was selected based on a comparison of maximum significant wave height, H s , and water elevation with the NACCS 100-year spectral parameters at the study area. Let us note that these 100-year conditions do not necessarily lead to the statistical 100-year beach erosion event, as this would also be dependent on other factors such as storm duration and succession of storms (Callaghan et al., 2008).

Methods 2.2.1 XBeach Model
XBeach is a 2-D nearshore numerical model developed to assess the natural coastal morphological response during time-varying storm conditions. The response includes dune erosion, overwash and breaching. The model conceptual approach is to simulate processes that are occurring in the four erosion regimes defined by Sallenger (2000): swash, collision, overwash and inundation.
The model employs a 2-D depth averaged description of the wave groups and accompanying infragravity waves to resolve the swash dynamics. The wavegroup forcing is derived from the time-varying wave-action balance (e.g., Phillips, 1977). A dissipation model  is used in combination with wave groups. A roller model (Svendsen, 1984) is used to represent momentum stored in surface "rollers", leading to a shoreward shift in wave forcing. The model has two hydrodynamic modules, a short wave and a depth-averaged flow module, as well as two morphodynamic modules, a morphology change and a sediment transport  The governing equations for each module is summarized hereafter:

Short Wave Module
Short wave propagation is based on conservation of wave action, coupled with a roller energy balance. The directional distribution of the action density is taken into account, but the frequency spectrum is parameterized and reduced to a single frequency parameter (Holthuijsen et al., 1989).
The conservation of wave action, A, defined as the ratio of the wave energy density, S w , to the intrinsic frequency, σ, is expressed as and, with D w , the total wave energy dissipation due to wave breaking; c x , c y , are the components of the wave action propagation speed propagating perpendicularly to the wave crest and at an angle θ from the underlying current u, defined as: with u L and v L , the cross-shore and alongshore depth-averaged Lagrangian velocities of the current respectively, and c g the group velocity obtained from linear wave theory.
The roller energy balance is coupled to the wave-action energy balance with the dissipation of wave energy, D w , acting as a source term for the roller energy balance: with S r , the roller energy, and D r , the total roller energy dissipation. The rollers contribute to the radiation stresses and are added to the wave standard radiation stresses calculated with linear wave theory to create the wave forcing (F x and F x ) for the mean flow equations: F x (x, y, t) = − ∂S yy,w + S yy,r ∂y + ∂S xy,w + S xy,r ∂y (2.7)

Flow Module
In the flow module, low-frequency and mean flows are modeled using depthaveraged shallow water equations. To account for the wave-induced mass flux and the subsequent return flow, these are cast into a depth-averaged Generalized Lagrangian Mean (GLM) formulation (Andrews and McIntyre, 1978;Walstra et al., 2000). Continuity and momentum equations are formulated in terms of the Lagrangian velocity (u L ,v L ) defined as the distance a water particle travels in one wave period, divided by the period.
The GLM continuity and momentum equations are therefore given by: Where η is the water level, τ bx and τ by are the bed shear stresses, ν h is the horizontal viscosity, and f is the Coriolis coefficient. The bottom shear stress terms are calculated with the Eulerian velocities experienced by the bed, u E , related to the Lagrangian velocities by subtracting the Stokes drift in the x-and y-directions respectively. The bed shear stress is calculated with: The bottom friction coefficient, c f , in this study is specified using a Manning formulation. Manning's n value is related to the dimensionless bed coefficient friction value by:

Sediment Transport Module
A depth-averaged advection-diffusion equation is used to model sediment transport (Galappatti and Vreugdenhil, 1985). The entrainment or deposition of sediment is controlled by the difference between the actual sediment concentration, C, and the equilibrium concentration, C eq . This difference is therefore the source term in the sediment transport equation: with h, the water depth, D h , the sediment diffusion coefficient, and T s , an adaptation time. T s is given by a simple approximation based on local water depth and the sediment fall velocity, w s : where a small value of T s corresponds to a nearly instantaneous sediment response.
Since it is well known that wave nonlinearity affects sediment transport (e.g., Janssen et al., 1998), and that the wave modules do not include nonlinearity, the sediment transport equation is modified to include a parametric form of wave nonlinearity based on wave skewness and asymmetry (Van Thiel de Vries, 2009).
The wave skewness and asymmetry is used to approximate the part of the sediment advection velocity enhanced by nonlinear wave effects, u a . This velocity is added to the Eulerian velocity in the advection-diffusion equation, which is therefore modified as follows: The parameter u a is calculated as a function of wave skewness (S k ), a wave asymmetry parameter (A s ), root-mean-square velocity (u rms ) and two calibration factors (f Sk , f As ): A higher value of u a simulates a stronger onshore sediment transport component. Either the individual calibration factors can be adjusted, or the option to set both calibration parameters is done by adjusting the f acua calibration parameter.

Morphology Module
The morphology module includes bed level updating and avalanching. The bed level is updated based on the sediment transport formulation as well as avalanching. The bed level change is given by: where p is the porosity, f m or is the morphological acceleration factor of order (1-10) (e.g. Reniers et al., 2004) and q x and q y represent the sediment transport rates in x-and y-directions respectively.
Avalanching is introduced to represent the slumping of the sand during the storm induced dune erosion and the bed is updated accordingly: where m cr is the critical bed slope, with a similar expression for the y-direction.
There are separate critical slopes for dry and wet points, and inundated areas are considered more prone to slumping.

Boundary Conditions
The wave energy at the offshore boundary is prescribed in XBeach as a time series of JONSWAP spectral parameters. At the lateral boundaries, Neumann boundary conditions are applied, where the longshore gradient is set to zero.
For the flow boundary conditions, the offshore boundary is set to an absorbinggenerating boundary condition (Van Dongeren and Svedsen, 1997), where outgoing waves can leave the computational domain through the boundary with minimal reflection, while also specifying incoming waves. The water level gradient is set to zero for the lateral boundaries.

Model Calibration
The impact of Hurricane Sandy was modeled using XBeach in two different case studies, on the coast of New Jersey (Nederhoff, 2014), andFire Island, New York (De Vet, 2014). In both studies, the model sensitivity to the facua parameter and the bed friction was highlighted. De Vet (2014) showed an extreme overprediction of erosion when these two parameters were kept at default values, with the model performing with much better skill after they were adjusted. As a result of these studies, these parameters are listed in the XBeach tutorial as a two-step calibration process that can be used to calibrate the model when multiple regimes of Sallenger (2000) are expected. The facua parameter is suggested to calibrate the collision regime, and increasing the bed friction on the dune is suggested to calibrate the overwash regime. The sensitivity of the model to these two parameters is tested in this study.

XBeach Model Setup
The For all simulations the grain size parameters D50 and D90 were set to 0.58 mm and 1.31 mm respectively, based on a standard sieve analysis of sandy sediment taken from the study area. All other XBeach parameters were kept at their default values unless noted otherwise in the following sections where we describe the different simulation scenarios considered.
The locations and angles relative to North of the transects shown in Figure   2.4 are specified in

Hurricane Irene
Hurricane Irene was a tropical storm by the time it impacted the RI area on August 28, 2011. GSO's cross-shore profile measurements taken before and after the storm are given in Figure 2.5. Also shown in Figure 2.5 is the XBeach model input beach topography. From the available data, comparing pre-and post-storm GSO profiles, it seems that the storm stayed within the collision regime in terms of dune erosion.
Looking at Figure 2.5, three concerns are apparent: 1) the GSO profiles do not extend into the surf zone which is a very active area of erosion and accretion.
2) The 2011 RIGIS DEM shows a different profile than that observed before Irene, which adds significant uncertainty to any comparison. 3) Irene is a storm which only had a collision impact, therefore the calibration is possible for the collision regime only.
The XBeach offshore boundary conditions for Hurricane Irene were found using the hydrodynamic model ADCIRC coupled with the wave model SWAN (Booij et al., 1999) (2017). Figure 2.6 shows the significant wave height, peak period, and water elevation time series that were used to force the XBeach model for Hurricane Irene simulations. Figure 2.6: Hurricane Irene wave height, peak period, and water elevation time series at XBeach offshore boundary, obtained from an ADCIRC-SWAN regional grid simulation As discussed above, the collision regime calls for a calibration of the f acua parameter. The f acua parameter was varied from 0.1 to 0.3 in increments of 0.05 and bed friction was held constant at the default Manning's n of 0.02. The morphological acceleration factor was set to 1. The Hurricane Irene simulation scenarios are summarized in Results are shown in the next section.

100-Year Storm
To obtain a time series of wave and water level conditions to force the XBeach model for the 100-year storm, the NACCS database was searched at two save points, one slightly further offshore than the extent of the XBeach offshore boundary (Save Point ID 9136), and one in the nearshore close to the center of the model domain (Save Point ID 868), for a storm with peak wave height and water levels close to the NACCS 100-year statistical values determined at those same save points.    In order to get a feel for the erosion impact regime expected for the 100-year storm at different locations along the dune in the study area, the equations and definitions given in Stockdon et al. (2012), based on empirical parameterization of wave setup, swash, and runup (Stockdon et al., 2006), are used and related to the storm impact model by Sallenger (2000)   η 98 , which is the extreme high water level attained during a storm, defined as the 98-percent non-exceedance level. If η 98 is larger than the dune toe but smaller than the dune crest, the collision regime is expected. If the η 98 is larger than the dune crest, the overwash regime is expected. The worst case scenario occurs when η 50 , the combination of the tide (η tide ), surge (η surge ) and wave setup (η setup ), is greater than the dune crest. Equations 2.20 through 2.24 are used to determine the relevant parameters (Stockdon, 2012).
With S, the total swash excursion about the setup level, H 0 , the deepwater wave height, L 0 , the deep water wave period, T , the wave period, g, the gravitational acceleration (9.81 m/s 2 ), and β m , the beach slope. The beach slope is defined as the slope of the line from the mean high water mark to the dune toe. Note that the equation for η 98 is slightly different from what is shown in Figure 2.8, with the factor 1.1 used to account for parameterization bias (Stockdon et al., 2012).
This analysis was performed at the seven cross-shore transects, previously introduced: the 5 GSO transect locations, and the 2 FEMA transect locations (        however, the simulations completely erode the dune at transect 19 to a level below zero NAVD88, and at transect 20 just slightly above zero NAVD88. The Oakley profile assumes a maximum elevation of about 1.6 m NAVD88, and is located more or less where the dune toe used to be.
A clear advantage of using the XBeach model is that the resultant dune profiles are known everywhere, rather than just at the site of specific transects. The Oakley and FEMA profiles do not indicate areas of breaching such as the XBeach model shows.

Conclusion
The 100-year storm erosion was simulated using the morphodynamic model XBeach. The NACCS database was used to find a synthetic storm time series at the offshore XBeach boundary to force the model. The storm chosen had peak wave height and water levels similar to the NACCS 100-year statistical values for each of those parameters at the study area location. The results are presented in map form as well as at cross-shore transects within the study area. The model shows sensitivity to the two parameters suggested by Nederhoff (2014) for XBeach calibration, the facua parameter and the bottom friction.

Results of modeling the historical storm Hurricane Irene demonstrate
XBeach's sensitivity to the facua parameter. In view of the uncertainty associated with the limited size of the observed transects, an exact calibration of this parameter was not possible and it was decided to approach the 100-year storm simulations from a sensitivity perspective, providing a range of possible dune profiles associated with a range of potential parameter values. As the erosion process was limited to the collision regime during Hurricane Irene, the calibration was therefore limited to this regime; the 100-year storm is expected to be in the overwash and inundation regimes and no data was readily available for bottom friction calibration.
Simulations of the 100-year storm for a range of realistic parameters show a range of results from a fully eroded dune to a dune progressing landward with a crest elevated above 2 m in some locations. Coastal geologist Oakley (2015) shows that, based on historical storms, the 100-year profile would provide a flattened dune with a crest of about 1.6 m NAVD88 located at about where the original dune toe was. The simulation with variable bottom friction results in the closest predicted elevation to this scenario, although there are significant differences where the dune is located as well as dune shape, when comparing the simulation results to the generalized barrier profile.
The high sensitivity of the model to both the facua and bed friction parameters demonstrates the necessity to have reliable field measurements covering not just the dune profile but also the full equilibrium beach profile, ideally extending underwater up to 300m offshore. The recent practice of collecting pre-and post-storm LiDAR data when extreme storm events are expected will aid in the calibration of such a model so that it can be used to make predictions with more certainty.