Parameter sensitivity and uncertainty analysis for a storm surge and wave model

Development and simulation of synthetic hurricane tracks is a common methodology used to estimate hurricane hazards in the absence of empirical coastal surge and wave observations. Such methods typically rely on numerical models to translate stochastically generated hurricane wind and pressure forcing into coastal surge and wave estimates. The model output uncertainty associated with selection of appropriate model parameters must therefore be addressed. The computational overburden of probabilistic surge hazard estimates is exacerbated by the high dimensionality of numerical surge and wave models. We present a model parameter sensitivity analysis of the Delft3D model for the simulation of hazards posed by Hurricane Bob (1991) utilizing three theoretical wind distributions (NWS23, modified Rankine, and Holland). The sensitive model parameters (of 11 total considered) include wind drag, the depth-induced breaking γB, and the bottom roughness. Several parameters show no sensitivity (threshold depth, eddy viscosity, wave triad parameters, and depth-induced breaking αB) and can therefore be excluded to reduce the computational overburden of probabilistic surge hazard estimates. The sensitive model parameters also demonstrate a large number of interactions between parameters and a nonlinear model response. While model outputs showed sensitivity to several parameters, the ability of these parameters to act as tuning parameters for calibration is somewhat limited as proper model calibration is strongly reliant on accurate wind and pressure forcing data. A comparison of the model performance with forcings from the different wind models is also presented.


Introduction
We present a parameter sensitivity analysis of the Delft3D computer model under extreme storm conditions using Hurricane Bob (1993) as the underlying event.The analysis allows for an evaluation of the model's ability to reproduce observed values of water surface elevation and wave height which are relevant for storm surge hazard predictions.In addition, because publicly available wind observations do not provide sufficient information to drive the model, we evaluate the influence on the model performance of three widely used formulations to derive wind fields: NWS23, Holland, and Rankine.Finally, an assessment of the impact of the model grid resolution is also presented.
Specifically, this paper aims to (1) demonstrate the importance of model parameter selection in storm surge and wave modeling and (2) reduce the computational demand for producing surge and wave model parameter-related uncertainty estimates.

Uncertainty in storm surge hazard predictions
Hurricane hazards are commonly estimated from historical catalogs of coastal surge and wave characteristics.Walton (2000) provides a thorough review and discussion of these methods.The accuracy of this hazard analysis approach relies heavily on having an extensive continuous record of storm surges and waves.In many locations records of coastal surges exist only for durations much shorter than the return periods of interest.This shortfall of necessary data makes the development of hazard estimates of infrequent surges through this methodology unreliable.L. A. Bastidas et al.: Parameter sensitivity and uncertainty analysis for a storm surge and wave model In cases where empirical surge evidence is limited we may utilize alternate methods of estimating hurricane hazard.Irish et al. (2009Irish et al. ( , 2011a, b) , b) and Resio et al. (2009) demonstrate an approach which incorporates historical hurricane tracks and parameters to provide additional insight into hurricane surge hazard.They demonstrate that the joint probability method (JPM) of hurricane surge estimation may produce more reliable return period surge estimates than those based solely on empirical observations.Joint probability distributions of hurricane parameters are used to formulate synthetic probable hurricanes which are used to force numerical surge models.Research by Resio et al. (2013) and Tonkin et al. (2000) demonstrates that uncertainty with respect to hurricane parameterizations has a significant effect on hazard estimates as these methods still rely on empirical observations.Emanuel (2006), Emanuel et al. (2006a), and Vickery et al. (2000), among others, present methodologies which utilize physically based deterministic atmospheric models to simulate stochastically generated hurricane tracks.In this way deterministic models may estimate the feasible hurricane strength based on sea surface temperatures (SSTs) and atmospheric forcing at the storm boundaries.These synthetic tracks are generated in such a way that the statistical properties of the historical hurricanes can be confirmed to follow the patterns of historically observed hurricanes for a region.This general methodology has the advantage of not having to rely directly on empirical hurricane observations to produce estimates of potential future hurricanes.
Approaches based on stochastically generated synthetic hurricanes, e.g., Emanuel (2006), Emanuel et al. (2006a), Resio et al. (2009), and Vickery et al. (2000), are a promising path towards estimating hurricane storm surge risk where hurricane landfalls are infrequent and historical records are incomplete.We may use numerical surge and wave models to translate stochastically generated hurricane tracks into coastal hazard estimates.These methods can be modified further to assess nonstationary risk by incorporating the impacts of changing climate forcing (Emanuel et al., 2008;Grinsted et al., 2013;Lin et al., 2010Lin et al., , 2012) ) which the JPM intrinsically lacks.
Estimation methods based on stochastic hurricane tracks have the advantage of being able to calculate the hazard posed by infrequent hurricanes which exceed our relatively brief historical records (Emanuel et al., 2006a;Resio et al., 2009;Vickery et al., 2000).As these methods require deterministic surge and wave modeling they have the distinct disadvantage of having to consider uncertainties associated with numerical modeling of surge and waves.Numerical surge and wave models inherently introduce some additional uncertainty as they are imperfect recreations of true storm surge and wave physics.In utilizing these models we must face the problem of model formulation and parameter value selection.We then must translate this uncertainty in the numerical sim-ulation of hurricane storm surges into additional hurricane hazard uncertainty.

Surge and wave model parameter sensitivity
Previous studies have demonstrated that hydrodynamic model parameter uncertainty has a significant effect on coastal simulations, e.g., on sediment transport (Briere et al., 2010), on water quality (Li et al., 2013), on nearshore currents and wave growth (Adrani and Kaihatu, 2012), on tidal propagation (Mayo et al., 2014), on tsunami generation and propagation (Knighton and Bastidas, 2015;Sraj et al., 2014), and on storm surge (Ferreira et al., 2014;Holt et al., 2015).Despite these findings, several recent studies on the validation of the Delft3D model have not considered potential effects of uncertainty in model parameter values, e.g., Elias et al. (2000) and Golshani (2011).
Model parameter sensitivity and related uncertainty analysis methodologies typically rely on Monte Carlo simulations which follow these generalized steps: (1) a number of samples are drawn from the feasible parameter space to produce unique parameter sets; (2) these parameter sets are then evaluated with the numerical model to produce a model output; and (3) some form of the variance of the model output is evaluated and potentially related back to the parameter value variations.
The higher the dimensionality (number of parameters) of a model, the greater the number of simulations which are required to determine the effects each parameter has on a particular model output.Coastal surge models are typically highly parameterized formulations for wind-wave and surge modeling.They require extensive determination of appropriate numerical and physical settings and parameter values.Further, these models typically include many model elements (cells or nodes) at which the flow equations must be solved at each time step.These two considerations result in a large computational overburden when employing probabilistic sensitivity and model uncertainty estimates which can make the effort somewhat infeasible or impractical in practice.
Delft3D (Deltares, 2014a, b) is a model commonly used for simulating meteorologically induced coastal surges and wave growth.Delft3D combines a hydrodynamic model for large-scale simulation of water surface elevations and currents, Delft3D-FLOW (Deltares, 2014b), with a spectral wave model for the simulation of surface waves, Delft3D-WAVE (Deltares, 2014a).The large number of parameters and high computational demand of Delft3D makes formal storm surge sensitivity and uncertainty analysis often difficult to undertake.Deterministic models such as Delft3D require a number of inputs; ideally all inputs have some significant effect on the model output.Delft3D is a detailed model which has been designed to simulate a wide range of physical phenomena, capable of simulating tidal propagation over hundreds of kilometers as well as water quality fate and transport over several meters (Deltares, 2008).
We evaluate the possibility that simulations of storm surge and hurricane-induced waves may not depend equally on all Delft3D input parameters.To achieve this, we perform a model parameter sensitivity analysis of Delft3D storm surge and wave computations and identify the primary parameters of importance through the simulation of Hurricane Bob over a North Atlantic domain (Fig. 1).In this way, parameterrelated uncertainty estimates can be developed from a restricted parameter set, thereby reducing the overall computational demand for developing model uncertainty estimates.
In order to reduce the computational demand of this sensitivity analysis, we apply the Morris method (Campolongo et al., 2007;Morris, 1991).The Morris method is an efficient algorithm for computing model parameter elementary effects or changes in an output as a result of a change to a single parameter.In addition to estimating the elementary effects of each model parameter, the Morris method can produce an estimate of the parameter interaction with other parameters.In applying the Morris method, we can identify which model parameters have a significant effect on simulated storm surge and wave characteristics and which parameters have dependencies with other parameters or demonstrate significant nonlinearity.

North Atlantic storm surge
We select the US North Atlantic coast to evaluate Delft3D model parameter sensitivity because this region is somewhat reliant on numerical simulations for accurate hurricane hazard estimates.Historical hurricane tracks show few landfalling hurricanes of significant strength within the region (Dailey et al., 2009;AOML, 2015) and few coastal surge observations at tidal stations (NOAA, 2015b).A qualitative review of this empirical evidence may imply that hurricane storm surge is not a concern; however, recent research suggests that analysis using only empirical surge and hurricane parameter records is, at best, inconclusive.Dailey et al. (2009) evaluate the record of historical hurricane tracks against historical SSTs and show that a purely statistical approach based on hurricane observations results in a wide uncertainty for hurricane hazard forecasts for the US North Atlantic coast.Donnelly et al. (2001Donnelly et al. ( , 2004) ) estimate that five category 3 or greater hurricanes have made landfall along the US North Atlantic coast within the last 700 years based on coastal sedimentary records of Rhode Island and New Jersey.These estimates suggest that the past 60 years of coastal surge observations likely do not contain an observed storm surge resulting from a hurricane near the physical upper threshold of hurricane intensity (i.e., the probable maximum intensity).A similar finding is presented in Tonkin et al. (2000), in which L. A. Bastidas et al.: Parameter sensitivity and uncertainty analysis for a storm surge and wave model hurricane minimum central pressure measurements for the North Atlantic are shown to correlate poorly with SST measurements.Tonkin et al. (2000) suggest that this finding is most likely due to undersampling the joint distribution of hurricane central pressures and SSTs for the North Atlantic region within the past 60 years.Lin et al. (2010) estimate hurricane risk to New York City, USA, through a statistical/deterministic hurricane risk assessment methodology described by Emanuel (2006) and Emanuel et al. (2006b).Their research shows significant storm surges for New York with return periods of less than 500 years, which further demonstrates the potential shortcomings of relying on empirical surge and hurricane records for hazard estimation.Similarly, Lin et al. (2014) propose that the Atlantic Ocean may presently be experiencing a period of reduced hurricane activity.They propose that highenergy hurricane landfalls may be more common than that estimated from the extant historical hurricane track and surge records.
The effects of climate change and sea level rise add additional uncertainty to North Atlantic storm surge estimates.Villarini et al. (2011) evaluate whether anthropogenic forcing could increase the frequency of land-falling hurricanes within the region.They conclude that projected increases in hurricane frequency are not necessarily supported by statistical projections and note that significant uncertainty between analyses methods exist.Alternately, Lin et al. (2012) utilize stochastic deterministic hurricane surge modeling and estimate an increase in hurricane hazard estimates due to future climate forcing.
As there is great uncertainty surrounding the hurricane hazard estimates for the North Atlantic region, stochastic deterministic hurricane simulations are a promising path towards developing reliable hazard estimates (Lin et al., 2012(Lin et al., , 2014)).As such, we must acknowledge that numerical surge model parameter uncertainties will affect these estimates.To facilitate model parameter uncertainty estimates, we present a storm surge model parameter sensitivity analysis for Delft3D.As stated before, this paper aims to (1) demonstrate the importance of model parameter selection in storm surge and wave modeling and (2) reduce the computational demand for producing surge and wave model parameter-related uncertainty estimates.

Delft3D model description
We simulate two-dimensional, depth-averaged, unsteady flow characterizing hurricane wind and pressure setup with Delft3D-FLOW (Deltares, 2014b).The Navier-Stokes equations for incompressible flow are solved under the shallow water and Boussinesq assumptions.These equations are reduced to an implicit finite difference approximation through the Crank-Nicolson numerical scheme (Deltares, 2014b).The Delft3D-FLOW model was developed on a spherical grid at a 5 km spatial resolution and simulated at a time step of 60 s to satisfy the Courant-Freidrichs-Lewy condition of the Delft3D-FLOW solution technique.Though Delft3D-FLOW gives the users control over the solution technique, all simulations were performed with the Cyclic (Stelling and Leenderste, 1992) solution for the momentum equation.We perform all simulations with depth forced boundary conditions for open boundaries to reproduce tidal propagation, with 12 tidal forces components.
We simulate surface wind waves with Delft3D-WAVE, a derivative of the Simulating WAves Nearshore (SWAN) model (Deltares, 2014a).SWAN is a spectral wave model that evaluates the refracted wave height and wave angle based on linear wave theory (Booij et al., 1999;Deltares, 2014a).The SWAN model accounts for (refractive) propagation due to current and depth and represents the processes of wave generation by wind, dissipation due to white-capping, bottom friction, depth-induced wave breaking, and nonlinear wavewave interactions (both quadruplets and triads) (Booij et al., 1999;Deltares, 2014a).The SWAN model is based on the discrete spectral action balance equation and is fully spectral (across all directions and frequencies) (Dietrich et al., 2012).We use the same spatial grid for Delft3D-WAVE computations as was applied to the Delft3D-FLOW model.The spectral wave energy is computed at a 15 min time step using the nonstationary computational model.
We couple the Delft3D-FLOW and Delft3D-WAVE models for hurricane surge simulation at a 30 min time step.The wave forces computed in Delft3D-WAVE enhance the energy dissipation at the bed boundary layer in the storm surge model and generate a net mass flux affecting the current.These effects are accounted for by passing the radiation stress gradient determined from the computed wave parameters from Delft3D-WAVE to the Delft3D-FLOW model.The water levels and currents computed by the Delft3D-FLOW model are then passed back to the Delft3D-WAVE model for more accurate wave estimates (Deltares, 2014a).
Delft3D-FLOW and Delft3D-WAVE allow for considerable control of the hydrodynamic processes.Each model is highly parameterized.This allows the user to vary physical settings (e.g., wind drag coefficients, water density, gravitational constant, horizontal eddy viscosity, bottom roughness) as well as numerical settings (e.g., numerical solution technique, numerical convergence criteria, wetting drying thresholds).We evaluate the sensitivity of hurricane surge simulations to model parameters which have been considered to be classic calibration parameters as well as parameters which previous studies have demonstrated exert a significant effect on model uncertainty (Table 1).Each parameter is described in detail in the following sections.

Wind drag
The wind drag relationship defines the air water boundary condition for surge modeling.Surface winds exert a shear stress on the water surface which accelerates the water col-Nat.Hazards Earth Syst.Sci., 16,[2195][2196][2197][2198][2199][2200][2201][2202][2203][2204][2205][2206][2207][2208][2209][2210]2016 www.nat-hazards-earth-syst-sci.net/16/2195/2016/2009) present wind drag formulations as a function of surface wind speed.These studies suggest that wind drag increases linearly up to some wind speed termed the breakpoint velocity.Beyond this breakpoint wind speed, the drag coefficient reaches some limiting value or decreases slightly.Further research has demonstrated additional complexity suggesting the wind drag coefficient is also a function of the sea state (Andreas et al., 2012;Reichl et al., 2014) global location and temperature (Kara et al., 2007) and has some dependence on the bottom friction formulation (Johnson and Kofoed-Hansen, 2000;Zijlema et al., 2012).The considerable research that has been applied to estimating the proper wind drag coefficients to reproduce historical hurricanes demonstrates that there is some general agreement on the significance of this model input for accurate surge simulations (Bacopoulos et al., 2012;Cheung et al., 2007;Huang et al., 2013;Vatvani et al., 2012;Zachry et al., 2013).
Hereinafter, we consider the wind drag formulation to be a three-point function of the wind velocity, as described in Deltares (2014b).This results in a three-parameter model where we must determine the breakpoint wind speed (U B ), the breakpoint wind drag coefficient (C B ), and the terminal wind drag coefficient (C C ).The wind speed for the terminal wind drag (C C ) is fixed at 100 m s −1 .

Horizontal eddy viscosity
The horizontal eddy viscosity is a concept that primarily attempts to reproduce small-scale horizontal turbulent eddies and shear losses which cannot be simulated with a hydrodynamic model by utilizing a large computational grid size (Deltares, 2014b).These additional hydraulic losses are accounted for within Delft3D simulations through modification of a horizontal eddy viscosity term (ν H ). The larger the model grid, the more the smaller losses are neglected.The horizontal eddy viscosity term is considered a calibration parameter for Delft3D-FLOW, which is commonly a function of the model grid size (Deltares, 2014b).As we have selected a model grid resolution of 5 km, the horizontal eddy viscosity should be a significant consideration and is included in the sensitivity analysis.

FLOW bottom friction
The bottom friction formulation determines the frictional energy loss at the ocean bed boundary condition.Delft3D-FLOW and Delft3D-WAVE each require a separate selection of bottom friction formulation and parameter values.The formulation chosen for this research within Delft3D-FLOW is the spatially homogenous Manning's roughness.Delft3D-FLOW internally converts Manning's roughness values to a depth-dependent Chézy roughness for all computations (Deltares, 2014b).Previous research has demonstrated the Manning's roughness formulation is appropriate for simulation of the ocean bed boundary condition hydraulic losses and that this parameter has some effect on simulation results of long-wavelength wave propagation (Dao and Tkalich, 2007;Knighton and Bastidas, 2015;Sraj et al., 2014).

Threshold depth
The threshold depth parameter (D T ) is a numerical setting for Delft3D-FLOW which controls the wetting and drying of model cells.The threshold depth term specifically describes the depth below which a model cell will be considered dry and therefore excluded from the simulation.Medeiros and Hagen (2013) review different wetting and drying algorithms for hydrodynamic simulations including Delft3D.The threshold depth approach to cell wetting can result in artificial resistance to water propagation across cells and therefore may affect the model results in coastal areas.Selection of a threshold depth which is too small may result in numerical issues within a simulation.Horstman et al. (2013) demonstrate that the threshold depth within Delft3D-FLOW is a consideration for the simulation of tidal propagation.

WAVE bottom friction
We simulate wave energy dissipation by the ocean bed with the JoNSWAP (Hasselmann et al., 1973;Siadatmousavi et al., 2010) bottom friction formulation with a spatially homogenous friction coefficient (C JON ).
Several studies have identified the JoNSWAP parameter value as a significant consideration for the simulation of wave propagation within shallow water (Cialone and Smith, 2007;Johnson and Kofoed-Hansen, 2000;Mortlock et al., 2014;Padilla-Hernández and Monbaliu, 2001;Zijlema et al., 2012).The JoNSWAP bottom friction formulation has been historically considered to vary between two values representing swell conditions (0.038 m 2 s −3 ) and local wind-driven wave growth (0.067 m 2 s −3 ) (Hasselmann et al., 1973).Recently, van Vledder et al. (2010) suggested that the potential range of feasible bottom friction values may be more constrained than previously assumed.They demonstrate that the coefficient previously used to represent swell conditions may also more accurately reproduce bed dissipation for locally generated wind waves.

Depth-induced breaking
The depth-induced breaking model of Battjes and Janssen (1978) is used within Delft3D-WAVE spectral model to simulate the dissipation of waves within shallow water due to wave breaking (Booij et al., 1999).The depth-induced breaking along with the wave bed friction model determines the point of wave breaking and wave energy dissipation.The parameterization of this model requires estimates of the depthinduced breaking alpha (α B ) and gamma (γ B ) parameters.The α B parameter controls the rate of dissipation, whereas the γ B parameter controls the ratio of wave height to water depth at which wave breaking occurs.
It is acknowledged that more detailed depth-induced breaking models have been proposed which may represent an improvement over the current Delft3D-WAVE formulation.Filipot and Cheung (2012), Smit et al. (2013), and van der Westhuysen (2010) demonstrate potential limitations of the application of the SWAN model to coral reefs related to the reproduction of energy dissipation.

Nonlinear triad interactions
Wave triads simulate nonlinear wave-wave interactions.Wave-wave interactions occur when resonant wave frequencies exchange energy.This exchange transfers energy across the wave spectrum.The proportionality coefficient, α T , is a tunable parameter to modify the wave-wave interactions.The maximum frequency considered for wave-wave interactions is controlled by the β T parameter.
Nonlinear triad (three-wave) interactions have been shown to have a significant effect within shallow water (Beji and Battjes, 1993).Delft3D-WAVE incorporates nonlinear triad interactions through the lumped triad approximation (Eldeberky and Battjes, 1996).Akpınar et al. (2012) demonstrate that the parameterization of the triad model as an important consideration for simulation of waves over the Black Sea.

Hurricane Bob simulation
In this paper, Hurricane Bob (1991) is used as the primary model forcing data to estimate model parameter sensitivity.We chose Hurricane Bob for the following reasons: 1.The use of a historical hurricane allows us to compare model results with observed surges and waves.In this way we can determine not only the sensitivity of model outputs to parameter values but also which parameters enable Delft3D to accurately reproduce observations (i.e., serve as useful calibration parameters).
2. Hurricane Bob was a recent hurricane.The best-track data for this storm system are available at a higher temporal resolution than other historical category 2 landfalling hurricanes for the region (NOAA, 2015a).
3. Hurricane Bob is one of six hurricanes since 1950 to maintain a category 2 strength within 400 km of Boston, MA, USA (NOAA, 2015a).Hurricane Bob then quickly lost strength dropping to a tropical storm near Portland, Maine, USA (Fig. 1).This range of wind speeds within the study area allows a better exploration of the wind drag model parameterization of Delft3D (see Sect. 2.2.1).
4. Hurricane Bob traveled in a northeasterly direction along the US Atlantic coast (NOAA, 2015a) (Fig. 1).The track of this hurricane resulted in data being recorded at many tidal water level stations (NOAA, 2015b) and wave buoys (NOAA, 2015c).Hurricane Gloria (1985) had a similar strength and direction; however, Gloria made landfall in Connecticut, USA (NOAA, 2015a).
5. Cheung et al. (2007) show that an idealized Rankine wind field model of Hurricane Bob provides a reasonable representation of the storm.Their idealized wind model accurately reproduces observed wind velocity and pressure at land stations as well as coastal surge and wave characteristics based on the HURDAT (NOAA, 2015a) best-track observations.
Best-track data for Hurricane Bob were obtained from HUR-DAT (NOAA, 2015a).As noted in previous studies (Ling and Chavas, 2012;Resio et al., 2013;Taflanidis et al., 2011;Zhong et al., 2010), the recreation of hurricane wind and pressure fields from hurricane parameterizations can have a significant impact on model simulation results.Wind and pressure fields were developed for hurricane Bob using the NWS 23 (NOAA, 1979), the modified Rankine wind field as described in Cheung et al. (2007), and the Holland wind field (Holland, 1980).The NWS23 vortex model (NOAA, 1979) is an analytical formulation for reproduction of spatially distributed hurricane wind and pressure fields.The Holland vortex model is a modification to the analytical vortex model (Holland, 1980).The Holland B parameter, determined by the maximum wind speed and Coriolis forces, is used to modify the shape of the wind and pressure profiles of the hurricane (Holland, 1980(Holland, , 2008)).The modified Rankine model (Depperman, 1947) is based on a Rankine vortex which assumes solid body rotation within the radius of maximum winds (RMW) and a decaying wind speed inversely proportional to distance beyond the RMW.The modified Rankine model contains a tuning parameter, X, which we choose as 0.5 based on recommendations in Cheung et al. (2007) for Hurricane Bob.An adjustment for asymmetry of the wind field is applied to each model based on methods described by Jelesnianski (1966).

Observed surge and wave height
Hourly storm surge records from Atlantic City, Bar Harbor, Point Judith, Sandy Hook, and Woods Hole tidal stations (NOAA, 2015b) were used to evaluate the Delft3D-FLOW ability to reproduce coastal water surface elevations by varying model parameter values.
Hourly measurements from buoys 44007, 44008, 44013, and 44025 (NOAA, 2015c) were used to evaluate Delft3D-WAVE reproduction of significant wave heights.As noted in Table 2, the buoys available contain no measurements in true shallow water.In order to explore the depth dependence of wave parameter sensitives, we evaluate the model parameter sensitivity at the tidal gage stations.Though no measure can be given for reproduction of observed wave characteristics at these locations, we evaluate the effect of model parameter values on peak significant wave heights.Table 2  stations selected for model parameter sensitivity evaluation within this study.

Parameter sensitivity analysis
The Morris method (Campolongo et al., 2007;Morris, 1991) is a sensitivity analysis method that is particularly well suited to a model with significant computational overburden, as is the case here.The method does not need simplifying assumptions about the input/output behavior (Campolongo et al., 2000).The design is an efficient algorithm composed of individual randomized one-at-a-time designs, in which the impact of changing the value of each of the chosen parameters is evaluated in turn.A number of trajectories is initialized at a random position within the parameter space hypercube.Each move along the trajectory represents a change to one randomly selected parameter value.An estimate of the elementary effect of each model parameter is computed for each trajectory.Although different sampling schemes can be used, we follow the original Morris design (Morris, 1991).Overall, we used 50 trajectories, each one comprising 12 parameter sets, as we analyze 11 parameters, for a total of 600 simulations for each of the three wind models considered.Morris (1991) proposes two metrics that may be computed from the results of all trajectories.The mean of the elementary effects (µ) and the standard deviation of the elementary effects (σ ).Campolongo et al. (2007) suggest the use of the mean of the absolute elementary effects instead (µ * ).The µ and µ * parameters give an indication of the analyzed output sensitivity to a specified parameter.The σ parameter indicates nonlinearity in the model output response to changes in the model parameter or interdependencies between parameters.Hereinafter, we will refer to them as the Campolongo indices.For details of the method the reader is referred to Morris (1991) and Campolongo et al. (2007).
The output functions evaluated for the sensitivity analysis are chosen to allow for an evaluation of the hurricane hazard estimates which are commonly concerned with the peak flood elevation.For that reason we evaluate the sensitivity of peak wave height and peak surge elevation at each buoy and tidal station respectively.We also evaluate the parameter sensitivity for the entire simulation period by means of the root mean square error (RMSE) and the mean absolute error (MAE) with respect to the observed data.The RMSE represents an overall model error which emphasizes periods of large magnitude values (e.g., peak surge and wave heights).The MAE does not ascribe more weight to high values of model error as the RMSE.As shown in Fig. 2, the NWS23 (NOAA, 1979), Holland (Holland, 1980), and Rankine (Cheung et al., 2007) wind field models based on the hurricane best-track parameterization result in different wind forcing model inputs.The Rankine wind field model provides a more consistent match to wind speed observations as demonstrated by the RMSE at buoys 44013, 44008, 44007, and 44025.The Rankine model minimizes the error introduced by the forcing wind field at three of these specific locations.
The predicted wind directions are consistently similar for all three models.They are deemed an adequate representation of wind direction, which implies the best-track hurricane data are generally accurate.The peak winds at buoy 44025 arrive several hours earlier than the observed peak for all the models.It is assumed that this discrepancy may be related to an inaccurate position along the hurricane track from the best-track data.These incorrect forcing data impose some limitation on the model's ability to reproduce observed values at this location.

Delft3D-FLOW parameter sensitivity
Delft3D-predicted water surface elevations and significant wave heights show sensitivity with respect to the wind drag terms (U B , C B , C C ) and the bottom friction (n) (Figs. 3, 4).The bottom friction parameter has a significant influence only at the Bar Harbor station as this location is subject to large tidal oscillations.Stations with smaller tidal oscillations (< 1 m) show lesser effect of the bottom friction on peak surge elevation, RMSE, or MAE.Bottom friction formula- tion of Manning's n also had a significant effect on the wave height at Buoys 44007 and 44025.This effect is likely due to the wave buoys being located in shallower water than the other buoys and therefore more influenced by the bed friction (Table 2).
The wind drag parameters reveal significant sensitivity at Sandy Hook, Woods Hole, and Point Judith for peak surge elevation.The importance of the wind drag terms scales with proximity of the hurricane track to the tidal gage station and resulting surge elevation (Fig. 1).These same locations showed some sensitivity of the wind drag parameters to RMSE and MAE, but the sensitivity was somewhat reduced.These results suggest that the ability to properly calibrate these model parameters is more reliant on the quality of the wind forcing data applied as opposed to appropriate model parameter selection.The lack of sensitivity of the wind drag demonstrated at Atlantic City and Bar Harbor is ascribed to the Hurricane Bob causing only a minor surge at these locations.
The sensitive FLOW parameters all showed a significantly large value of σ (Figs. 3 and 4).Per the Morris method, this suggests a strong interaction among model parameters.This result is similar to the findings of Johnson and Kofoed-Hansen (2000) and Zijlema et al. (2012), who reported that the wind drag and bottom friction formulations have a shared dependency.Our results show that the dependency of these parameters must be considered when evaluating the effects of model parameter uncertainty and selecting appropriate parameter values for model calibration.The threshold depth parameter (D T ) and the horizontal eddy viscosity parameter (ν H ) have no discernable effect on the model output.We suggest that these parameters may be safely neglected in future hurricane hazard uncertainty studies, thereby reducing the computational demand.It should be noted that the D T parameter has numerical implications (Deltares, 2014b) and should still be carefully selected to avoid improper calculation of water surface elevations in areas with strong tidal oscillations.Within the present study any value within the numerically allowable range produced similar quality results.

Delft3D-WAVE parameter sensitivity results
The Delft3D-WAVE model parameterization is primarily related to shallow water processes where wave energy is dissipated due to wave-bed interactions.As such we see a spatially distributed set of model parameter sensitivities.At each NOAA wave buoy the simulated waves are primarily deep  water waves where the bed influence is minimal.At these locations the Delft3D-WAVE model predictions are insensitive to model parameter values.This finding implies that the existing NOAA buoys do not supply useful calibration infor-  Nat.Hazards Earth Syst.Sci., 16,[2195][2196][2197][2198][2199][2200][2201][2202][2203][2204][2205][2206][2207][2208][2209][2210]2016 www.nat-hazards-earth-syst-sci.net/16/2195/2016/mation for hurricane wave modeling.Along the coast at the tidal stations the predicted waves experience wave-bed interactions and therefore show greater sensitivity to the model parameters.
The wave parameters C JON and γ B had some minor effect on the peak surge elevation (Fig. 3).It has been previously shown that the wave setup can have some effect on storm surge predictions (Weaver and Slinn, 2004); however, these results demonstrate that the parameterization of the wave model does not play a significant role in predicting the peak surge elevation.The primary consideration here is that the wave model was coupled with the surge model to impart the appropriate wave stresses.
Delft3D-WAVE model predictions at NOAA buoys within deep water show significant sensitivity with respect to the C JON parameter, Manning's roughness, and the wind drag parameters (Fig. 4).Here we observe an almost identical parameter sensitivity with respect to wind-wave simulations.The depth-induced breaking γ B parameter showed some minor sensitivity.The WAVE model predicted that peak wave height is almost exclusively a function of the C JON parameter and the wind drag parameters.The additional parameters affecting model output only show up when evaluating the entire time series with RMSE and MAE.Within shallow water at the tidal stations, the predicted wave heights are primarily sensitive to the C JON and γ B parameters (Fig. 5).
Wave parameters showing sensitivity do not show an interaction among the wave parameters (Figs. 3,4 and 5).These feasible space of these parameters can be treated as marginal parameter spaces independent of other model parameters.
The β T , α T , and α B wave parameters had no significant effect on the simulated wave height.Selection of any parameter values within the allowable range for these parameters produced similar results.We therefore suggest that these parameters may be neglected for model calibration and uncertainty analysis.

Delft3D-FLOW simulation uncertainty for 5 km resolution model
As stated in Sect.2.5, in order to assess the model sensitivity, we ran Delft3D with 600 different parameter sets for each of the three wind models, i.e., a total of 1800 runs.The 600 samples provide a thorough coverage of the feasible parameter space, specified in Table 1, and can be used to assess the overall model performance and the associated parameterrelated uncertainty.
On Fig. 6 we present the entire set of 600 water surface elevations (ensemble) obtained from the simulations with 5 km resolution for each wind model at five tidal gage locations.The mean, the 50 and 95 % quantiles of the corresponding distribution are highlighted.They are picked to show the response from south to north over the domain following the track of Hurricane Bob.The ensemble results for wave height at the buoy locations are presented on Fig. 7.The error statis- tics for the mean at all the locations are also presented in Table 3 and on Fig. 8.
The results highlight that the model has a somewhat high level of precision, i.e., the bounds of the simulations are quite tight.The accuracy of the simulations, i.e., the bracketing of the observations, has some problems.For all the three wind models, at Bar Harbor the tidal amplitude during the simulation period is larger than the observed with an overestimation of the peak water surface elevation.There are also some timing errors on the peak value, particularly at the Point Judith location.Interestingly, the model shows some surge not observed in the data at the Sandy Hook location.It appears that the NWS23 yields a superior performance simulation at the Woods Hole location.The Holland model overestimates the peak value almost by a factor of 2.
Based on the error measures computed (Fig. 8, Table 3) the overall performance of the model with the NWS23 wind model seems to yield simulations that more closely resemble the observations at the Bar Harbor location by a significant margin.This is mostly related to the timings.The accuracy of the Rankine model outperforms the other two, except at Bar Harbor.This is most likely related to the best fitting of the wind fields using the Rankine model (Fig. 2).
The wave height simulations show a better performance for the Rankine model, with the Holland significantly overpredicting at buoys 44007 and 44013.Overall, it appears that for the chosen event and locations the Holland model shows the less accurate performance.

Delft3D-FLOW simulation uncertainty for multiple resolution model
A model with nested finer resolutions (∼ 500 m) around the location of the tidal gages was also setup to evaluate   polongo et al. (2007).Two additional locations were considered for the evaluation: Newport, Rhode Island, and Portland, Maine.The results of this ensemble of simulations are shown on Fig. 9 and the error statistics on Table 4 and Fig. 10.The only location with a significant improvement, over the coarse resolution, in model performance is Bar Harbor.The RMSE and MAE are reduced by almost a factor of 2. At this location, a significant increase in the precision of the simulations is also observed.This may be related to significant changes in the bathymetry.At the other locations, somewhat unexpectedly, there is actually a deterioration in the precision of the model.Improvements in the accuracy are also location dependent.For example, a deterioration in accuracy is observed at Sandy Hook.The improvement at the other locations, in terms of the errors, is marginal.
As for the coarse-resolution model, the Holland wind field shows the least accurate performance.It seems that the Holland model used here needs some tuning to improve the model responses.

Summary and Conclusions
In the present study we have used a sensitivity analysis methodology that is particularly suited for models with large computational overburden to determine the model parameter sensitivities for the case of hurricane-induced storm surges.An evaluation of the overall model performance, using a large ensemble, has been conducted which allowed for the determination of the overall model precision and accuracy.The results from the sensitivity analysis will allow for the reduction in the required number of simulations to calibrate the models.
Selection of the appropriate theoretical wind field model is a significant consideration for surge and wave modeling.The model parameters demonstrate similar sensitivity with different wind and pressure field forcing data; however, the ability of Delft3D parameters to function as calibration parameters for successful reproduction of storm surge and wave characteristics is largely dependent on proper wind and pressure forcing.
The Delft3D-FLOW model can be reformulated to a four parameter model for hurricane storm surge hazard simulation.The primary parameters of interest are U B , C B , n, and  The threshold depth parameter (D T ), horizontal eddy viscosity parameter (ν H ), nonlinear triad interaction parameters (α T , β T ), and depth-induced breaking alpha parameter (α B ) had no significant effect on the hurricane surge or wave hazard model output.The dimensionality, and therefore the computational overburden, of Delft3D storm surge and wave simulations can be reduced considerably.This is particularly important for probabilistic hazard estimates which require a significant number of simulations.
The sensitive model parameters showed significant nonlinearity in the model response and interactions among model parameters.Calibration of a Delft3D storm surge model should therefore consider the dependency of model parameters on each other.A traditional "one-at-a-time" calibration methodology may oversimplify the task of model calibration and could arrive at incorrect parameter value combinations.
Overall, Delft3D shows an ability to reproduce the water surface observations with reasonable precision and accuracy at most of the locations considered.However, the performance in terms of the wave height is of a lesser accuracy, with the precision significantly decreasing at the tail of the simulated event.As expected, the simulations are dependent on the wind fields driving the model.
For the specific locations used, the specific storm (Hurricane Bob), and with the pre-specified parameters for the L. A. Bastidas et al.: Parameter sensitivity and uncertainty analysis for a storm surge and wave model wind models, the Holland model produced an overall less accurate and less precise set of simulations.This suggests that some fine tuning of the wind field model parameters should be required in order to improve the quality of the simulations associated with a specific wind model.
We are currently working on the use of optimization algorithms for Delft3D calibration and identification of parameter value distributions, making use of the results presented here.

Data availability
All the data are available from the NOAA websites: NOAA Re-Analysis Project, Hurricane Research Division, NOAA Tide and Currents website, and National Buoy Data Center (NOAA, 2015a, b, c).The location identification of the data used is listed in the paper.See Table 4 for a listing of all the data locations used.

Figure 1 .
Figure 1.Delft3D model domain showing Hurricane Bob (1991) track, tidal stations, and wave buoys along the US North Atlantic Coast.

3
Results and discussion 3.1 Comparison of NWS23, Holland, and Rankine wind field forcing data We first present a comparison of the Delft3D storm surge and wave wind model forcing data.Though not a Delft3D model parameter, but rather an input forcing, selection of the wind field representation of a historical storm is a significant choice faced by modelers.Errors and uncertainty in the primary forcing data have a significant effect on model outputs.

Figure 3 .
Figure 3. Campolongo sensitivity indices for water surface elevations at tidal gage locations for different wind models and error function: blue upward pointed triangle is Max diff ; red left pointed triangle is RMSE; yellow right pointed triangle is MAE.Triangle size proportional to σ parameter.

Figure 4 .
Figure 4. Campolongo sensitivity indices for wave height at buoy locations for different wind models and error function: blue upward pointed triangle is Max diff ; red left pointed triangle is RMSE; yellow right pointed triangle is MAE.Triangle size proportional to σ index.

Figure 5 .
Figure 5. Campolongo sensitivity indices for wave height at buoy locations for different wind models: blue upward pointed triangle is NWS23; red right pointed triangle is Holland; yellow left pointed triangle is Rankine.Triangle size proportional to σ index.

Figure 6 .
Figure 6.Water surface elevation simulation results at tidal locations for different wind models for the 5 km resolution model.Error measures in meters.

Figure 7 .
Figure 7. Wave height simulation results at buoy locations for different wind models for 5 km resolution model.Error measures in meters.

Figure 8 .
Figure 8. Error performance measures, in meters, for mean of simulations with 5 km resolution.Water surface elevations at tidal gages (left panel) and wave heights at buoys (right panel).

Figure 9 .
Figure 9. Water surface elevation simulations results at tidal locations for different wind models with multiple resolutions.Error measures in meters.

Figure 10 .
Figure 10.Error performance measures for water surface elevation, in meters, for mean of multiple resolution simulations.

Table 1 .
Feasible parameter space for Delft3D model.
Powell et al. (2003), and Vickery et al. (n a wind set up (where wind setup is a component of the total surge) along coastal areas.Additionally, the wind stress applied over a fetch results in the growth of surface waves which are simulated through the spectral Delft-WAVE (SWAN) model.Surface waves shoal as they propagate into shallow coastal areas and can pose a flood hazard due to wave run-up and overtopping.Andreas et al. (2012),Donelan et al. (2004),Makin (2005),Powell et al. (2003), and Vickery et al. (

Table 2 .
Hurricane surge and wave observation locations.

Table 3 .
Error performance measures for mean of 5 km resolution simulations.

Table 4 .
Error performance measures for mean of multiple resolution simulations.