Effect of Bottom Friction , Wind Drag Coefficient , and Meteorological Forcing in Hindcast of Hurricane Rita Storm Surge Using SWAN + ADCIRC Model

An evaluation of the effect of bottom friction, wind drag coefficient, and meteorological forcing is conducted using a tightly coupled wave and circulation model, SWAN + ADCIRC (i.e., Simulating WAves Nearshore + ADvanced CIRCulation), to hindcast the storm surge of Hurricane Rita (2005). Wind drag coefficient formulations of Powell, Zijlema, and Peng & Li are used to calculate wind stresses. Bottom friction and wind drag coefficients are systematically increased and decreased to quantify their impacts on the hindcast. Different meteorological forcing options are applied to study the effect of wind fields on storm surge development and propagation. Simulated water levels are compared with observed data collected from about 150 locations. It is evident that a lower bottom friction causes higher and faster surge propagation, and earlier arrival of inundation peak at locations far from the land fall. Drag coefficients of Powell, with or without a cap of 0.002, and Zijlema produce similar results, while that of Peng & Li slightly overpredicted the surge. Wind fields may cause overprediction or underprediction of the surge, depending on the choice of the wind model. A good agreement is found between Zijlema’s findings and this study; that simultaneously decreasing or increasing both bottom friction and wind drag essentially provides the same hindcast results.


Introduction
Storm surge due to tropical cyclones is an above-normal coastal water level caused by strong winds, atmospheric pressure differences, and wind-driven wave setup [1].Impacts of storm surges are exacerbated by astronomical tides, precipitation runoff, and geostrophic effects.As tropical cyclones approach land, they encounter varied geographical features, such as rivers, barrier islands, levees, raised roads and railroads, and low-lying topography [1][2][3].Some of these features amplify the storm surge, and others attenuate it through bottom friction and surface roughness.Storm surge magnitude is also substantially affected by the shallowness and orientation of the water body relative to storm path, the storm's intensity and size, and the timing of tides.
Recent occurrences of major hurricanes, such as Charley, Frances, and Ivan (2004); Katrina, and Rita (2005); Gustav and Ike (2008); and Isaac and Sandy (2012) and their large impacts call for better understanding of the sensitivity of coastal storm surges to different forcings and parameters that are needed for accurate surge predictions [4][5][6][7].Despite recent advancements, there are areas of storm surge modeling that can be improved, such as bathymetry and topographic data, especially accounting for erosion and sedimentation due to storm surges, waves, and winds [3]; wind field generation, especially accounting for the asymmetric nature of the storm [8][9][10][11][12]; drag coefficients that account for sea state and surface roughness [2,13]; directional bottom friction formulation [14]; and wind-sea interaction formulation [10, [15][16][17].
The parameterization of bottom friction, surface roughness, and surface canopy is typically done as a function of the physical landscape data [18], which contain some level of measurement uncertainties [19].Furthermore, some coastal areas and wetlands are subject to continuous erosion and sedimentation due to frequent storm surges, waves, and winds.These activities are likely to alter the bathymetry and topography of coastal areas over time.Increased bottom friction may attenuate the storm surge propagation in wetlands [3], while the opposite is true for decreased bottom friction.A correct bottom friction based on actual sea floor material is needed to help modeling some unique physical aspects of storm surge, such as a forerunner, which is a gradual rise of the water level along the coast preceding the arrival of hurricane [20].While the importance of bottom friction is well known, a systematic study on the quantitative effect of low and high bottom friction on storm surges is not available in open literature.Therefore, the first objective of the present study is to systematically vary the bottom friction and quantify its effect on the hindcast of Hurricane Rita's surge.
Hurricane winds are the primary driver of storm surge.Hurricane track forecasting has improved significantly over the years, but the forecasting of hurricane intensity still contains significant uncertainties [21].A variety of numerical and statistical models have been developed for forecasting (e.g., [22][23][24]) and hindcasting hurricane wind fields (e.g., [25,26]).The most common approaches used in surface wind modeling for tropical cyclones may be categorized as: (1) Simple analytical parametric models, such as the original Holland [23]; (2) Dynamical models such as the Planetary Boundary Layer (PBL) model of Chow [27] as later implemented by Cardone et al. [28], Shapiro [29], Thompson and Cardone [30], and Vickery et al. [31]; the Dynamic Holland Model [32] and the Asymmetric Holland Model [22,33]; (3) Full scale physics-based dynamical models such as Fifth-Generation Penn State/NCAR Mesoscale Model (MM5) [34], Geophysical Fluid Dynamics Laboratory (GFDL) model [35], and Weather Research and Forecasting (WRF) model [36]; and (4) Kinematical methods, most notably Hurricane WIND (HWIND) [37] Interactive Objective Kinematic Analysis (IOKA) by Oceanweather Inc (OWI) [11,38].The second objective of the present study is to use some of the popular wind models in hindcast of Hurricane Rita's surge and assess their performance against the observed data.
The storm surge model needs wind stress, a measure of the air-sea momentum exchange, which is derived from the wind field generated through one of the above approaches.Unless an explicit atmospheric, wave, and ocean coupled model is used [10], calculation of wind stress requires specification of a wind drag coefficient (Cd), which is an active area of research [16].Most drag coefficients are calculated based on wind speed only, although extensive research suggests that sea conditions, fetch shoaling, wave breaking, and sea spray and foam may significantly influence the drag coefficient ( [16] presents a review of some of those studies).
A popular linearly varying drag coefficient is provided by Garratt [39], which has been the default option in ADCIRC [40][41][42] for a long time.Among many other formulations, the quadratic profiles of Zijlema et al. [15] and Peng & Li [43] are reported to perform well in hurricane wind speed ranges.Powell [44] analyzed 2664 Global Positioning System (GPS) dropsonde profiles collected between 1997 and 2005 that included some extreme hurricane datasets.He found that drag coefficient depends on the azimuthal angle of the hurricane translational motion, and proposed formulations that depend on sectors of hurricanes (e.g., right, left, rear, etc.).Powell's drag coefficient is arguably the most sophisticated one available in open literature.The third objective of the present study is to assess the performance of different wind drag options in the hindcast of Hurricane Rita's surge.In addition, Zijlema et al. [15] argued that both lower wind drag and bottom friction and higher wind drag and bottom friction essentially provide the same wave and storm surge hindcast results.The present study re-examines the validity of this argument.
The rest of the paper is organized as follows: Section 2 presents model details, Section 3 presents simulation cases and capability criteria, Section 4 presents results and discussion, and finally, Section 5 presents concluding remarks.

Model Details
The SWAN + ADCIRC model is used to hindcast the storm surge of Hurricane Rita (2005).Hurricane Rita formed on 18 September and hit the coast of southwestern Louisiana as a Category 3 at 0740 Coordinated Universal Time (UTC) 24 September 2005 [45].On arrival into the Gulf of Mexico, Rita rapidly strengthened to a Category 5. Gradually, it weakened to Category 4 and then maintained a Category 3 status up to the time of landfall between Sabine Pass, Texas and Holy Beach, Louisiana [45].The spread of the surge extended from the Mississippi coast in the east and Port O'Connor region (a small unincorporated village between Galveston and Corpus Christi in Texas) in the west.
Rita hindcasts done by Kerr et al. [1] using SWAN + ADCIRC are taken as the default case in the present study.However, the semi-implicit solver option of ADCIRC is adopted in the present study, instead of the lumped-explicit one used by Kerr et al.The summary of the default model settings is reviewed here (details can be found in Kerr et al. [1]).The ULtralite-Levee-Removed (ULLR) mesh (417,642 nodes, 826,866 linear triangular elements) is used to run SWAN + ADCIRC (see Figure 1a).The region of interest along with Hurricane Rita's track is shown in Figure 1b.A spatially varying distribution of Manning's roughness coefficient derived from land-use databases [1] is used for computing the bottom friction in ADCIRC.The same spatially varying Manning's n coefficients are exported to use in SWAN.The rest of the paper is organized as follows: Section 2 presents model details, Section 3 presents simulation cases and capability criteria, Section 4 presents results and discussion, and finally, Section 5 presents concluding remarks.

Model Details
The SWAN + ADCIRC model is used to hindcast the storm surge of Hurricane Rita (2005).Hurricane Rita formed on 18 September and hit the coast of southwestern Louisiana as a Category 3 at 0740 Coordinated Universal Time (UTC) 24 September 2005 [45].On arrival into the Gulf of Mexico, Rita rapidly strengthened to a Category 5. Gradually, it weakened to Category 4 and then maintained a Category 3 status up to the time of landfall between Sabine Pass, Texas and Holy Beach, Louisiana [45].The spread of the surge extended from the Mississippi coast in the east and Port O'Connor region (a small unincorporated village between Galveston and Corpus Christi in Texas) in the west.
Rita hindcasts done by Kerr et al. [1] using SWAN + ADCIRC are taken as the default case in the present study.However, the semi-implicit solver option of ADCIRC is adopted in the present study, instead of the lumped-explicit one used by Kerr et al.The summary of the default model settings is reviewed here (details can be found in Kerr et al. [1]).The ULtralite-Levee-Removed (ULLR) mesh (417,642 nodes, 826,866 linear triangular elements) is used to run SWAN + ADCIRC (see Figure 1a).The region of interest along with Hurricane Rita's track is shown in Figure 1b.A spatially varying distribution of Manning's roughness coefficient derived from land-use databases [1] is used for computing the bottom friction in ADCIRC.The same spatially varying Manning's n coefficients are exported to use in SWAN.In this study, a wind drag coefficient due to Powell [44] with a cap of Cd ≤ 0.002 [1] is used in ADCIRC, and this is the default setup.Typically, SWAN applies a wind drag coefficient due to Wu [46] with a cap of Cd ≤ 0.0035 [6]; but for this study, SWAN uses whatever wind drag coefficient used in ADCIRC [1].Drag coefficient formulations of Zijlema [15] and Peng & Li [43] are used to investigate how these affect hurricane storm surges.Figure 2a shows how Powell, Zijlema and Peng & Li drag coefficient profiles vary with respect to wind speed.The drag coefficient profile of Zijlema typically resides in the lower envelope of most profiles available in the literature [16], while that of Peng & Li resides at the upper envelope of the lot. Figure 2b displays the 75% and 125% of Powell's sector based drag coefficient profiles.The main tidal constituents, K1, O1, Q1, P1, M2, K2, N2, and S2 are generated from the EC2001 tidal database [47,48].In this study, a wind drag coefficient due to Powell [44] with a cap of Cd ≤ 0.002 [1] is used in ADCIRC, and this is the default setup.Typically, SWAN applies a wind drag coefficient due to Wu [46] with a cap of Cd ≤ 0.0035 [6]; but for this study, SWAN uses whatever wind drag coefficient used in ADCIRC [1].Drag coefficient formulations of Zijlema [15] and Peng & Li [43] are used to investigate how these affect hurricane storm surges.Figure 2a shows how Powell, Zijlema and Peng & Li drag coefficient profiles vary with respect to wind speed.The drag coefficient profile of Zijlema typically resides in the lower envelope of most profiles available in the literature [16], while that of Peng & Li resides at the upper envelope of the lot. Figure 2b displays the 75% and 125% of Powell's sector based drag coefficient profiles.The main tidal constituents, K1, O1, Q1, P1, M2, K2, N2, and S2 are generated from the EC2001 tidal database [47,48].The OceanWeather Inc. (OWI) structured and data-assimilated wind and pressure fields are used as the default wind model, similar to Kerr et al. [1].To study the effect of other meteorological forcing options, the following meteorological setups are used in instead of the OWI model: (1) NOAA National Hurricane Research Division (NHRD) HWIND Model; (2) Planetary Boundary Layer (PBL) Model; (3) Dynamic Holland Model (HM); and (4) Asymmetric Holland Model (AHM).Some technical details of these models are presented below.
The OWI wind field for Rita was generated by blending of TC96 (short for Thompson and Cardone, 1996) mesoscale model [30], an inner core wind field transformed to 30 min averaged sustained winds with gulf scale winds using the Interactive Objective Kinematic Analysis (IOKA) system [11,38].The hindcast winds use measurements from anemometers, airborne and land-based Doppler radar, microwave radiometers, buoys, ships, aircrafts, coastal stations, and satellite measurements.
HWIND wind field "snapshots" are generated at 6-h intervals after aircraft reconnaissance missions in a given system is done [49].HWIND incorporates a comprehensive set of observational data sources spanning more than thirty land, air, and sea-based platforms, such as aircraft reconnaissance, dropsondes, airborne Doppler radar wind measurements, buoys, satellites, and coastal observation networks.HWIND datasets include integrated kinetic energy, which is considered to be the index for measuring tropical cyclone size and intensity, and an alternative to the Saffir-Simpson Scale.One drawback of HWIND is the lack of wind field "snapshot" data soon after the hurricane makes landfall, which is needed for storm surge hindcasting.For the present study, the latest available HWIND snapshot is extrapolated in time with adjustments in storm position and a gradual reduction factor of the hurricane intensity.For Hurricane Rita, the last snapshot available is at 24 September, 10:30 AM.For the next 24 h, this snapshot is used by replacing its eye location with the current eye location and assigning an intensity reduction factor less than 1 (e.g., 0.9, 0.6, 0.5, and 0.4 for 24 September noon, 6 PM, 25 September midnight, and 6 AM, respectively, are used in the present study).These reduction factors are approximated from the ratios of maximum wind speeds at those times reported in best track data in comparison to the maximum wind speed on 24 September, 10:30 AM.
The PBL model is an application of a theoretical model of the horizontal airflow in the boundary layer of a moving vortex [27].The model solves the vertically averaged equations of motion that govern boundary layer subject to horizontal and vertical shear stresses.The center of the storm translates with a constant velocity.The pressure gradients and velocities of the storm in the radial direction are calculated from empirical equations that are devised based on numerous field data and analysis of wind profiles.
Hurricane best track or forecast data contains some basic information, including eye location and time, maximum wind speed and radius, and central pressure.The Dynamic Holland model [23] calculates some parameters from those data to apply in empirical equations to calculate the atmospheric pressure and gradient wind velocity, from which wind velocity at 10 m height is calculated [32].This model assumes symmetric distribution of the wind field.However, hurricanes often have asymmetric wind field distributions.The Asymmetric Holland Model (AHM) uses either R64, R50, or R34 distance to strongest wind isotach (i.e., 64 kt, 50 kt, or 34 kt away from the eye) to solve for a different maximum radius in each storm quadrant (NE, NW, SW, SE) [22,33].In other words, AHM uses the single strongest isotach in each quadrant.The AHM is a significant improvement over the Dynamic Holland Model.

Simulation Cases and Capability Criteria
A total of 12 simulation cases are performed (Table 1).Cases 1-5 are run with SWAN + ADCIRC, varying the bottom friction and wind drag formulations.As mentioned earlier, Case 1 is the default simulation which has the same parameter attributes as the Rita storm surge hindcast reported by Kerr et al. [1].Case 2 has similar parameters as Case 1 except that Manning's n is reduced to 60%, while for Case 3, Manning's n is increased to 150%.Note that, an attempt to reduce Case 1 to 50%, the simulation was unsuccessful due to unrealistic water elevation.Cases 4 and 5 use drag formulations by Zijlema, and Peng & Li, respectively.Cases 6-9 use the same setup as the Case 1, except that wind fields are generated by HWIND, PBL, AH, and AHM wind models, respectively.To examine combined effect of wind drag and bottom friction, the default case parameters are carefully reduced to 75% in Case 10 and increased to 125% in Case 11.Case 12 is set up similar to Case 1, but with a non-capped drag coefficient.Each simulation is cold started on 0000 UTC 13 August 2005, with a 36-day tides-only period that allows the tides to reach a dynamic equilibrium.This is followed by a 7.2-day Rita simulation from 0000 UTC 18 September 2005 to 0500 UTC 26 September 2005.Qualitative and quantitative comparisons using water level time series and high-water mark (HWM) of simulated results and observed data are presented for the station locations (see Table 2).These observation station locations with respect to Hurricane Rita track are shown in Figure 1b.Simulation capability is quantified using the capability criterion used in Kerr et al. [1].These are repeated here for clarity: Coefficient of Determination (R 2 is a number that indicates how well data fit a regression line, with an ideal value of one), Root Mean Square Error (E RMS is a measure of the magnitude of error, with an ideal value of zero), Mean Error (E), Slope of the best fit line (m, with an ideal value of one), Mean Normalized Bias (B mn is a measure of the model's magnitude of overprediction or underprediction normalized to the observed value, with an ideal value of zero): Scatter Index (SI, which is the standard deviation normalized by the mean observed value, with an ideal value of zero): Mean Absolute Error (MAE), and Mean Normalized Error (E NORM , which is the mean error normalized by the mean observed value, with an ideal value of zero): where O is the observed value, E is the error in terms of simulated minus observed and N is the number of data points.Throughout the simulations, some of the HWM's locations are predicted to be dry, which varies from case to case.In the statistical analyses of HWMs, only wet locations predicted by simulations are considered.

Results and Discussion
The hindcast simulations of Hurricane Rita storm surge using different model setups (see Table 1) are performed and compared against each other and with observed data.Since Case 1 is the default case in this study and is reported by Kerr et al. [1] as well, it is taken as a benchmark for model comparisons.Figure 3a shows maximum wind track created using the OWI model that is used in Case 1.The track line is also shown on the color plot.It is clear that the right sector of Hurricane Rita had higher wind velocity than the left one.Figure 3b shows the wind vector at 8:30 AM on 24 September 2005.The anticlockwise rotation of the hurricane is obvious.Wind at the right of the track moves northward, while wind at the left moves southward.Wind is the primary force here that creates the water motion and surge.Unique geographical features around the landfall of Hurricane Rita produce interesting water elevation contours and velocity vectors.Water is pushed against the shoreline near Pecan Island to cause about 5 m surge at the right of landfall, while water at the left side could slide past the coast, near the Galveston area, very easily, causing not more than 3 m surge there.Figure 3c shows water elevation and velocity vector at 8:30 AM on 24 September 2005.

Effect of Bottom Friction and Drag Coefficient Formulation
The effects of reducing and increasing bottom friction are shown in the maximum water elevation and velocity contour plots in Figure 4. To assess differences, both maximum water elevation and velocity were subtracted from those of the default case, Case 1. Figure 4a shows that Case 2 with 60% of default bottom friction caused a higher water elevation (i.e., negative differences), as expected.
On the eastside of Rita landfall near Pecan Island, water is pushed against the shoreline, and Case 2 produced 1.5 m higher water elevation than Case 1.On the other hand, since west of the landfall location water was able to flow southward uninterrupted, elevations did not reach more than 0.5 m higher than in Case 1.However, the water velocity of Case 2 reached 1.5 m/s higher than that of Case 1 near the Galveston area, as shown in Figure 4b.Increased bottom friction caused decreased water elevation due to increase in energy dissipation at the sea floor.As shown in Figure 4c,d), Case 3 with 150% bottom friction of Case 1, produced 1 m lower water elevation near Pecan Island, and 1 m/s slower water velocity near the Galveston area.Far from the shallow water regions, the effects of lower or higher bottom friction were distinguishable.

Effect of Bottom Friction and Drag Coefficient Formulation
The effects of reducing and increasing bottom friction are shown in the maximum water elevation and velocity contour plots in Figure 4. To assess differences, both maximum water elevation and velocity were subtracted from those of the default case, Case 1. Figure 4a shows that Case 2 with 60% of default bottom friction caused a higher water elevation (i.e., negative differences), as expected.On the eastside of Rita landfall near Pecan Island, water is pushed against the shoreline, and Case 2 produced 1.5 m higher water elevation than Case 1.On the other hand, since west of the landfall location water was able to flow southward uninterrupted, elevations did not reach more than 0.5 m higher than in Case 1.However, the water velocity of Case 2 reached 1.5 m/s higher than that of Case 1 near the Galveston area, as shown in Figure 4b.Increased bottom friction caused decreased water elevation due to increase in energy dissipation at the sea floor.As shown in Figure 4c,d), Case 3 with 150% bottom friction of Case 1, produced 1 m lower water elevation near Pecan Island, and 1 m/s slower water velocity near the Galveston area.Far from the shallow water regions, the effects of lower or higher bottom friction were distinguishable.Wind stress is calculated from the wind speeds and wind drag coefficients.The drag coefficient formulations are predominantly wind speed-dependent empirical equations.To assess differences produced by different drag coefficients, both maximum water elevation and velocity were subtracted from these of the default case, Case 1. Figure 5a,b shows the maximum water elevation and velocity differences between Cases 1 and 4, using Powell's drag coefficient with a cap of 0.002 and Zijlema's formulation, respectively.East of landfall, Powell predicts a 0.3 m higher water elevation than Case 1, due to its slightly higher drag coefficient values.Powell produces a 0.5 m/s higher water velocity on and at the left side of the track.However, the difference in the drag coefficients near the shallow water regions is not significant.Peng & Li has consistently higher drag coefficients than Powell's.Figure 5c shows that Peng & Li produces an approximately 0.5 m or higher inundation near Pecan Island areas.However, the differences in velocity field between Cases 1 and 5 are indistinguishable, as shown in Figure 5d.Since water can easily flow southward on the west side of the track, Peng & Li's slightly higher drag coefficients do not produce any significant differences in elevation and velocity in comparison to those of Case 1.
water regions is not significant.Peng & Li has consistently higher drag coefficients than Powell's.Figure 5c shows that Peng & Li produces an approximately 0.5 m or higher inundation near Pecan Island areas.However, the differences in velocity field between Cases 1 and 5 are indistinguishable, as shown in Figure 5d.Since water can easily flow southward on the west side of the track, Peng & Li's slightly higher drag coefficients do not produce any significant differences in elevation and velocity in comparison to those of Case 1.The SWAN + ADCIRC water level time series are compared to observed water levels at stations (see Figure 1b).Away from landfall at station (J), Figure 6a indicates that the Case 1 simulation accurately predicted the surge, but with an early peak arrival.Kerr et al. [1] indicated a possible error in the hydraulic connectivity in the marsh as the reason.The station seemed to be inundated hours before the arrival of the peak, whereas the simulation predicted zero water elevation at the beginning and a sudden rise of water at a later time.With the reduction of the bottom friction to 60%, the Case 2 simulation result shows the peak arriving an hour earlier than Case 1.This finding is in agreement with the fact that the simulation of forerunner surges requires accounting for low bottom friction due to sand/gravel and mud in the ocean floor [20].In addition to early peak arrival, the peak is overpredicted significantly for Case 2. The result of Case 3 with a 150% increase in bottom friction shows that there is a phase lag, i.e., the peak arrives slightly later than the observed peak.The peak The SWAN + ADCIRC water level time series are compared to observed water levels at stations (see Figure 1b).Away from landfall at station (J), Figure 6a indicates that the Case 1 simulation accurately predicted the surge, but with an early peak arrival.Kerr et al. [1] indicated a possible error in the hydraulic connectivity in the marsh as the reason.The station seemed to be inundated hours before the arrival of the peak, whereas the simulation predicted zero water elevation at the beginning and a sudden rise of water at a later time.With the reduction of the bottom friction to 60%, the Case 2 simulation result shows the peak arriving an hour earlier than Case 1.This finding is in agreement with the fact that the simulation of forerunner surges requires accounting for low bottom friction due to sand/gravel and mud in the ocean floor [20].In addition to early peak arrival, the peak is overpredicted significantly for Case 2. The result of Case 3 with a 150% increase in bottom friction shows that there is a phase lag, i.e., the peak arrives slightly later than the observed peak.The peak is reduced significantly as well.All indicate that bottom friction, mesh resolution, and other features may have contributed to the mismatch between simulated results and observed data.
In the region of landfall, Figure 6b-e at stations (H, D, I, and F), Case 1 overpredicts water the elevation peak by 0.5 m to 1 m at most stations, some of which Kerr et al. [1] attributed to improper road resolution.Lowering the bottom friction, as in Case 2, shows a 1 to 2 m higher water peak than the observed peak.Increasing bottom friction, as in Case 3, shows almost the same water peak level as the observed data, but with about one-hour phase lag at several stations.
At inland stations (E, C, and G) (Figure 6f-h), the default Case 1 overpredicted the peak water levels with a phase lead.With the bottom friction reduced to 60%, Case 2 overpredicts the maximum water peak even more, along with a larger phase lead, which is expected.When modeled with friction increased to 150%, as in Case 3, the peak and phase mostly match with the observed data.
The choice of wind drag formulation (i.e., Cases 4 and 5) in all stations appears to have a very negligible effect, and these results are similar to those of Case 1.However, station locations (look for black dots) superimposed on Figure 5a,c suggest that Case 4 (Zijlema) should produce about 0.3 m lower elevation, and Case 5 (Peng & Li) should produce about 0.5 m higher elevation, than Case 1.This apparent anomaly is due to the fact that output files used for these time series plots (i.e., Figure 6) are created intermittently every 10 min.The maximum elevation and velocity data used in Figure 4 are filtered from all time steps during the actual simulation.The peak elevation durations in Cases 4 and 5 probably are due to short lived waves, which are not captured in time series.
to sand/gravel and mud in the ocean floor [20].In addition to early peak arrival, the peak is overpredicted significantly for Case 2. The result of Case 3 with a 150% increase in bottom friction shows that there is a phase lag, i.e., the peak arrives slightly later than the observed peak.The peak is reduced significantly as well.All indicate that bottom friction, mesh resolution, and other features may have contributed to the mismatch between simulated results and observed data.
In the region of landfall, Figure 6b-e at stations (H, D, I, and F), Case 1 overpredicts water the elevation peak by 0.5 m to 1 m at most stations, some of which Kerr et al. [1] attributed to improper road resolution.Lowering the bottom friction, as in Case 2, shows a 1 to 2 m higher water peak than the observed peak.Increasing bottom friction, as in Case 3, shows almost the same water peak level as the observed data, but with about one-hour phase lag at several stations.
At inland stations (E, C, and G) (Figure 6f-h), the default Case 1 overpredicted the peak water levels with a phase lead.With the bottom friction reduced to 60%, Case 2 overpredicts the maximum water peak even more, along with a larger phase lead, which is expected.When modeled with friction increased to 150%, as in Case 3, the peak and phase mostly match with the observed data.
The choice of wind drag formulation (i.e., Cases 4 and 5) in all stations appears to have a very negligible effect, and these results are similar to those of Case 1.However, station locations (look for black dots) superimposed on Figure 5a,c suggest that Case 4 (Zijlema) should produce about 0.3 m lower elevation, and Case 5 (Peng & Li) should produce about 0.5 m higher elevation, than Case 1.This apparent anomaly is due to the fact that output files used for these time series plots (i.e., Figure 6) are created intermittently every 10 min.The maximum elevation and velocity data used in Figure 4 are filtered from all time steps during the actual simulation.The peak elevation durations in Cases 4 and 5 probably are due to short lived waves, which are not captured in time series.The peak simulated water levels are compared for each of the simulations to measured high water marks (HWMs) [50] in Figure 7. Just as in Kerr et al. [1], only wet stations are used for all simulated HWM analysis.Since all plots look similar, except for Case 2, only Cases 1 and 2 are presented in Figure 7.The summarized statistics are given in Table 3.It is worth noting that Kerr et al. [1] obtained an R 2 value of 0.663 using the lumped-explicit solver option in the SWAN + ADCIRC model to hindcast Rita's storm surge.The present study, which uses the semi-implicit solver option, obtained 0.704 for SWAN + ADCIRC (Case 1) models, as shown in Table 3.With the bottom friction decreased, Case 2 has the lowest correlation coefficient, highest mean error, Bias, SI and MAE, although most of the HWM spots are inundated due to low bottom frictions.With 20 dry locations, Case 3 has the lowest mean error, Bias, SI, MAE, ENORM and a high correlation coefficient.Drag coefficient of Zijlema (Case 4) shows slightly better statistics than the Powell's (Case 1).Peng & Li's drag formulation (Case 5) shows lower correlation coefficient; higher mean error, Bias, SI, MAE, and ENORM than those of Powell's (Case 1), but it has a higher number of wet stations.It is obvious that the correlation coefficient alone may be a misleading criterion to judge HWM comparison.A simulation may produce less inundation extent, yet score a high correlation coefficient since dry nodes are typically ignored from the statistics.The peak simulated water levels are compared for each of the simulations to measured high water marks (HWMs) [50] in Figure 7. Just as in Kerr et al. [1], only wet stations are used for all simulated HWM analysis.Since all plots look similar, except for Case 2, only Cases 1 and 2 are presented in Figure 7.The summarized statistics are given in Table 3.It is worth noting that Kerr et al. [1] obtained an R 2 value of 0.663 using the lumped-explicit solver option in the SWAN + ADCIRC model to hindcast Rita's storm surge.The present study, which uses the semi-implicit solver option, obtained 0.704 for SWAN + ADCIRC (Case 1) models, as shown in Table 3.With the bottom friction decreased, Case 2 has the lowest correlation coefficient, highest mean error, Bias, SI and MAE, although most of the HWM spots are inundated due to low bottom frictions.With 20 dry locations, Case 3 has the lowest mean error, Bias, SI, MAE, E NORM and a high correlation coefficient.Drag coefficient of Zijlema (Case 4) shows slightly better statistics than the Powell's (Case 1).Peng & Li's drag formulation (Case 5) shows lower correlation coefficient; higher mean error, Bias, SI, MAE, and E NORM than those of Powell's (Case 1), but it has a higher number of wet stations.It is obvious that the correlation coefficient alone may be a misleading criterion to judge HWM comparison.A simulation may produce less inundation extent, yet score a high correlation coefficient since dry nodes are typically ignored from the statistics.

Effect of Wind Field
SWAN + ADCIRC is used to hindcast Rita's storm surge with different meteorological forcing options, such as HWIND (Case 6), PBL (Case 7), Dynamic Holland Model (Case 8), and Asymmetric Holland Model (Case 9). Figure 8 shows the maximum wind track generated by each wind model and their wind speed difference from that of OWI (Case 1) side by side.All models predicted reduced peak wind speed before the hurricane landfall, correctly indicating the weakened status of Rita before landfall.The OWI wind model captures the hurricane extent reasonably well.The HWIND model (Case 6), seems to do well in the ocean, but shows a somewhat stronger wind field than Case 1.As seen in Figure 8b, HWIND wind speed is higher by 5 m/s at some locations along the track.Note that HWIND does not have snapshots after landfall, and an extrapolation is needed for subsequent snapshots from the last available one.The PBL model (Case 7), overpredicts the intensity by 5 m/s or higher at several locations, especially at landfall (see Figure 8d).The dynamic Holland model (Case 8), captures the maximum winds reasonably well, predicting better than HWIND at some locations.However, HM does not capture the radial extent of the vortex.More specifically, in the middle of the Gulf of Mexico its wind field is underpredicted by 10 m/s or more (see Figure 8f).The Asymmetric Holland model (Case 8) underpredicts the wind speed by about 5 m/s along the track in comparison to OWI.Its radial extent is a matter of concern, although it is better than that of the of Holland model.

Effect of Wind Field
SWAN + ADCIRC is used to hindcast Rita's storm surge with different meteorological forcing options, such as HWIND (Case 6), PBL (Case 7), Dynamic Holland Model (Case 8), and Asymmetric Holland Model (Case 9). Figure 8 shows the maximum wind track generated by each wind model and their wind speed difference from that of OWI (Case 1) side by side.All models predicted reduced peak wind speed before the hurricane landfall, correctly indicating the weakened status of Rita before landfall.The OWI wind model captures the hurricane extent reasonably well.The HWIND model (Case 6), seems to do well in the ocean, but shows a somewhat stronger wind field than Case 1.As seen in Figure 8b, HWIND wind speed is higher by 5 m/s at some locations along the track.Note that HWIND does not have snapshots after landfall, and an extrapolation is needed for subsequent snapshots from the last available one.The PBL model (Case 7), overpredicts the intensity by 5 m/s or higher at several locations, especially at landfall (see Figure 8d).The dynamic Holland model (Case 8), captures the maximum winds reasonably well, predicting better than HWIND at some locations.However, HM does not capture the radial extent of the vortex.More specifically, in the middle of the Gulf of Mexico its wind field is underpredicted by 10 m/s or more (see Figure 8f).The Asymmetric Holland model (Case 8) underpredicts the wind speed by about 5 m/s along the track in comparison to OWI.Its radial extent is a matter of concern, although it is better than that of the of Holland model.Differences in maximum water elevation in comparison to the default case, Case 1, are shown in Figure 9a,c,e,g.Comparative vector plots of wind fields at 8:30 AM on 24 September 2005 from different models are shown in Figure 9b,d,f,h.Case 6 surge generated by HWIND wind field (Figure 9a) shows that the inundation near landfall is within 0.35 m of Case 1, but at the Louisiana coast the surge is underpredicted by 1 m or more.The reason is that HWIND does not have the temporal details of the wind field beyond its snapshot further into the ocean, as seen in Figure 9b.The water elevation of Case 7, as shown in Figure 9c, is at least 0.5 m higher than that of Case 1, which is expected, since the PBL wind field is stronger.The PBL snapshot is large enough to capture the details of Louisiana, as seen in Figure 9d.As shown in Figure 9e, the dynamic Holland model (Case 8), does not predict the storm surge accurately, and underpredicts the surge by 0.5 m near coastal regions.The most likely reason for the underprediction is the symmetric structure of its wind vortex.As can be seen in Figure 8f, the spiral vortex pattern of OWI is conducive to pushing water to land, while the curvature effect of the circular vortex by the Holland model allows more water to go to the east, pushing less water toward the land near landfall.The extent of the wind field is also very limited, and the inundation in the Louisiana area is underpredicted by 1 m or more.As seen in Figure 9g, the Holland Asymmetric model, Case 9, shows a better inundation near landfall than that of the Holland model.The asymmetric wind vortex of the Asymmetric Holland Model, as shown in Figure 9h performs better than the Holland model, but the extent of the wind field is still a concern since the Louisiana area is underpredicted by more than 1 m.Differences in maximum water elevation in comparison to the default case, Case 1, are shown in Figure 9a,c,e,g.Comparative vector plots of wind fields at 8:30 AM on 24 September 2005 from different models are shown in Figure 9b,d,f,h.Case 6 surge generated by HWIND wind field (Figure 9a) shows that the inundation near landfall is within 0.35 m of Case 1, but at the Louisiana coast the surge is underpredicted by 1 m or more.The reason is that HWIND does not have the temporal details of the wind field beyond its snapshot further into the ocean, as seen in Figure 9b.The water elevation of Case 7, as shown in Figure 9c, is at least 0.5 m higher than that of Case 1, which is expected, since the PBL wind field is stronger.The PBL snapshot is large enough to capture the details of Louisiana, as seen in Figure 9d.As shown in Figure 9e, the dynamic Holland model (Case 8), does not predict the storm surge accurately, and underpredicts the surge by 0.5 m near coastal regions.The most likely reason for the underprediction is the symmetric structure of its wind vortex.As can be seen in Figure 8f, the spiral vortex pattern of OWI is conducive to pushing water to land, while the curvature effect of the circular vortex by the Holland model allows more water to go to the east, pushing less water toward the land near landfall.The extent of the wind field is also very limited, and the inundation in the Louisiana area is underpredicted by 1 m or more.As seen in Figure 9g, the Holland Asymmetric model, Case 9, shows a better inundation near landfall than that of the Holland model.The asymmetric wind vortex of the Asymmetric Holland Model, as shown in Figure 9h performs better than the Holland model, but the extent of the wind field is still a concern since the Louisiana area is underpredicted by more than 1 m.performs better than the Holland model, but the extent of the wind field is still a concern since the Louisiana area is underpredicted by more than 1 m.Simulated water level time series with different meteorological forcings in SWAN + ADCIRC are evaluated against observed and published data [1], as shown in Figure 10.As expected, the PBL model (Case 7), overpredicts peaks by at least 1 m in comparison to the observed ones, especially for stations far from the landfall (e.g., J, E, C, and G).Interestingly, the Holland model (Case 8) predicts arrival of the peak a few hours early at almost all stations, which is worse than the results by other models.This peculiarity is probably due to its significantly different wind vortex structure, compared to the rest of the wind models.Simulated water level time series with different meteorological forcings in SWAN + ADCIRC are evaluated against observed and published data [1], as shown in Figure 10.As expected, the PBL model (Case 7), overpredicts peaks by at least 1 m in comparison to the observed ones, especially for stations far from the landfall (e.g., J, E, C, and G).Interestingly, the Holland model (Case 8) predicts arrival of the peak a few hours early at almost all stations, which is worse than the results by other models.This peculiarity is probably due to its significantly different wind vortex structure, compared to the rest of the wind models.Simulated maximum water elevations are compared for each of the cases against measured HWMs [50] (not shown for brevity).The coefficients of correlation are similar in all the cases, and range between 0.66 and 0.71.The HWM statistical comparisons are presented in Table 3. Case 1 has the best correlation coefficient, RMSE, Mean Error, Bias, MAE and ENORM.Although the inundation extents are better for the PBL (Case 7) model, it overpredicts the inundation.Some cases have good statistics, but these have a number of dry stations as well.Based on the reasoning that most HWM locations are actually flooded with reasonable statistical capabilities, OWI, HWIND, Asymmetric Holland, PBL, and Holland models performed best, in that order, among all wind models studied here for the Rita storm surge hindcast.

Combined Effect of Wind Drag and Bottom Coefficient
This study re-examined the earlier findings by Zijlema et al. [15] that reduced drag coefficient with reduced bottom friction performs as well as a large drag coefficient with large bottom friction.As stated earlier in Table 1, Case 12 is similar to Case 1, but does not have a cap of 0.002 for Powell's drag coefficient.Taking Case 12 as the standard, differences of maximum water elevation and Simulated maximum water elevations are compared for each of the cases against measured HWMs [50] (not shown for brevity).The coefficients of correlation are similar in all the cases, and range between 0.66 and 0.71.The HWM statistical comparisons are presented in Table 3. Case 1 has the best correlation coefficient, RMSE, Mean Error, Bias, MAE and E NORM .Although the inundation extents are better for the PBL (Case 7) model, it overpredicts the inundation.Some cases have good statistics, but these have a number of dry stations as well.Based on the reasoning that most HWM locations are actually flooded with reasonable statistical capabilities, OWI, HWIND, Asymmetric Holland, PBL, and Holland models performed best, in that order, among all wind models studied here for the Rita storm surge hindcast.

Combined Effect of Wind Drag and Bottom Coefficient
This study re-examined the earlier findings by Zijlema et al. [15] that reduced drag coefficient with reduced bottom friction performs as well as a large drag coefficient with large bottom friction.As stated earlier in Table 1, Case 12 is similar to Case 1, but does not have a cap of 0.002 for Powell's drag coefficient.Taking Case 12 as the standard, differences of maximum water elevation and velocity between Cases 12 and 10, and Cases 12 and 11 are presented in Figure 11.When both drag coefficient and bottom friction are reduced to 75% of Case 12, as in the Case 10, the maximum water elevation does not differ more than 0.25 m in comparison to that of Case 12 in the shoreline, as seen in Figure 11a.Away from the shoreline the differences are indistinguishable.Similarly, water velocity does not differ by more than 0.25 m/s near the landfall areas as seen in Figure 11b.A very similar conclusion can be made from Figure 11c,d for Case 11, when bottom friction and drag coefficients are increased to 125% of Case 12.These findings are in agreement with Zijlema et al. [15].The water level time series are evaluated against observed and published data [1], and the comparison of some selected stations is shown in Figure 12.Overall, all three cases show similar performance for stations close to the landfall.It is noted that both Cases 10 and 11 predict slightly higher peaks than the observed data at some stations.Since both bottom friction and drag are involved, it is difficult to ascertain the reason for such apparently random behavior.Case 10 predicts early peak arrival compared to other cases for stations away from landfall.Note that it was established earlier that lower bottom friction (e.g., Case 9) causes earlier peak arrival than cases with higher friction (e.g., Case 11).These results indicate a need for further study to determine more accurate bottom friction parameterization, and other physical features that may affect the speed of storm surge away from landfall.
The water level time series are evaluated against observed and published data [1], and the comparison of some selected stations is shown in Figure 12.Overall, all three cases show similar performance for stations close to the landfall.It is noted that both Cases 10 and 11 predict slightly higher peaks than the observed data at some stations.Since both bottom friction and drag are involved, it is difficult to ascertain the reason for such apparently random behavior.Case 10 predicts early peak arrival compared to other cases for stations away from landfall.Note that it was established earlier that lower bottom friction (e.g., Case 9) causes earlier peak arrival than cases with higher friction (e.g., Case 11).These results indicate a need for further study to determine more accurate bottom friction parameterization, and other physical features that may affect the speed of storm surge away from landfall.The peak simulated water levels are compared for each of the simulations with measured HWMs [50] (see Table 3).All simulations are only analyzed for wet HWM locations.Case 11 has overall the best correlation coefficient, RMSE, Mean Error, Bias, MAE and ENORM values, but has 9 dry stations.Note that Case 1 has 7 and Case 10 has 6 dry stations.Interestingly, Case 10 has lower maximum water elevations (see Figure 11), yet it has more wet stations than Case 11.It appears that inundation spread far when lower wind drag and bottom friction are used.

Effect of Cap in Powell's Wind Drag Coefficient
To check the effect of cap in Powell's drag coefficients, differences of maximum water elevation and velocity for Cases 1 and 12 are plotted in Figure 13.The differences are very small, less than 0.05 m, and can be neglected for all practical purposes.At the landfall, Hurricane Rita had stronger wind at the right sector than the left one.The right sector wind pushed the water against the shoreline, which has the potential to increase the water elevation.However, Powell's drag coefficient at the right sector is averaged at 0.00174.The weaker left wind pushed water southward, which was able to freely flow past the shoreline without causing any accumulation.The impact of the high drag coefficient is not fully realized in the Hurricane Rita case.The results could be different if the hurricane structure and geographical features were different at the landfall location.The peak simulated water levels are compared for each of the simulations with measured HWMs [50] (see Table 3).All simulations are only analyzed for wet HWM locations.Case 11 has overall the best correlation coefficient, RMSE, Mean Error, Bias, MAE and E NORM values, but has 9 dry stations.Note that Case 1 has 7 and Case 10 has 6 dry stations.Interestingly, Case 10 has lower maximum water elevations (see Figure 11), yet it has more wet stations than Case 11.It appears that inundation spread far when lower wind drag and bottom friction are used.

Effect of Cap in Powell's Wind Drag Coefficient
To check the effect of cap in Powell's drag coefficients, differences of maximum water elevation and velocity for Cases 1 and 12 are plotted in Figure 13.The differences are very small, less than 0.05 m, and can be neglected for all practical purposes.At the landfall, Hurricane Rita had stronger wind at the right sector than the left one.The right sector wind pushed the water against the shoreline, which has the potential to increase the water elevation.However, Powell's drag coefficient at the right sector is averaged at 0.00174.The weaker left wind pushed water southward, which was able to freely flow past the shoreline without causing any accumulation.The impact of the high drag coefficient is not fully realized in the Hurricane Rita case.The results could be different if the hurricane structure and geographical features were different at the landfall location.underpredicts the surge in Louisiana by 1 m or more in comparison to OWI.Time series plots reveal that the Holland model predicts early arrivals of the peak, most likely due to a significantly different wind vortex pattern.Based on statistical analysis, OWI performs the best, followed by HWIND, Asymmetric Holland, PBL, and Holland models for the hindcast of Rita's storm surge.
When both bottom friction and Powell's wind drag (without cap) are decreased or increased, the overall maximum water elevations pattern did not vary noticeably.Current study further supports the findings of Zijlema that simultaneously decreasing or increasing both bottom friction and wind drag essentially provides almost similar hindcast results.The only noticeable departure is the early arrival (by a few hours) of the surge peak at some stations far from landfall where reduced bottom friction is expected (e.g., Case 10).
The wind drag coefficient of Powell with and without cap did not make much difference in the Rita hindcast results.At landfall, Hurricane Rita had stronger wind at the right sector than the left one.Powell's formulation has high drag coefficient in the left sector, where Rita had a weak southward wind.The weaker left wind pushed water southward, which could freely flow past the shoreline without causing any accumulation.The results could be different if the hurricane structure and geographical features were different at the landfall location.

Figure 1 .
Figure 1.Unstructured ULtralite-Levee-Removed (ULLR) [1] mesh (a) resolution; (b) observation station locations with respect to Hurricane Rita track.Note that the bathymetry is manually set between 1 m and −1 m to show the shoreline.Here a negative bathymetry represents above sea level.

Figure 1 .
Figure1.Unstructured ULtralite-Levee-Removed (ULLR)[1] mesh (a) resolution; (b) observation station locations with respect to Hurricane Rita track.Note that the bathymetry is manually set between 1 m and −1 m to show the shoreline.Here a negative bathymetry represents above sea level.

Figure 2 .
Figure 2. Drag coefficient profiles as functions of wind speed, (a) popular drag formulations used in the present study; (b) 75% and 125% bounds of Powell's [44] sector based drag coefficients.The OceanWeather Inc. (OWI) structured and data-assimilated wind and pressure fields are used as the default wind model, similar to Kerr et al. [1].To study the effect of other meteorological forcing options, the following meteorological setups are used in instead of the OWI model: (1) NOAA National Hurricane Research Division (NHRD) HWIND Model; (2) Planetary Boundary Layer (PBL) Model; (3) Dynamic Holland Model (HM); and (4) Asymmetric Holland Model (AHM).Some technical details of these models are presented below.

Figure 2 .
Figure 2. Drag coefficient profiles as functions of wind speed, (a) popular drag formulations used in the present study; (b) 75% and 125% bounds of Powell's [44] sector based drag coefficients.
Figure 3d displays maximum water elevation, which inundates the coastal areas from Louisiana to Texas.The maximum water velocity is shown in Figure 3e.The highest velocity of 4 m/s is obtained on or near the track.Water flows from the landfall area southward along the Texas coastline.Different aspects of the results are discussed in the following subsections.

Figure 4 . 4 .
Figure 4. Contour plots of maximum water elevation and velocity differences for hindcast of Hurricane Rita storm surges.(a) maximum elevation difference: Case 1-Case 2; (b) maximum velocity Figure 4. Contour plots of maximum water elevation and velocity differences for hindcast of Hurricane Rita storm surges.(a) maximum elevation difference: Case 1-Case 2; (b) maximum velocity difference: Case 1-Case 2; (c) maximum elevation difference: Case 1-Case 3; (d) maximum velocity difference: Case 1-Case 3.

JFigure 7 .
Figure 7. Scatter plots of Rita HWM's.Red, yellow, and light green points indicate over-prediction by the model; blue, dark blue and purple points indicate underprediction.Green points indicate a match within 0.5 m.The black line represents the best fit lines.(a) Case 1; (b) Case 2.

Figure 7 .
Figure 7. Scatter plots of Rita HWM's.Red, yellow, and light green points indicate over-prediction by the model; blue, dark blue and purple points indicate underprediction.Green points indicate a match within 0.5 m.The black line represents the best fit lines.(a) Case 1; (b) Case 2.

Figure 7 .
Figure 7. Scatter plots of Rita HWM's.Red, yellow, and light green points indicate over-prediction by the model; blue, dark blue and purple points indicate underprediction.Green points indicate a match within 0.5 m.The black line represents the best fit lines.(a) Case 1; (b) Case 2.

Figure 8 .
Figure 8. Rita maximum wind track plots (i.e., 'max wind') and difference of wind speed between OWI and other wind models (i.e., 'max wind diff') (a and b) Case 6; (c and d) Case 7; (e and f) Case 8; (g and h) Case 9.

Figure 8 .
Figure 8. Rita maximum wind track plots (i.e., 'max wind') and difference of wind speed between OWI and other wind models (i.e., 'max wind diff') (a and b) Case 6; (c and d) Case 7; (e and f) Case 8; (g and h) Case 9.