Storm Surge and Wave Impact of Low-Probability Hurricanes on the Lower Delaware Bay—Calibration and Application

Hurricanes pose major threats to coastal communities and sensitive infrastructure, including nuclear power plants, located in the vicinity of hurricane-prone coastal regions. This study focuses on evaluating the storm surge and wave impact of low-probability hurricanes on the lower Delaware Bay using the Delft3D dynamically coupled wave and flow model. The model comprised Overall and Nested domains. The Overall model domain encompassed portions of the Atlantic Ocean, Delaware Bay, and Chesapeake Bay. The two-level Nested model domains encompassed the Delaware Estuary, its floodplain, and a portion of the continental shelf. Low-probability hurricanes are critical considerations in designing and licensing of new nuclear power plants as well as in establishing mitigating strategies for existing power facilities and other infrastructure types. The philosophy behind low-probability hurricane modeling is to establish reasonable water surface elevation and wave characteristics that have very low to no probability of being exceeded in the region. The area of interest (AOI) is located on the west bank of Delaware Bay, almost 16 miles upstream of its mouth. The model was first calibrated for Hurricane Isabel (2003) and then applied to synthetic hurricanes with very low probability of occurrence to establish the storm surge envelope at the AOI. The model calibration results agreed reasonably well with field observations of water surface elevation, wind velocity, wave height, and wave period. A range of meteorological, storm track direction, and storm bearing parameters that produce the highest sustained wind speeds were estimated using the National Weather Service (NWS) methodology and applied to the model. Simulations resulted in a maximum stillwater elevation and wave height of 7.5 m NAVD88 and 2.5 m, respectively, at the AOI. Comparison of results with the U.S. Army Corps of Engineers, North Atlantic Coastal Comprehensive Study (USACE-NACCS) storm surge values at the AOI demonstrates that the estimated elevation has an annual exceedance probability of less than 10 − 4 .


Introduction
Tropical cyclones, hurricanes, and typhoons can be categorized among the most dramatic and fatal natural phenomena. Hurricanes have been the costliest and deadliest natural hazards in U.S. history [1]. The damages to the Washington D.C. area during Hurricane Isabel (2003) were estimated to cost 3.3 billion dollars [2,3]. Similarly, the National Hurricane Center (NHC) of the National Oceanic & Atmospheric Administration (NOAA), reported an economic cost of 27 billion dollars to New York City and New Jersey as a result of Hurricane Sandy (2012). More than 650,000 homes were destroyed and 8 million people lost power in these two densely populated metropolitan areas. The combination of those events and of the 11 March 2011 accident at the Fukushima Daiichi nuclear power plant resulting from the Great Tohoku Earthquake and subsequent tsunami necessitated the reevaluation of A qualitative review of the available historical hurricanes and their impact might imply that hurricane storm surge and wave with very low probability is of no concern; however, recent research, such as Dailey et al. and Donnelly et al. [12][13][14], suggests that analysis using only historical hurricane records is, at best, inconclusive. Delaware Bay and its floodplain received less attention than their surrounding estuaries and coastal regions, especially with the impact of low-probability hurricanes. In this study, care is taken to filling this gap and addressing the practical aspects of storm surge and wave modeling as it relates to vulnerable sensitive infrastructures, low-probability events, optimized parameter selection, wind-field generation, calibration, and computational needs. The goal was to develop a numerical model with the least computational requirement and capable of being used for facilities located near coastal regions.

Deterministic vs. Probabilistic
Two broad approaches in performing studies associated with flood hazards are deterministic and probabilistic. In the deterministic approach, the probable maximum storm surge is determined by modeling the synthetic hurricane tracks generated using the National Weather Service (NWS) Report 23 [15] guideline. NWS Report 23 provided the criteria for determining wind field for the most severe hurricane with highest sustained wind speed at a specified coastal location, yet reasonably probable. The wind field is then run in a hydrodynamic model, such as ADvanced CIRCulation (ADCIRC), Delft3D, TELEMAC, etc., to determine the probable maximum water surface elevation (i.e., low probability storm surge). Meteorological parameters evaluated in NWS Report 23 are central pressure (P c ), peripheral pressure (P a ), radius of maximum winds (R max ), forward speed (V f ), track direction (θ), and inflow angle (φ).
Probabilistic flood hazards are used to estimate the frequency of external flooding scenarios. For sensitive infrastructures, including nuclear power plants, the extrapolation of the hazard to extremely low frequency (less than 10 −4 ) is required. There are various frequency estimation methods, such as gage analysis, Monte Carlo, Empirical Simulation Technique (EST), Joint Probability Method (JPM), and JPM-OS (Optimal Sampling) [16,17]. EST is a bootstrapping, resampling-based method that assumes similar magnitude and frequency between past and future events. The JPM and JPM-OS approaches are used more commonly today. JPM utilizes a parametric representation of tropical cyclones characterizing intensity, size, track, and speed [18]. It then assigns a joint probability to each combination of storm parameters. JPM relies on a numerical method to calculate coastal flood elevations that would be generated by those storms. One advantage of the probabilistic versus deterministic approach is the ability to assign a frequency to the event. In this study, the deterministic approach is followed to establish the parameters required for developing a very low-probability wind field. The probability of estimated surge values is found indirectly by comparing them to the probability of other studies.

Study Area
The area of interest (AOI) is located on the west bank of Delaware Bay (Bay) at approximately 32 km north of its confluence with the Atlantic Ocean and 12 km northeast of Milford, Delaware.
The Bay has a length of 210 km starting from its riverine portion at Trenton, New Jersey to its connection with the Atlantic Ocean at the bay side [19]. The Bay width gradually increases from upstream to downstream. The mouth of the Bay extends between Cape May and Cape Henlopen with a width of 18 km. The estuary reaches its maximum width of 45 km (see Figure 1) approximately 11 km northwest of the mouth. The Bay has an average bathymetry of 7 m [20] with a channel (deeper than 28 m) passing in the middle of the Bay. The channel acts as a conduit for connecting the Atlantic Ocean to the ports and facilities upstream. The estuary has an axis that is tilted 35 degrees counterclockwise west of due north.  The AOI is within a wetlands area with a saline fringe geographic setting. Reed et al. [21] concluded that with approximately 2 mm/year of sea level rise, the wetlands around the AOI would be marginal and transformed into open water. This is very important as the wetlands tend to act as a buffer in attenuating the storm surge [22][23][24]. Shorelines of the AOI are undisturbed sand beaches.
Since 1749, 108 tropical cyclones, including some that became extratropical, have affected the area [25]. Most have occurred in September, which coincides with the peak of the Atlantic hurricane season. The majority of the tropical cyclones that have occurred in Delaware have brought only rainfall or strong waves (higher than 1.5 m), though a few have caused deaths in the state. Six of the most deadly hurricanes with their general track are graphed in Figure 1. Hurricane Isabel was selected for calibration purposes as is described in Section 2.4.5. For Hurricane Isabel, the hourly position of the eye is shown in yellow circles.

Computational Models
Delft3D (Version 4, Deltares, Delft, The Netherlands) is a highly developed computational model with several modules and utilities that are fully integrated for solving multi-disciplinary problems in coastal, river, and estuarine environments [26]. Three Delft3D components, FLOW, WAVE and WES, are dynamically coupled in this study for the automatic transfer and exchange of relevant data. Major features of each module relevant to the hurricane modeling are described in following subsections.

Delft3D-FLOW Module
The FLOW module is the heart of Delft3D and a multidimensional hydrodynamic and transport simulation program that calculates non-steady flow and transport phenomena resulting from riverine, tidal, and meteorological (i.e., atmospheric, winds, and waves) forcing [27].
The Delft3D-Flow module in two-dimensional (2D), depth-averaged, hydrostatic, and unsteady form is used for characterizing the hydrodynamics of the historical and synthetic hurricanes. In this format, the continuity [Equation (1)] and Reynolds-averaged Navier-Stokes [Equations (2) and (3)] equations for an incompressible fluid are solved over an orthogonal curvilinear grid in spherical coordinates: ∂ζ ∂t ∂v ∂t where u and v are depth-averaged velocities along the x and y directions, ζ is the water surface elevation above the still water level, f is the Coriolis parameter, g is the gravitational acceleration, ρ is the water density, d is the water depth above the bottom of and below the still water level, ν H is the horizontal viscosity coefficient, U 10 and V 10 are wind velocities at 10-m height above the still water level along the x and y directions, n is the Manning coefficient, and C d is the wind-drag coefficient dependent on U 10 and V 10 . Manning's "n" and wind-drag values (C d ) are two important calibration coefficients in all coastal flood modeling efforts. Manning's "n" is a factor of land cover and changes from 0.02 in open water to 0.15 in high-density and developed areas [28]. Manning's roughness has been used in similar hurricane storm surge studies. Garzon and Ferreira [3] study results suggested that the Manning's "n" coefficient can significantly impact the model accuracy in estimating peak of maximum elevation, currents, flood duration, and extension in overland areas.
Wind-drag translates the wind speed into wind setup along the coastline [27] and growth of surface waves [29]. Wind-drag coefficient increases with increasing wind speed and then is capped at a higher breakpoint wind speed (e.g., 30 m/s based on Vatvani et al. [30]). This reflects the dependence of surface roughness on wind speed. Wind-drag functionality has been studied and modified by various research, including [31][32][33][34][35]. They showed that C d ranges from 0.0015 to 0.005. Delft3D-FLOW module introduces the three breakpoint piece-wise formulation [27] of the wind speed and drag relation. Such functionality provides flexibility in defining several relations.

Delft3D-WAVE Module
The Delft3D-WAVE module is the Simulating WAves Nearshore (SWAN) model. SWAN is a fully spectral third-generation wave model that can realistically simulate generation, propagation, transformation, breaking, and associated hydrodynamic forcing of random, short-crested, wind-generated waves [36]. The SWAN model is properly formulated to consider most important physical aspects of wave propagation, generation, and dissipation [37]. Coupling FLOW and WAVE modules is essential in storm-surge modeling as water levels, currents, and waves are interconnected and affect each other through radiation stress gradients. SWAN has traditionally been used for nearshore and shallow waters. Dietrich et al. [38] provided a comprehensive examination of applying SWAN to deep waters with hurricane modeling applications. The SWAN model has been used extensively to simulate waves in both shallow and deep waters [3,30,36,39,40]. SWAN calculates the evolution, in 2D geographical x-space and time t, of the wave action density N using the action balance equation: The left side of Equation (4) is the kinematic part, which represents, respectively, the temporal variation of the action density, the propagation of action in geographical space (with c g representing group velocity and U the ambient current), shifting of the relative frequency due to variation in depth and current (with c σ as propagation velocity in σ-space), depth-induced and current-induced refraction (with c θ as propagation velocity in θ-space). The term S tot on the right-hand side is the source term in terms of energy density representing the effects of generation, dissipation, and nonlinear wave-wave interactions [36].
Holland [45] introduced a scaled hyperbolic pressure profile used in gradient wind equations. The main hurricane-related parameters required in this model are radius of maximum wind (R max ), maximum wind speed V max , central pressure (p c ), peripheral pressure (p a ), and track angle (θ). The Holland model has been used in numerous studies, including [47] and [48].
The NWS Report 23 vortex model [15] analytically recreates the spatially distributed hurricane wind and pressure fields. The Rankine model [49] assumes solid body rotation inside the R max with decaying wind speed beyond the R max . TC96 is a highly-refined mesoscale-moving vortex formulation based on the equation of horizontal motion, vertically averaged through the depth of the planetary boundary layer. It was first developed by Chow [50] and then modified by Cardone et al. [51]. The WRAMS model is an implementation of the Regional Atmospheric Modeling System (RAMS) by WeatherFlow. RAMS [52,53] is a mesoscale, nonhydrostatic, 3D modeling system. It has recently been used to evaluate the impact of the areas affected by Hurricane Sandy [9]. MM5 [46] is a nonhydrostatic scheme that simulates the atmospheric circulation.
WES is a modification of the Holland scheme to introduce asymmetry features of a real hurricane. This asymmetry is brought about by applying the translation speed of the cyclone center displacement as steering current and by introducing rotation of wind speed due to friction. Delft3D uses the WES technique in generating the wind field, which is the applied model in this study. In Delft3D, wind fields are generated on a moving spider web grid using Equation (5). Spider web is a grid in polar coordinate that is centered on the center of the tropical cyclone. At the center, this circular grid system has a high-grid resolution, which enables capturing strong variation of wind and pressure gradients in the center area of the typhoon. The grid resolution decreases proportionally with increase of the distance from the center of the grid: where r is the distance from the center of the cyclone, f is Coriolis parameter, ρ a is density of air (i.e., 1.10 kg m −3 ), p drop is difference between peripheral and central pressure, V max is the maximum sustained wind speed, R max is the radius of maximum wind, and e is the base of natural logarithm (i.e., 2.71828). Laknath et al. [54] compared the wind field generated by the MM5 model and by the WES techniques. They reported that the WES technique underestimates wind for areas located far from the center of the storm. However, in general, a hurricane contribution would be minimal to those remote areas and, therefore, is not considered a flaw of the WES scheme.

Model Domain
The physics-based numerical models, such as Delft3D, are valuable for accurately predicting the physical processes involved in hurricane surge response over large regions; however, they are computationally intensive when used at the spatial resolution required for accurate surge estimation [55].
Spectral, temporal, and spatial resolution are major decision points in each modeling mechanisms as they impose limitation on circulation and wave analysis. Nesting of structured meshes is a well-known technique in improving the negative impact of resolution. In this technique, the size of the model mesh decreases in a progressive fashion around the AOI. It is established that the storm surge models are sensitive to the open-boundary conditions, spatial resolution, and the size of model domain [38]. However, any model should also be economically justifiable and efficient with limiting the model error within an acceptable range or in making conservative assumptions.
Li et al. [56] investigated the impact of the model domain size on storm surge. They identified a threshold where increasing domain size, both cross-shore and alongshore, has no further impact on simulation results. They also showed that the domain size is not impacted by hurricane intensity, but has linear relation with the hurricane radius of maximum wind (R max ). This finding enables a modeler to establish a model domain that is calibrated for a significant hurricane and then apply it for a lower-probability event.
Another consideration in defining a cost-effective and efficient model domain is the location of the offshore boundary. A study by Garzon and Ferreira [3] suggested that the location of the open boundary might not impact the simulation results as long as the model resolution is sufficient to properly model the movement of water. However, Blain et al. [57] showed that locating the offshore boundary of the model in an area with significant surge, causes the near-shore surge to be underestimated. This indirectly indicates that the offshore boundary should be located beyond the continental shelf. Li et al. [56] reported that the slope of the Atlantic continental shelf ranges from 0.0003 to 0.01. Their study, combined with Irish et al. [58], showed that the storm surge increases with decrease of the slope of the continental shelf. They concluded that, for a shallow bathymetry, the domain size should be increased to avoid artificial impact of the slope. The width of the shelf near our study area is approximately 154 km. Based on Li et al. [56], storm surge is not sensitive to the change in shelf width for width higher than 100 km. According to their analysis, the appropriate domain size for this study, considering the width of the shelf, is 100 km (cross-shore) by 150 km (along-shore). In addition, considering the hurricane categories (higher than 5), the minimum domain size is 200 km (cross-shore) by 250 km (along-shore).
The above findings and preliminary sensitivity analysis resulted in one Overall and two Nested domains (see Figure 2) in this study. All domains are set up in a 2D fashion. RGFGRID tool of Delft3D is used to generate the curvilinear grid. Primary grid properties including curvature, smoothness, length of edge, and orthogonality are used to increase grid quality. Figures 3 and 4, respectively, show the bathymetry of the Overall and Nested domains. Coastline and bathymetry are from the ADvanced CIRCulation (ADCIRC) Federal Emergency Management Agency (FEMA) Region III model [59]. The bathymetry is developed by the Coastal Process Branch of the Flood and Storm Protection Division, U.S. Army Engineering Research and Development Center (ERDC) using best available data. The QUICKIN tool of Delft3D combined with ArcGIS from Esri are used to assign bathymetry data to model grids. In addition, high-resolution Light Detection and Ranging (LiDAR) data were used to complement the topography near the study area.  All model domains are oriented along the axis of the Delaware Bay, which is southeast towards northwest. Considering their size, all model domains are in a spherical coordinate system. Other major properties associated with each model domain are presented in following subsections.

Overall Domain
The Overall domain includes the Delaware Bay extending from the mouth to the Commodore Barry Bridge on the Delaware River ( Figure 2). The offshore portion of the domain extends approximately 2100 km toward the Atlantic Ocean and beyond the continental shelf. The model domain has a total length of 2700 km and width of 2400 km. The grid contains 87,963 cells, covering an area of 5.1 million km 2 . The Overall model domain has grids with an averaged constant resolution of 7.5 km (i.e., 0.06 • ). The Overall domain is considered large enough to ensure that the wind field gradient and offshore winds are properly captured. This resolution was determined to be optimal in resolving offshore surge and wave characteristics. Sensitivity analysis was performed with coarser and finer resolution models, but resulted in similar results. For example, increasing the number of grid points from 87,963 to 160,000, almost twice, improved the model results by only 3 cm. However, the model run time was increased by a factor of 5.  . These resolutions were determined to be optimal in resolving surge and wave characteristics. The variation in water level of less than 50 cm was considered as the criteria. Sensitivity analysis was performed with coarser and finer resolution domains, but resulted in similar results.

Standard Setting
Delft3D-Flow uses the standard Alternating Direction Implicit (ADI) for solving the shallow water equations described above. Sensitivity tests are carried out in order to determine the largest time step for which the ADI-method still yields accurate results. Temporal differencing for the simulation is set at three minutes for the Overall grid and one minute for the Nested grids to avoid numerical instabilities. Delft3D-Wave (i.e., SWAN) is based on implicit numerical schemes, and, therefore, it is not limited by the Courant stability criterion as a balance between time and space. In this sense, the time step in SWAN is not restricted. However, the accuracy of the SWAN results is obviously affected by the time step. The Sheng et al. [60] study showed that, for slow-moving hurricanes, the stationary mode of SWAN yields reasonable wave conditions as compared to the non-stationary mode. In simulating storm surge, coupling between flow and wave models as well as coupling between Overall and Nested grids is common. The flow and wave simulations are performed in one-way coupling mode for the Overall domain, meaning that only currents and water levels from the flow model are fed into the wave model. After establishing the open-boundary conditions for the Nested1 domain, surge and wave are coupled in two-way fashion, where results calculated by one (i.e., water levels and current by Delft3D-FLOW and radiation stresses gradient by Delft3D-Wave) are fed back and forced to the other within a specified time steps. All domains always use the same wind field and drag relations.
Flow and wave models are run on the same grids to produce a seamless coupling. SWAN is a slow-running model and is typically run only once for a few flow model time steps. Sensitivity analysis and literature review showed that a 20-min coupling provides satisfactory results. This means that, after twenty 60-s Delft3D-Flow time steps, flow and wave models mutually exchange information. In this fashion, Delft3D-Flow receives radiation stresses; in return, Delft3D-Wave updates water level and current.
The main flow model parameters include initial water level of 0.0 m, fluid density of 1025 kg m −3 , and constant horizontal eddy viscosity of 10 m 2 s −1 . Drying and flooding are included in all models with the criterion set to 10 cm.
The spectral resolution in SWAN was configured with 36 frequency bins over a range of 0.04-0.6 Hz and directional space bins of 10 • . The frequency range is found through extensive literature review and sensitivity analysis. Results showed that the lowest frequency must be somewhat smaller than 0.7 times the value of the lowest expected frequency. The highest frequency must be at least 2.5 to 3 times the highest expected frequency. For example, hurricane Isabel frequency varied from 0.06 to 0.2 Hz.
The 2nd generation of SWAN, with bottom friction of 0.067 m 2 s −3 parameterized in the JONSWAP formulation, is used. Sensitivity analysis on using both 3rd and 2nd generation schemes were performed. The 2nd generation resulted in improvement in the wave results and reduced the model run time. Findings are compatible with other researchers such as [61]. In order to set the bottom friction in JONSWAP model, a range of frictions from 0.038 m 2 s −3 to 0.067 m 2 s −3 , as suggested by Vledder et al. [62], were tested. The depth-induced process of wave breaking is not entirely understood yet in spectral mode. In contrast, the dissipation is simulated using the bore model that conserves the spectral shape [63]. In this study, a disputation coefficient of α = 1.00 and a breaking parameter of γ = 0.73 is applied. Whitecapping, nonlinear triad, refraction, and frequency shift were not considered. Sensitivity analysis showed inclusion of these formulations has minimal impact on the ultimate result, which is the maximum probable storm surge at the AOI.
The default drag coefficient (C d ) of 0.0012 underestimated the generated wave heights at the buoy locations. Therefore, after sensitivity analysis, the drag coefficient was set to vary linearly from 0.0012 at 0 m/s wind speed to 0.004 at 30 m/s wind speed. Such modification has been reported in other recent studies, such as [48,64].

Model Forcing
Tides and winds are the primary forces considered in this study as described in following subsections.

River Discharges
Kerr et al. and Resio et al. [44,65] studies showed that river discharge can significantly impact storm surge along confined rivers; however, they also stated that the river discharge can be neglected. In this study, the AOI is located in the bay side of the Delaware Estuary, which is significantly wider than the upstream portion (i.e., beyond the confluence of Chesapeake and Delaware Bay Canal and Delaware Bay). Simulating storm surges using a range of river flows applied to Trenton, New Jersey resulted in negligible difference in water levels in the study area. Therefore, no river discharge is considered in this analysis.

Tides
As shown in Figure 2 In the equation above, H(t) represents temporal variation of water level, k the number of constituents, i the constituent index, A i the local amplitude of tide, φ i the constituent phase, F i the nodal amplitude factor, T i the period, and V i the astronomical argument. The A i and F i are input variables for each selected constituent. Delft3D calculates T i , F i , and V i for each constituent.
In Delft3D, the phase and amplitude are required on both sides of each open boundary. ADCIRC is used to create these parameters from the western North Atlantic, Caribbean Ocean and Gulf of Mexico tidal databases [67]. Tides were applied only to the calibration portion of the study and not the performance runs. In order to estimate the probable maximum storm surge, 10% exceedance high tide (EHT) is used as initial condition in all performance runs. Based on ANSI standard [68], this has a value of 1.7 m above Mean Low Water (MLW) near the study area at the Breakwater Harbor, Delaware. Converting the MLW to NAVD88 (same as bathymetry data) makes the 10% EHT equal to 0.95 m NAVD88.

Meteorological Forcing-Wind Field
Precise prediction of the maximum water level requires the accurate computation of the storm surge, which, again, depends on the accuracy of the wind forcing it receives. Wind fields were developed for both historical and synthetic storms for input to the Delft3D-Flow and Wave modules using the methodology described in Section 2.2.3, with the following inputs.

Delft3D-WES Setting
As mentioned above, the Delft3D-WES module computes the hurricane wind and pressure fields on a so-called spiderweb, which provides a moving snapshot of the hurricane wind field. Several parameters, including radius of the tropical cyclone (2000 km), number of rows (60), and number of columns (1200), define the size of the spiderweb grid. These parameters were carefully selected through a series of sensitivity analyses and calibration. The resolution in the radial direction then becomes equal to 1666 m. A resolution between 1000 and 2000 m is required especially when analyzing compact cyclones with high wind gradients [69]. Sensitivity analysis of Hurricane Isabel showed that increasing cyclone radius results in increased generated winds until a threshold, where its increase has no more impact. Increasing the number of rows and columns in the spiderweb did not have a significant impact on the results.
Considering the available input parameters, seven methods are implemented inside the WES module. Method 3, requiring location of cyclone eye (i.e., latitude and longitude), maximum sustained wind speed (V max ), pressures drop (P drop , i.e., difference between peripheral pressure and the minimum central pressure), and the radius of maximum winds (R max ) was used. A peripheral pressure of 1013 mb was applied.

Historical Storm
Based on U.S. Nuclear Regulatory Commission (NRC) guidelines, models being used for forecasting low-probability events should be validated using the most significant and recent event relevant to the region. In addition, considering the counterclockwise rotation of hurricanes, only those with a track line located west of the study area have the capability of producing sustained wind speed along the estuary. The major historical hurricanes that impacted the region were extracted from NOAA database [70]. The combination of National Hurricane Center (NHC), also known as HURDAT [72], and extended best-track datasets [73] were used for establishing the model inputs. Data analysis showed that Hurricane Isabel can be divided into seven segments with similar characteristics (see Table 1). The hourly position of the hurricane, along with their characteristics, was then interpolated from these seven segments. The Hurricane Isabel track with hourly position of the eye is presented in Figure 1. Velocities shown in Table 1 are the 1-min maximum sustained surface winds. Probabilistic storm surge studies require large sets of storm tracks, on the order of thousands, to adequately cover the parameter and probability spaces. For example, more than 17,000 storm tracks were used in Hanson et al. [10], more than 8000 tracks in Lin et al. [2], and more than 1050 synthetic tropical cyclones in Nadal et al. [74]. Lack of sufficient historical records makes the development of synthetic hurricane tracks a critical step in estimating hazards from coastal flooding [17]. Such tracks should follow similar statistical properties as recorded historical hurricanes within the study area [6,75,76].
Numerical surge and wave models are then used to translate stochastically generated hurricane tracks into coastal hazard estimates. However, a majority of these studies rely on fast-running models, such as Sea, Lake, and Overland Surges (SLOSH) developed by NOAA [77], in order to perform screening and, thereby, reducing the number of tracks required on the order of hundreds. However, SLOSH simulation may not be able to capture some unusual water responses to storms at locations with complex geophysical features. Other computationally intensive hydrodynamic models, such as ADCIRC, Delft3D, or TELEMAC3D, may be used to evaluate tracks selected by SLOSH. The probability of each storm is then combined with the modeling surge values to create a cumulative distribution function (CDF) to estimate the flood risk.
Researchers have shown that only hurricanes passing through a narrow band of landfall locations that generate maximum wind speed near the study site area will result in the largest surges at that specific location [44,55]. Furthermore, in the Northern Hemisphere, the Coriolis effect causes the peak of the storm surge to fall to the right of the storm track at landfall [78,79] In this study, the range of meteorological hurricane parameters is estimated combining the evaluation of site-specific features and the NWS Report 23 [15]. This guideline provides criteria for determining wind fields along the Gulf of Mexico and East Coast of the U.S., for the most severe hurricane reasonably possible at the study area known as Probable Maximum Hurricane (PMH). NWS Report 23 provides a single value for the peripheral and central pressure, along with upper and lower limits for other hurricane parameters, such as radius of maximum winds, forward speed, track direction, and maximum sustained wind speed. The procedure starts with calculating the distance between the study area and the U.S.-Mexico border in nautical miles (nm). The AOI has a NWS mile post of 4074.4 km (2200 nm). Other parameters are found based on this location index. The parameters estimated using NWS Report 23 are adjusted after evaluating the major historical hurricanes within the study area, ensuring that the most severe hurricane range is enveloped. Historical hurricane parameters are accessible from [80,81]. The range of the hurricane parameters extracted from major hurricanes within the region is presented in Table 2. After performing sensitivity analysis on numerous hurricane tracks, 41 categorized in 25 common landfalls and angles were used for the analysis (see Figure 2). Table 3 presents the range of parameters used to construct the synthetic hurricanes. Another important parameter in estimating the maximum sustained wind speed is the Holland B. Ref. [82] shows that the Holland B parameter is a factor of sea surface temperature (SST) and varies from 0.75 to 2. For the Delaware Bay, the SST during the hurricane season (August-October), varies from 20 to 25 degrees Celsius. Referring to [82], the average Holland parameter becomes one for the region. Sensitivity analysis on forward speed showed that it only changes the timing of the peak and not its magnitude; therefore, forward speed is used only to control the hurricane approach toward the AOI.

Model Calibration and Sensitivity
Several studies have been performed to understand the impact of various model parameters on storm surge estimation [56,57,83].
Bastidas et al. [17] performed a comprehensive sensitivity analysis for storm surge and wave modeling using Delft3D within the hurricane hazard framework. Findings of such studies are appropriate; however, it is noted that parameter sensitivity suited for one hurricane or one region may not be appropriate for others. Therefore, in addition to considering the findings of previous research, further sensitivity analyses were performed as part of this modeling effort. In calibrating the models, it should be emphasized that the sole goal was to determine the maximum values and not the timing, as the intent would be for power plant designs. Therefore, only the peak values were used for calibrating and not the timing aspects. Further calibration would be required if used for flood timing assessment.

Observations
Hourly measurements of wind speed from buoy 44009 [84] were used to evaluate Delft3D-WES reproduction of the wind speed time series. Peak significant wave height and period during Hurricane Isabel were extracted from buoy 44009. This station is located 48 km southeast of Cape May, New Jersey, at a water depth of 43 m, which is the closest to the study area. They were used to evaluate Delft3D-Wave reproduction of significant wave height and period.
Hourly water level records from several NOAA tidal stations [85], including Atlantic City, Ocean City, Cape May, Lewes, and Brandywine Shoal Light, were evaluated as part of the calibration and are presented in Figure 5. The figure shows that the Brandywine Shoal Light tide gauge has the highest range of water level and is the one closest to the study area. Therefore, it was used to evaluate the Delft3D-Flow ability to reproduce water level elevations during Hurricane Isabel.

Calibration Model Comparison
Delft3D-FLOW and WAVE models are both calibrated. The calibration procedure is further discussed in detail in this section. A 15-day ramping period, including only tidal forcing, was completed prior to the validation run to ensure that water levels were correctly represented at the start of the simulations. The model was then run for eight days, from 12 September 2003, to 20 September 2003, with forcing as discussed above. In this period, Hurricane Isabel propagated from the southern boundary of the Overall model domain, followed a northwest track, made landfall on 18 September 2003, and dissipated over land. Wind field, water level, significant wave height, and wave period were all compared with the observed data at their associated stations. Model performance was tested by calculating root-mean-squared error (RMSE), bias, and scatter index (SI) using the following: where O i are observations, S i are model predictions and N is the total number of points. Comparison of simulated and observed wind at NDBC station 44009 is presented in Figure 6. Model results show that the simulated peak wind speed and trend are matching the observed values. The wind field is over-predicted with RMSE of 1.1 m/s, particularly when the peak of the storm reaches the targeted buoy, near the mouth of the Delaware Bay. Figure 7 shows  Several model simulations were undertaken to calibrate the Delft3D-WAVE model, with each run varying specific parameter or friction factor. Calibration was conducted against measured significant wave height and peak wave period at NDBC Buoy 44009 described earlier. The final parameters set after this process were described in Section 2.3.4. Wave calibration results are presented in Table 4. Results are provided for all three levels of the model domain to show the improvement brought about with implementation of the nesting technique. In general, the wave model tends to overestimate the larger wave heights (i.e., waves beyond the Continental shelf) during the peak of the hurricane. Results best matched for the Nested2 grid, with an RMSE of 0.27 m and 1.0 s for significant wave height and peak period, respectively. Similar to previous studies, such as Dietrich et al. [38], the significant wave height timing deviation between simulated and recorded values was observed. The model captures well the peak significant wave height and storm surge, with only slight overestimation of 0.04 m and 0.3 m, respectively. Overall, given the magnitude of errors, the model performance was accepted, especially around the study area.  As described earlier, observed waves are not available inside the bay. Therefore, in order to ensure proper calibration of the Delft3D-WAVE model inside the bay and close to the site, results are compared with their associated values generated using FEMA Region III [10]. Comparison is made for Hurricane Isabel at a location near the NOAA tidal gage at Brandywine Shoal Light, DE.
The peak significant wave height by FEMA Region III and Delft3D-WAVE models are 2.3 m and 2.37 m, respectively. The difference is only 7 cm, which is an indication of great model performance.
Some models use a cluster of hundreds of computers to handle high-resolution computations over a large domain. The usage of cluster is avoided in this study; however, high-resolution computations are provided near the area of interest (AOI). The simulations for this modeling effort were performed on a Dell OptiPlex 7010 PC Workstation with Windows 7 Enterprise (64-bit edition) and an Intel Xeon E5-2690 Processor (2.99 GHz) with 16 GB RAM. Combined storm surge and wave model run for 1 h on this computer to simulate each three days of real time.

Results and Discussion
Upon satisfactory calibration of the model for wind, wave, and water level, the model is used to further simulate synthetic hurricane tracks (Section 2.4.6) and estimate the probable maximum hurricane. Of all the synthetic hurricane tracks, the one with a central pressure of 920 mb, radius of maximum wind of 45 nm, and a maximum sustained wind speed of 108 knots resulted in the highest water level of 7.54 m NAVD88 at AOI. The track, as highlighted in Figure 2, is oriented 60 • counterclockwise due north. It makes landfall approximately 45 nm west of the Delaware Bay mouth near Thurf Marsh Islands. The left panel of Figure 10 shows the colored pattern of water surface elevation, in meters NAVD88; vectors show depth-averaged current velocities in m/s (length of the vector is proportional to current speed). The right panel shows a colored pattern of wind speed in m/s; wind vectors are also shown.
The primary force of generating surge and surface waves is the wind [86]. In this hurricane track, a constant wind blows from the southeast, east, and south, across the Continental Shelf for several days. These prevailing directions are due to counter-clockwise rotation of the storm around the center as it moves northwest through the Atlantic Ocean. Wind speeds increased exponentially from 10 m/s 16 h before the storm's first landfall to 37.6 m/s 2 h after landfall, where the maximum wind speed and surge occurs, and then subsides after landfall, with approximate image mirror of the upward slope. The water surface elevation at the AOI rapidly increases to 2 m NAVD88 8 h before the storm's first landfall, when the center of the hurricane is 300 km offshore. It then grows exponentially to the maximum of 7.54 m NAVD88 2 h after landfall. Site grade elevation is 1.5 m NAVD88, and, therefore, there will be approximately 6 m of water at the site during the peak of storm. The AOI stays inundated for 20 h. The peaks of the storm surge and wind speed are synchronized; however, there is a 9-h lag between the growth of wind speed and water level. A comparison between storm surges generated by Hurricane Isabel and Synthetic Track 12 shows that the synthetic hurricane resulted in approximately six times higher surge than the historical event of Isabel.
Storm surge regionally built up to levels of 7-10 m NAVD88 in and around the AOI (see Figure 10). Computed storm surge was rather uniform along the axis of the lower Delaware Bay. Current vectors show water being driven into the region from the southeast by wind. This southeast-to-northwest path was the predominant pattern prior to landfall. Storm surge that had built up near the AOI propagated 18 km beyond the shoreline. As the storm continued to move northwest, surge surrounding the AOI continued to decrease along its entire periphery.
The pattern of the significant wave height (Hs) for the Nested2 model at the time of maximum wave along with associated time series is presented in Figure 11. The synthetic hurricane that resulted in maximum storm surge is used. Hs started to grow at the AOI 10 h before landfall and increased exponentially up to a maximum of 2.54 m 2 h after landfall. The associated maximum peak wave period was 7.4 s. Hs varied from 2 m to 4 m around the site with larger waves in the range of 10 m at about 48 km offshore. The peak of Hs, storm surge, and wind speed are synchronized; however, there is an 8-h lag between the growth of wind speed and significant wave height. The combined Overall, Nested1, and Nested2 Delft3D models' application provided accurate calculations of, and extremely valuable insights into, the development and propagation of hurricane-induced storm surge into the region.
Comparing the impact of various synthetic hurricanes on the water surface elevation (WSEL) at the AOI (see Table 5), the following results are yielded: 1. Of the tracks with the same landfall location, central pressure (c p ), radius of maximum winds (R max ), and maximum wind speed (v max ) those that rotated counterclockwise from the axis of the Delaware Bay resulted in the highest WSEL (see Figure 12a). Tracks passing from the axis of the Delaware Bay (i.e., landfall equal to zero) were the only exceptions. Such orientation pushes the water in a more straightforward path toward the site. 2. For tracks with the same angle or orientation, c p , R max , and v max ; those with a landfall location spaced 1 R max west of the study area resulted in highest WSEL (see Figure 12b). Increasing the distance beyond that reduces the water levels. This could be attributed to locating the maximum winds closest to the AOI and is similar to findings from others, such as Irish et al. [58]. 3. Keeping all the parameters constant except R max , if the hurricane makes landfall on the west side of the site at a distance 1 R max away, the WSEL increases with increasing the R max (see Figure 12c).
However, this would not be the case for tracks located farther from the AOI. The same explanation as given in the above article again applies. 4. Results show a linear relation between the central pressure and storm surge (see Figure 12d).
Decreasing the central pressure increases the pressure drop, which eventually enhances maximum wind speeds and results in higher impacts during its passage. Each 1 millibar drop in the central pressure results in approximately 0.6 m increase in storm surge. Figure 12. Relationship between major hurricane parameters including central pressure, radius of maximum wind, landfall location, and orientation with water surface elevation (WSEL) at the area of interest (AOI). The legend of (a,c,d) shows the closest distance between the AOI and synthetic hurricane track normalized with Rmax. Same applies to horizontal axis of (b). Assigning a specific probability to the modeled water surface elevation was out of the scope of this study. However, the North Atlantic Coastal Comprehensive Study (NACCS, [87]) is used to approximately find the probability of the storm surge values presented here. NACCS is comprised of several models, including ADCIRC, WAve Model (WAM), and STeady state spectral WAVE (STWAVE). A total of 1050 synthetic tropical storms and 100 extratropical historical storms are simulated as part of the NACCS. The resolution of the unstructured grid varies from 30 to 50 m near the coast and 100 km offshore. Statistical analyses of past storms are employed in developing the synthetic storms. NACCS reports storm responses at 19,000 locations along the Atlantic Coast from frequent storms (e.g., 1-year AEP) to extremely rare (10 −4 AEP). NACCS reports a mean surge of 5.62 m NAVD88 at AOI with a 10 −4 AEP, which is approximately 35% lower than the surge estimated in this study. Therefore, using the same trend as those reported in NACCS, a surge of 7.52 m NAVD88 has, with high confidence, an AEP between 10 −6 and 10 −5 .

Conclusions
In this study, high-resolution numerical simulations of coupled storm surge, tides, and waves, forced by historical and synthetic tropical cyclones (TC) were used to examine the overall impact of low-probability hurricanes on a specific location within the lower Delaware Bay and estimate the probable maximum storm surge (PMSS). First, the combination of Overall, Nested1, and Nested2 Delft3D models were calibrated and validated with data from Hurricane Isabel (2003). WES, FLOW, and WAVE modules of the Delft3D software were calibrated by comparing the simulated values against observed ones. Then, they were used to simulate synthetic TCs that led to estimating a PMSS that has a very low probability of between 10 −6 and 10 −5 .
The greatest surge levels were experienced along the western boundary of the lower Delaware Bay during the passage of a synthetic hurricane that makes landfall west of that area. Model results showed that the combination of placing maximum winds over the area and pushing the water towards the area results in highest impact. It should be clarified that the percentage impact of waves on storm surge are considered out of the scope of this study.
For storms making landfall west of the AOI at a distance equal to R max , a positive logarithmic relation was found between increasing the distance and the maximum surge at the site. Li et al. and Irish et al. [56,88] also reported that the maximum storm surge increases with increasing hurricane R max .