Storm surge and wave simulations in the Gulf of Mexico using a consistent drag relation for atmospheric and storm surge models

To simulate winds and water levels, numerical weather prediction (NWP) and storm surge models generally use the traditional bulk relation for wind stress, which is characterized by a wind drag coefficient. A still commonly used drag coefficient in those models, some of them were developed in the past, is based on a relation, according to which the magnitude of the coefficient is either constant or increases monotonically with increasing surface wind speed (Bender, 2007; Kim et al., 2008; Kohno and Higaki, 2006). The NWP and surge models are often tuned independently from each other in order to obtain good results. Observations have indicated that the magnitude of the drag coefficient levels off at a wind speed of about 30 m s −1, and then decreases with further increase of the wind speed. Above a wind speed of approximately 30 m s −1, the stress above the air-sea interface starts to saturate. To represent the reducing and levelling off of the drag coefficient, the original Charnock drag formulation has been extended with a correction term. In line with the above, the Delft3D storm surge model is tested using both Charnock’s and improved Makin’s wind drag parameterization to evaluate the improvements on the storm surge model results, with and without inclusion of the wave effects. The effect of waves on storm surge is included by simultaneously simulating waves with the SWAN model on identical model grids in a coupled mode. However, the results presented here will focus on the storm surge results that include the wave effects. The runs were carried out in the Gulf of Mexico for Katrina and Ivan hurricane events. The storm surge model was initially forced with Hwind data (Powell et al., 2010) to test the effect of the Makin’s wind drag parameterization on the storm surge model separately. The computed wind, water levels and waves are subsequently compared with observation data. Based on the good results obtained, we conclude that, for a good reproduction of the storm surges under hurricane conditions, Makin’s new drag parameterization is favourable above the traditional Charnock relation. Furthermore, we are encouraged by these results to continue the studies and establish the effect of improved Makin’s wind drag parameterization in the wave model. The results from this study will be used to evaluate the relevance of extending the present towards implementation of a similar wind drag parameterization in the SWAN wave model, in line with our aim to apply a consistent wind drag formulation throughout the entire storm surge modelling approach.


Introduction
Storm surge is largely determined by a combination of inverse barometric effects and air-sea interactions.Under hurricane conditions, the transfer of momentum from the atmosphere into the ocean represents the dominant process.This transfer of momentum results in the generation of the surge and waves which then interact.The generated wave field introduces radiation stress gradients which affect the water level (and currents).
For the computation of the surface momentum flux, traditionally NWP and storm surge models use a common bulk relation, which involves a formulation for the drag coefficient C D .The magnitude of this coefficient increases to a good approximation linearly with increasing wind speed (Charnock, 1955;Large et al., 1981;Smith et al., 1975).Often, the drag coefficient is capped at low and high ranges of wind speed (see Fig. 1).This type of wind drag parameterization has proven to work well for simulations of storm surge under extra-tropical weather conditions (de Vries et al., 1995;Gerritsen et al., 1995), albeit after applying some fine tuning.
The validity of Charnock's relation, or other similar wind drag relation, has limitation when applied to hurricane conditions.Observations and laboratory experiments (Powell et al., 2003;Donelan, 2004;Jarosz et al., 2007) indicate that for hurricane wind speeds, the momentum flux starts to saturate.At a wind speed of ∼30 m s −1 , the magnitude of the drag coefficient reaches its peak value and it decreases when the wind speed further increases.Makin (2005), Bye and Jenkins (2006) and Kudryavtsev (2006), among others, have come up with an explanation for this drag reduction.By describing the impact of spray droplets on the atmospheric flow in the wave boundary layer, Makin (2005) was able to explain the observed reduction in the surface drag and proposed a drag parameterization that reflects this reduction.To test the implementation of the improved drag parameterization in the NWP HIRLAM model and the Delft3D storm surge model in a consistent manner, a number of numerical experiments are carried out.Due to frequent occurrences of hurricanes and availability of wind and observational data, these experiments were performed in the Gulf of Mexico.The HIRLAM winds for hurricanes Ivan (2004) and Katrina (2005) were simulated on a grid with a horizontal resolution of 0.05 degrees in order to resolve the strong gradients in the sea level pressure and large (local) gradients in the hurricane wind field.The HIRLAM model and the impact of the wind drag parameterization on the simulated hurricane intensity and wind field have been described in detail by Zweers et al. (2010) and it will not be repeated here.They showed that hurricanes cannot attain their intensity in fore-casts with the Charnock relation and that the new parameterization improved the simulated hurricane winds significantly (see Fig. 2).
This paper focuses on the results of the coupled storm surge and wave models using HIRLAM wind data that was produced using the Charnock relation as well as the improved drag parameterization by Makin.The wave and storm surge models are subsequently forced with these data.To check the behaviour of wind drag parameterization separately on computed surge levels, H * wind data (Powell et al., 1998(Powell et al., , 2010) ) for hurricanes Ivan and Katrina have been used to force the coupled storm surge-wave models using both Charnock's and Makin's drag formulation.The H * wind data used will be described and discussed in detail later in this article.By combining the different wind forcing and drag parameterisations we have simulated four scenarios: -H * wind data using Charnock's wind drag parameterization in the storm surge model; -H * wind data using Makin's wind drag parameterization in the storm surge model; -HIRLAM wind data (simulated using Charnock's parameterization, will be denoted as HIRLAM1) and Charnock's wind drag parameterization in the storm surge model; and -HIRLAM wind data (simulated using Makin's parameterization, will be denoted as HIRLAM2) and Makin's wind drag parameterization in the storm surge model.

The model grids and boundary condition
The coupled storm surge and wave simulations are performed on two different model grids and schematisations, defined on a spheroid, in an off-line nested configuration (see Fig. 3): -  (Egbert and Erofeeva, 2002).Due to the size of the GoM model, tide-generating   (Booij, 1999) uses identical grids as the storm surge model.

Model schematisation and bathymetry
The topography and the bathymetry of the Panhandle model (see Fig. 3) were constructed in 2008 using data downloaded from the National Geophysical Data Centre (NGDC) and United States Geological Survey (USGS) web sites.The quality and the resolution of the topography and bathymetry data have, since then, been improved.However, given the relative coarseness of the Panhandle model and given the aim of the present study is to gain insight in the overall behaviour of the wind drag parameterizations on the induced surges, it was deemed not essential to update the models with new data as of yet.

Important model input parameters
The Delft3D storm surge models are run with an integration time of respectively 10 min for the GoM and 2 min for the higher resolution Panhandle model.The bottom friction is computed using the Manning's formulation with the value of the coefficient equal to 0.024 on sea and 0.035 on land. of Hasselmann (1974) or Van der Westhuysen et al. (2007).
Depth-induced breaking is computed through a model by Battjes and Janssen (1978) with the breaking index γ = 0.73.The bottom friction is based on the JONSWAP formulation (Hasselmann et al., 1973) with the friction coefficient set equal to 0.067 m 2 s −3 .The effect of wave induced forces in the storm surge model is taken into the account by coupling SWAN wave to the Delft3D storm surge model.The wave forces will, among others, enhance the energy dissipation near the bottom in the storm surge model and generate a net mass flux affecting the current, especially in the cross-shore direction.These effects are accounted for by passing on radiation stress gradient determined from the computed wave parameters every 10 min to the storm surge model during the next 10 min.The water levels and currents computed by the storm surge model are then passed on to the wave model after this interval, which is then used by the wave model to compute the wave parameters.The choice of this interval of 10 min has been made after experimenting with different coupling times and was found to be optimal in view of the required computation time.

Monitoring stations
For comparison with observations, water level time series are stored and compared at a number of stations.In this paper, the results at four representative coastal stations are presented.Similarly, the time series for the wave parameters have been compared at a number of on-shore and offshore stations.Here, the results are presented for four representative coastal and offshore stations.An overview of the station locations is given in Fig. 4.

Wind data
As mentioned earlier, the implementation of Makin's improved wind drag parameterization in Delft3D is initially tested by applying the 3 hourly gridded H * wind data to the coupled storm surge and wave models after transforming the data into a moving circular grid (Vatvani et al., 2002).H * wind data, which is constructed through analysis of observation data from multitude sources and defined as winds prevailing at 10 m above surface level (U 10 ), has often been applied to simulate storm surges (Dietrich et al., 2011).The originally 1 min averaged H * wind data has been adjusted to 10 min averaged wind speed, a time scale more appropriate for determining the ocean response and commonly applied for storm surge models, using the recommended conversion factor mentioned in the WMO report (Harper, 2010).However, the data is provided on a moving grid, centred on the hurricane eye, covering a rectangle area with a size of approximately 1000 by 1000 km.If the coverage of the data is not extended, the wind speed in almost all stations would have been equal to zero when the station concerned is outside the domain of the original H * wind data.This effectively introduces a discontinuity in the wind field, which is physically unrealistic, and would cause the SWAN model to generate unrealistic wave fields with long periods (>10 s) along the discontinuous wind field front.Initially, those wave heights are small, which then propagate through the entire model area.However, these spurious waves can grow in height when the actual wind intensifies.To prevent the generation of these unrealistic waves, the data gap outside the area covered by the original H * wind data has been merged with synthetically computed hurricane winds covering a larger area using Holland's method (1980).As we are not interested in generating the entire hurricane wind field and want only to supplement the data gap far away from the hurricane centre, the synthetic hurricane winds have been generated using the following simple but commonly applied values for the Hollands' parameter A and B: Where V max equals the reported maximum wind speed, ρ equals the air density, P drop equals the reported pressure drop and R max equals the radius of maximum wind speed.In our computation, R max has been set to 25 km.The generated synthetic winds have been adjusted to account for the interaction with the steering flow (Chan and Gray, 1982) inducing the asymmetry in the wind field.
The wind speed time series of the extended and the original wind field at station W42003 are shown, as an example, in Fig. 5.In the same figure, observed wind speed is also presented.The 10 min averaged observed wind speed (anemometer height equals 5 m) has been converted to 10 m above surface level by applying a height correction factor using the Power Law Method. 2 The figure shows good agreement between the extended data and observed wind speed for both hurricane events at this station.Similar improvement is found for other stations (W42001, W42002, W42007 and W42040).When these extended wind fields are applied, the spurious wave fields are not generated anymore.The HIRLAM winds, which are already defined as 10 min averaged wind and as winds prevailing at 10 m above surface level (U 10 ), are then subsequently applied to the models without changing other model input parameters.The results for the storm surge and wave computations for H * wind and HIRLAM wind are discussed below.

Results
In the following sections, the computed and observed surge and wave parameters will be presented and discussed for both hurricane events.

Storm surge
The computed surge levels at different stations for both hurricanes and for all the scenarios are compared to the observed surge levels in Fig. 6.Also presented is a table containing quantitative analysis of the computed surge in the form of RMS error (Table 1).Besides the time series of the computed surge levels, the time series of the contribution of the wave induced surge is also depicted.The wave induced surge is computed by running the same model suite with identical wind scenarios but without taking into account the effect of wave and wave-surge interaction and subtracting these results from the coupled surge-wave model results.
With respect to the drag parameterization, we observe that Charnock's parameterization tends to overestimate the surge levels considerably, especially when the computed surge level is large (>1 m), which can be directly associated with periods of strong wind conditions (wind speed >30 m s −1 ).This can be expected as the values of the wind drag coefficient based on Charnock's parameterization is significantly larger in this high wind speed regime (see Fig. 1).The over-estimation also still occurs when HIRLAM1 winds are applied, despite lower computed wind speeds.When the surge levels computed using Makin's drag formulation are compared with those using Charnock's drag relation, then it can be observed that the average RMS errors for the stations considered, especially in the Katrina case, have been reduced by 10 to 30 %.With the exception of station Pensacola in the Katrina case, where the RMS error is increased, overall we can conclude that a storm surge model produces better results when run with Makin's drag formulation.
For scenarios where the surge compares best with data, the contribution of the wave towards the setup in water level can reach up to 30 % of the total surge.Directly on the coast line, on very shallow areas, this contribution may occasionally even be higher, reaching a setup value as high as 1.50 m, as shown in Fig. 7.This finding is consistent with results from Resio and Westerink (2008), Wolf (2008) and Dietrich et al. (2011).This shows that wave contribution on surge can be significant and should not be neglected.
Comparing the different wind sources used, we see that HIRLAM2 winds for Katrina event is shown to be capable of reproducing the storm surge during for both events with almost similar quality (RMS error = 0.19 cm) as the H * wind data (RMS error = 0.17 cm).The surge computed at station Grand Isle East is an exceptional case where according to Zweers et al. (2010) the track computed by HIRLAM is slightly off the actual track.In the Ivan case, HIRLAM2 winds did not perform as well which can already be expected when we compare the computed and observed maximum winds depicted in Fig. 2.

Wave results
The observed and computed significant wave heights, wave direction and mean wave period for both hurricane events for all the scenarios at a number of stations are depicted in Fig. 8.The RMS errors at eight different station for all scenarios are presented in Table 2. Since the drag parameterization in the wave model has not been changed and the winds used in these experiments are comparable in speed and direction, the effect of different drag parameterization in the storm surge model on the computed waves is expected to manifests primarily through flow and water levels computed by the storm surge model.Because water level attribute strongest to the wave conditions in relative shallow areas, their potential effects on the waves generated are expected to be the largest in shallow areas and only when the computed water levels (relative to the depth) for the scenarios simulated differs substantially.When the differences in the computed water levels are small, then the simulated wave parameters will be comparable, which is confirmed by our computational results (see Fig. 8 and Table 2).The computed wave heights, wave directions and wave periods, except for local differences to be discussed below, do not deviate strongly from each other.This is also confirmed by the computed average RMS errors are almost equal for all the scenarios computed.
For Ivan, the computed significant wave heights at all stations resembles the observed data reasonably well.However, at station CSI06 there is a time lag of approximately 12 h before the significant wave height increases to its peak value.At stations W420003, W420039 and W420040 the computed wave direction follows the trends of the observed data, but the absolute differences can be as large as 90 degrees (at W42003 around the 16th of September).The computed mean wave period resembles the observed data well but again at CSI-06 the wave period shows similar time lag of approximately 12 h.In the case of Katrina, the computed wave parameters are significantly better.However, even here, we still see some discrepancies.At stations W42007 and CSI-06 the computed mean wave period is underestimated by approximately 4 s.
Comparison of the computed spatial pattern of the significant wave height and wave mean period using H * wind data and Makin's drag formulation (Fig. 9) with modelled results from Dietrich et al. (2011) for Katrina at 29th August 2005 10:00 UTC (depicted in Fig. 8 in the referred article) shows that, although along the coast, his results for the significant wave height values as well as the mean wave period values are slightly higher (i.e.wave penetrates slightly further towards the coast); the overall results are comparable, i.e. modelled significant wave height south, southeast and east of Birdfoot area equals approximately 21 m fanning out to a value of 9 m and the computed wave period values modelled in the same area lying between 12-15 s.Therefore, we believe that despite the local discrepancies encountered above, the modelled wave results are sufficient to represent the overall effect of wave on storm surge with reasonable accuracy.

Conclusions
This paper has presented the coupled storm surge and wave model results for a combination of four different wind and drag parameterization scenarios.When Charnock's wind drag parameterization is applied in the storm surge model, it consistently overestimates the storm surges for high wind speed values.This result is expected for the H * wind data, but the overestimation of the surge also occurs when the HIRLAM1 winds is applied, despite lower wind speed computed.
Based on quantitative and qualitative analysis of the results, we can conclude that by applying Makin's improved drag formulation to the storm surge model, the computed surge levels compare the observed surge levels very well, especially those simulated with H * wind data.However, it should be stressed that the coverage of H * wind data should be extended in case they are applied over a larger model area to avoid discontinuity in the wind field, which is physically unrealistic and would cause the SWAN model to generate unrealistic wave fields.
The influence of wave on storm surge levels results are significant and should not be neglected in storm surge forecast or hindcast.Results indicate that the choice of wind source (and wind drag coefficient used in the storm surge model) is less important, as the set-up generated by wave is somewhat less sensitive to the wind sources applied.
Finally, these results provide some confidence in our effort to apply a consistent drag parameterization throughout the entire modelling suite for storm surge simulations.In our view, adjusting the wind drag parameterization in the SWAN model and direct coupling the NWP models with storm surge and wave models would be a logical follow-up of this study.

Fig. 1 .
Fig. 1.Drag coefficient C D10 at 10 m height as a function of the 10-m wind speed U 10 according to Charnock's relation using the Charnock's constant value equal to 0.025 (solid, red) and the new parameterization following Makin (2005) (solid, blue) from Zweers et al. (2010).Dashed lines represent the approximations used in Delft3D.Observational data by Powell et al. (2003) are indicated by diamonds.

Fig. 3 .
Fig. 3.The Gulf of Mexico and Panhandle model grids (top) and bathymetry of the Panhandle model (bottom); open boundaries of the models are shown as white lines.forces inside this model are switched on.The tidal representations of these models have been validated with satisfactory results.As no wave forcing has been applied to the open boundary of the GoM model, all waves are generated inside the coarse grid model by the hurricane winds.The wave information is then prescribed along the open boundaries of the Panhandle model.The SWAN wave model 1(Booij, 1999) uses identical grids as the storm surge model.

Fig. 4 .
Fig. 4. Overview of monitoring stations where time series of the computed water levels and wave parameters are compared to observed data.

Fig. 5 .
Fig. 5. Extended H * wind data at buoy location W42003 compared to the original H * wind data and 10 min averaged observed data at 10 m above surface level.Left: results for hurricane Ivan and Right: results for hurricane Katrina.For the location of W42003, see Fig. 4. Discontinued gray line implies the gauge stopped recording the data.Discontinued crosses implies that the station lies outside the coverage of the original H * wind gridded data.

Fig. 6a .
Fig. 6a.Computed and observed surge level at a number of stations for hurricane Ivan (left) and Katrina (right).In each pane the contribution of the wave on the surge is depicted separately.For the location of stations see Fig. 4. Discontinued black line implies the gauge stopped recording the data.

Fig. 7 .
Fig. 7. Computed wave induced surge for Ivan and Katrina at two distinct times.

Fig. 8a .
Fig. 8a.Computed wave parameters (significant wave height, wave direction and mean wave period) compared to the observed data at a number of stations for hurricane Ivan (left) and Katrina (right).For the location of stations, see Fig. 4. Discontinued dots imply the gauge stopped recording the data.

Fig. 9 .
Fig. 9. Hurricane Katrina H * wind data and computed wave parameters along the Panhandle coastal area at 29th of August 2005 10:00 UTC.Top: wind contours and vectors, one minute averaged.Middle: significant wave height contours and wind vectors.Bottom: mean wave period contours and wind vectors.The high value of computed mean wave period in a small area around point 269.5 E, 29.25 N is a numerical artefact.

Table 1 .
RMS Error of the computed surge levels (wave effect included) at four stations depicted in Fig.6.Empty column implies lack of sufficient observed data.

Table 2 .
RMS Error of the computed significant wave height and average wave period at 10 stations depicted in Fig.8.Empty column implies lack of sufficient observed data.