Relative accuracy of HWRF reanalysis and a parametric wind model during the landfall of Hurricane Florence and the impacts on storm surge simulations

Prediction and reanalysis of storm surge rely on wind and pressure fields from either parametric tropical cyclone wind models or numerical weather model reanalysis, and both are subject to large errors during landfall. This study assesses two sets of wind/pressure fields for Hurricane Florence that made landfall along the Carolinas in September 2018 and appraises the impacts of differential structural errors in the two suites of modeled wind fields on the predictive accuracy of storm surge driven thereby. The first set was produced using Holland 2010 (H10), and the second set is the Hurricane Weather Research and Forecasting (HWRF) reanalysis created by the NWS National Centers for Environmental Prediction (NCEP). Each is validated using a large surface data set collected at public and commercial platforms and then is used as input forcing to a 2-D coastal hydrodynamic model (Delft3D Flexible Mesh) to produce storm surge along the Carolina coasts and major sounds. Major findings include the following. First, wind fields from HWRF are overall more accurate than those based on H10 for the periphery of the storm, though they exhibit limitations in resolving high wind speeds near the center. Second, applying H10 to the best track data for Florence yields an erroneously spike in wind speed on September 15th when the storm reduced to a tropical depression. Third, HWRF wind fields exhibit a progressively negative bias after landfall, likely due to deficiencies of the model in representing boundary layer processes, and to the lack of assimilation of surface product after landfall for compensating for these deficiencies. Fourth, using HWRF reanalysis as the forcings to Delft3D yields more accurate peak surges simulations, though there is severe underestimation of surge along the shoreline close to the track center. The peak surge simulations by Delft3D are biased low when driven by H10, even though over several locations the H10 model clearly overpredicts surface wind speeds. This contrast highlights the importance of resolving wind fields further away from the center in order to accurately reproduce storm surge and associated coastal flooding.


Introduction
Flooding in the coastal areas caused by rainfall and storm surges imposes devastating effects on lives, environment as well as economy for nations across the globe. Improving the resilience of coastal communities requires accurate predictions of coastal flooding in general and storm surge in particular. Despite the recent advances in modeling and computational techniques, accurate representation and prediction of storm surge remain challenging. Hydrodynamic models are often unable to adequately reproduce patterns of coastal flooding both in space and time. Reasons are manifold. These include uncertainties in atmospheric forcing (Cardone and Cox 2009;Dietrich et al. 2018;Mayo and Lin 2019;Abdolali et al. 2021); a lack of available high-quality and high-resolution bathymetric data in shallow areas as well as insufficient grid resolution to resolve the topo-bathy in complex coastal areas (Hell et al. 2012;Jacob and Stanev 2021;Acosta-Morel et al. 2021); a lack of available high-quality water level records during past tropical storms for model calibration (Asher et al. 2019); under-or misrepresentations of physical processes pertaining to flow and to atmosphere-ocean momentum exchange under strong wind regime and in shallow waters (Olabarrieta et al. 2012), a lack of coupling of different model components which have a direct effect on each other (i.e. atmospheric forcing, storm surge, wave and hydrology; Ma et al. 2020;Santiago-Collazo et al. 2019). Among these, abrupt changes in wind fields during landfall due to the increased drag have been recognized as a major source of errors in the simulation and prediction of wind fields by numeric weather models (Wang 2012;Leroux et al. 2018;Kim et al. 2020), and these errors have been cited as a key impediment to accurate storm surge simulations for locations along estuaries and coastal streams further away from the coast (Ferreira et al. 2014).
To date, post-analysis of storm surge and coastal risk assessment have often relied on parametric models of tropical cyclone (TC) wind and pressure fields constructed using semi-empirical relations (Mattocks and Forbes 2008;Vickery et al. 2009;Forbes et al. 2010). These models incorporate idealized assumptions of tropical storm structures and boundary layer processes (Holland 1980;Demaria et al. 1992;Houston and Powell 1994;Vickery et al. 2000;Phadke et al. 2003;Willoughby et al. 2006;Zhang et al. 2011;Chavas et al. 2015;Kepert et al. 2016), require only a few storm parameters and therefore simple to implement, and have the advantages of being able to directly integrate forecasted and observed TC track data. Among the contemporary parametric models, perhaps the most well-known is the model by Holland (1980) and its later variants (i.e., Holland 2008 andHolland et al. 2010), which now serve as the default forcing mechanism for coastal hydrodynamic models such as ADCIRC (Luettich et al. 1992) and Delft3D (Deltares 2014).
While these parametric wind models have been widely applied and, in many cases, yielded satisfactory storm surge simulations, a number of concerns remain. It is well known that these models are limited in their ability to resolve inter-storm variations in wind profiles, asymmetry in the storm structures, and changes of TC structures induced by enhanced friction drags during landfall (MacAfee and Pearson 2006;Fang et al. 2020). In addition, parametric models lack the ability to resolve background wind and pressure fields. These impair their ability to capture the waves, swells and surges generated by wind and pressure at longer distance from the TC center (Abdolali et al. 2021). Though practitioners often resort to the calibration of hydrodynamic models as a countermeasure to compensate for biases and errors in forcings as well as deficiencies in model structures (Lin and Chavas 2012), the efficacy of this practice, however, remains questionable.
Over recent years, high-resolution numerical weather prediction models (NWP) realtime or retrospective analyses of past landfalling TCs have become widely available, and these products have been increasingly applied in predicting/reconstructing storm surge events. Notable extant real-time analysis and reanalysis products include those from the Hurricane Weather Research and Forecasting (HWRF) model that is operational at the US National Weather Service (NWS; Ma et al. 2020), Hurricanes in a Multi-Scale Oceancoupled Non-hydrostatic model (HMON) running operationally at NCEP, GFDL model, and HiRES from European Center for Medium-range Weather Forecast (ECMWF; Molteni et al. 1996). These NWP models incorporate explicit representations of states of atmosphere and their interactions with ocean and land, and have the ability to assimilate a variety of surface and remotely sensed observations. HWRF forecast, for example, is produced by assimilating Doppler velocity from ground-based or air-borne Doppler radar (Tong et al. 2018;Lu and Wang 2020;Davis et al. 2021); upwelling microwave radiation measured by an airborne radiometer , and dropsonde observations (Powell et al. 2003;Franklin et al. 2003;Ryan et al. 2019). Owing to these strengths, these products are expected to offer physically more realistic depictions of TC wind and pressure fields during landfall than do those from parametric models. Nonetheless, the NWP models themselves are subject to biases and errors that arise from mis-or underrepresentation of processes. In particular, there have been reports of large departures of HWRF wind analysis and prediction of TCs after landfall (Kloetzke 2019;Ma et al. 2020), which likely reflect inadequate representations of boundary layer processes.
Heretofore, a plethora of studies have been undertaken with the purpose of illuminating the evolution of TC wind structures Wang 2012) over ocean and during landfall, assessing the skills (Resio et al. 2017;Annane et al. 2018) of NWP models in prognosticating TC tracks and structures, and predicting storm surge (Leroux et al. 2018;Bucci et al. 2021). Yet, very few of these have attempted to appraise the relative realism of wind and pressure fields produced by parametric models versus those from parametric models prior to, during and after landfall, or to assess the impacts of structural errors in these products on the surge simulations driven thereby. A notable exception is Dietrich et al., (2018), in which the authors examined the storm surge forecasts forced by a parametric model and predictions from a WRF model for Hurricane Isaac. The study, however, did not delve into the differential structural errors in the wind profiles that led to the contrasting predictive accuracy of ADCIRC.
The present study is motivated by the need of determining the relative strengths as well as deficiencies of parametric wind model and NWP reanalysis as forcings for storm surge simulations; it offers detailed, comparative analyses of wind and pressure fields from HWRF and the Holland (2010) wind model (henceforth referred to as H10) for Hurricane Florence, and the storm surge simulations driven by each. The specific objectives of the present study are twofold. The first is to determine the relative skill of H10 versus HWRF in reproducing the evolution of wind and pressure fields of Florence during and after its landfall, with a focus on identifying distance and quadrant dependent errors that can be related to the interactions between the storm and land. The second is to gauge the relative efficacy of the two forcing data sets in reproducing the inundation processes along the coast and major sounds. Hurricane Florence was chosen for the study for three reasons: 1) it is one of the most devastating landfall Hurricanes along the Carolinas in recent history; 2) it produces storm surge that penetrated far upstream (> 50 km) and speed and direction of wind over land may have strongly modulated the intensity of flooding along major sounds, and 3) a high-resolution HWRF reanalysis is available, so are a rich set of surface wind observations for validating the wind fields. Our working hypotheses include: a) accuracy 1 3 of wind fields from both H10 and HWRF reanalysis would deteriorate after landfall; b) quality of the H10 wind fields would decline more drastically due to its lack of explicit accounting for the increase in friction drags; and c) using HWRF reanalysis would yield more accurate storm surge simulations throughout the event.
The remainder of this paper is structured as follows. Section 2 provides descriptions of the methods, including the wind forcing (i.e., H10 and HWRF), and selected hydrodynamic model (i.e., Delft3D-FM), model parameters and inputs, and validation matrices. The objectives of this study are accomplished in Sect. 3 (i.e., Results), where the relative accuracy of two wind forcings is assessed and the impact of the storm surge simulations is determined. Section 4 discusses the findings of the study in details, and Sect. 5 summarizes the findings and offers recommendations for future works.

Hurricane Florence
According to NOAA's National Centers for Environmental Information (NCEI) and National Hurricane Center (NHC), Hurricane Florence of 2018 was by far one of the costliest (i.e., 12th) hurricanes that hit the mid-Atlantic region in recent history. The hurricane was first spotted as a tropical disturbance near Cape Verde Island off the West Africa coast in late August experiencing rapid intensification in early September and became a Category 4 hurricane around 5 September. The strength of the storm then declined to a tropical storm while traversing the Atlantic Ocean until 11 September, when it intensified again to a Category 1 hurricane. Florence made its first landfall on the south of Wrightsville Beach, NC on 14 September while retaining Category 1 strength, even though wind speed was mostly below 70 mph while approaching the coast (Fig. 1). The storm produced sizable storm surges across the NC coast which penetrated as far as 50-60 km inland through major rivers. Along the Neuse River near the city of New Bern, the storm surge exceeded 3-m and caused widespread flooding across the city that made the national headline. The major devastations by the storm were caused by its heavy rainfall-the maximum storm totals (up to 17 September) exceeded 900 mm in NC, shattering the previous record for the state. The storm degenerated into a depression on 15 September and underwent extratropical transition on 17 September before dissipating on the 18th, though flooding in many locations lingered on till the end of the month (Stewart and and Berg 2019).

Reanalysis product based on observation
The Hurricane Weather Research and Forecast (HWRF) system tropical-storm predictions that was by a consortium using the WRF Non-hydrostatic Mesoscale Model (WRF NMM; Janjic 2004) core maintained by the NWS National Centers for Environmental Prediction (NCEP). First introduced to operation in 2007, it has become one of the primary forecasting models for TC predictions in the NWS. The HWRF uses telescopic nesting: in its current operational setting, the parent domain of HWRF spans approximately 77.2° × 77.2° with a ~ 27 km mesh, the intermediate domain ~ 17.8° × ~ 17.8° with a ~ 9 km mesh, and the innermost domain ~ 5.9° × ~ 5.9° with ~ 1.5 km mesh (Biswas et al. 2018). The latest version has 75 vertical levels with 10 hPa increment. HWRF runs are coupled with Princeton Ocean Model (MPIPOM-TC) for all oceanic basins in the northern hemisphere. The HWRF model uses forecasts from the Global Forecast System (GFS) on the parent domain as lateral boundary conditions. In the forecast mode, it relies on a synthetic vortex to initialize a TC forecast (Biswas et al. 2018) and a bogus vortex to cold-start strong storms. The bogus vortex is created by smoothing the wind profile of the 2-D vortex until its radial maximum wind (RMW) matches the observed values. The HWRF uses a vortex relocation procedure in which a 6-h HWRF forecast from the previous cycle is used in determining the location of vortex (Nadimpalli et al. 2021), and which infuses position, structure, and intensity from the National Hurricane Center (NHC) storm message (Leslie and Holland 1995;Kwon and Cheong 2010;Zou et al. 2015;Biswas et al. 2018;Zhang et al. 2021). After identifying the vortex location, the initial conditions are further refined by assimilating a variety of observations including those by Stepped Frequency Microwave Radiometer (SFMR; Uhlhorn and Nolan 2012), airborne Doppler radar units, and dropsonde (Biswas et al. 2018).
In this study, we acquired the HWRF reanalysis wind and pressure fields for Hurricane Florence (2018) created at the NCEP EMC that are available from 0600 UTC on September 09, 2018, through 1200 UTC on September 18, 2018, or 9 days. The production of the reanalysis involved retrospectively running the HWRF using the best track from NHC and assimilating the surface and remotely sensed products. The HWRF wind and pressure reanalysis are on a grid mesh of approximately 1.5 km at 1-h intervals.

Parametric Wind Model
In this study, we adopted the Holland et al. (2010; referred to henceforth as H10) parametric pressure and wind model. H10 evolved from the original parametric model by Holland (1980, referred to henceforth as H80). In comparison with other contemporary analytical models (e.g., Emanuel 2010; Emanuel and Rotunno 2011;Chavas et al. 2015), it has the advantages of structural simplicity, relying on few assumptions regarding the structure of the hurricane boundary layer, and being widely used and tested. Holland et al. (2010), for example, demonstrated that H10 outperforms Emanuel (2010) model in capturing the surface wind profile as detected by aircraft reconnaissance. Lu et al. (2018) found that a parametric rainfall model based on H80, the predecessor of H10, better represents the rainfall fields for Hurricanes Isabel and Irene.
A brief review of the structure and evolution of the Holland models is provided here. The first Holland model, the H80, was formulated with the assumption that wind speed is invariant with height in the boundary layer (or the so-called slab model ;Holland 1980) and the surface pressure profile is approximated by a rectangular hyperbola (Schloemer 1954). The model relies on the assumption of cyclostrophic balance to estimate the wind speed close to the eye, and of gradient wind balance to derive the wind-pressure relations further away from the center. H80 uses two parameters, namely a, the scale parameter, and b, the shape parameter, and it requires as input center and environmental pressure, maximum wind speed, and the radius of maximum wind speed, Holland (2008) replaces the fixed b parameter by a time-variant b s that is estimated either from observed surface pressure and temperature when such observations are available, or alternatively from pressure drop at the center, change in pressure, translation speed and latitude. The resulting model, referred to as H08, produces wind profile directly at the surface without resorting to boundary-layer reduction relations (though in H80 the b parameter can be tuned to yield surface wind).
H08 further refines the model by introducing a radially variable exponent x to replace the constant ½ that arises from the cyclostrophic balance.
In the wind model of Holland (2010), herein denoted as H10, surface wind speed takes the following forms, where ΔP s is the pressure drop from a defined external pressure p ns to the cyclone center P cs , s is the surface air density, and e is the base of natural logarithms. The exponent b is a scaling parameter that defines the proportion of pressure gradient near the maximum wind radius. The equation can be rewritten as follows: where the subscript s refers to surface values at a nominal height of 10 m, v ms denotes maximum wind speed (Vmax); r v m is the radius of maximum wind speed (RMW), and x is a scaling parameter that adjusts the profile shape. Parameter b s is related to the original b by b s = bg s x , where g s is the reduction factor for gradient-to-surface winds. The b s parameter is estimated by following Holland (1980) using, In the absence of surface observations of pressure and temperature, b s for a given radius range is expressed as a function of incremental variation in pressure along the radius, temporal change in central pressure, translation velocity and latitude: where Δp s in hPa, P cs t is the intensity change in hPah −1 ; φ is the absolute value of latitude in degrees; and v t is the translation speed of cyclone in ms −1 .
The exponent x is related to incremental change in pressure.
And, the maximum wind speed is determined by surface pressure depression, vapor pressure and the revised Holland b s parameter: H10 recommends blending in a secondary wind maximum that has often been observed in major hurricanes (Willoughby et al. 1982;Wunsch and Didlake 2018). However, the authors caution that a large perturbation may result in a change of vorticity gradient and lead to barotropic instability (Holland et al. 2010). In order to avoid this instability, and out of the concern of uncertainties associated with respect to the locations and magnitude of the maxima, we chose not to implement the secondary maximum; rather, we adopted the central pressure-maximum wind relationship in this study that was described in Holland (2008), and our implementation uses radius for wind speeds of 34, 50, 64 and 100kts that are provided in the NHC Best Track (i.e., HURDAT2; Landsea et al. 2004).

Storm surge model
The Delft3D Flexible Mesh, also known as Delft3D-FM, is a fully integrated software suite by Deltares (2014), consisting of models that simulate a variety of processes in ocean, estuarine, tidal and inland streams, including those related to flow, sediment transport, morphodynamics, and water quality. Its hydrodynamic model, Delft3D-FLOW, is a finite-element model that uses unstructured (triangular) grid mesh. Delft3D-FLOW allows for both 2-and 3-dimensional representations of flows. The 3-D, depthresolving, version solves the full incompressible, non-hydrostatic Reynolds-Averaged Navier-Stokes (RANS) equations (Kim et al. 2016;Díaz-Carrasco et al. 2021). The 2-D version solves the shallow-water equations. In this study, we implemented a 2-D version of the model for its computational efficiency and our region of focus spans the Pamlico Sound and the Neuse River where severe flooding was reported. Many contemporary storm surge modeling efforts account for the impacts of nearshore waves, which have been shown to play major roles in amplifying the surge (Sheng et al. 2010;Weaver and Slinn 2005;Dietrich et al. 2011;Abdolali et al. 2021). During the landfall of Hurricane Katrina, wave set up was shown to bring more than 0.5 m of additional surge along the lower Mississippi River Delta Bunya et al. 2010). Over the study domain, however, we found the impacts of near-shore wave often muted due to the sheltering effects of the Barrier Island. As a result, we opted to exclude the modeling of near-shore wave ).

Model configuration
Several high-quality bathymetric and topographic data for this study have been combined in implementing the Delft3D FM. The entire western North Atlantic Ocean built using global SRTM15_PLUS bathymetry (Tozer et al. 2019), while NCEI Bathymetric Digital Elevation Model (30 m resolution) data has been utilized to enhance the representation of the bathymetry of the barrier island and Pamlico sound (Mulligan et al. 2019). These data are merged to create an unstructured grid mesh that spans much of the west Atlantic-the model domain extends from latitude 34.0° N and longitude 78.0° W to latitude 36.5° N and longitude 72.0° W, and mesh resolution ranges from 2000-m offshore to 30-m near the shoreline and up along major sounds and tributaries (Fig. 1).
The open water boundary, located in the deep ocean, was forced with tidal water level extracted from tidal model TPXO 9.0 (Egbert and Erofeeva 2002), whereas the upstream river boundary, located at Barnwell in the vicinity of Neuse River, was forced with observed discharge time series from USGS. However, the initial condition was fixed at mean sea level (MSL). Space-varying wind forcing from HWRF and H10 models were applied to the entire model domain. The dependency of the drag coefficient on the wind speed was specified according to Smith and Banke (1975), where a linearly varying two breakpoints, at 0 m/s (i.e., C d = 0.00063) and at 100 m/s (i.e., C d = 0.00723), were specified. The boundary smoothing time was fixed at 3600 s for the numerical parameters, and the dry cell threshold was set to 0.01 m to satisfy wetting and drying algorithm. Additionally, spatially varying bottom friction based on land use types, adopted from literature (Chow 1959;Kaiser et al. 2011), was defined using Manning's friction coefficients. In this study, the tide generating force due to earth's rotation was neglected as the model domain is not large enough to exhibit Coriolis effect. However, the horizontal eddy viscosity (i.e., 0.2 m 2 /s) and eddy diffusivity (i.e., 20 m 2 /s) were assumed to be constant for the entire domain. The default values of all other model parameters were used for this study. The simulation time of the storm event (i.e., Florence) was set from 09 September 2018 00:00:00 to 18 September 2018 00:00:00 local time using a user defined time step of 30 s with initial time step of 1 s and maximum time step of 600 s. However, due to the limitation of time window in wind forcing from HWRF, we used two days for model spin up. The model was simulated on a Linux cluster with 64 parallel processors, and the total clock time was around 8 h. Finally, the history output data was saved at 5 min interval, whereas the map outputs were saved at 3 h interval.

Atmospheric forcing and boundary conditions
The Delft3D model requires atmospheric forcings including surface wind and pressure, and tidal boundary conditions. As indicated earlier, in this study we employ two sets of wind and pressure fields to drive the hydrodynamic model, one based on the H10 parametric model and the other on the HWRF reanalysis. We interpolate the zonal and meridional wind speed (u and v, respectively), as well as surface pressure onto the Delft3D mesh using bilinear interpolation. As astronomical tide can be a crucial factor that contributes to coastal inundation during a storm event (Lai et al. 2021), in this study, we use the output from TPXO 9.0 global tidal model (Egbert and Erofeeva 2002) as outer boundary condition for our Delft3D model. We used only 8 main constituents (i.e., m 2 , s 2 , n 2 , k 2 , k 1 , o 1 , p 1 , and q 1 ) out of 37 to generate the tidal water level at the deep ocean.

Validation data sets
The wind and pressure fields from H10 and HWRF reanalysis will be validated against surface observations from five networks, namely National Data Buoy Center (NDBC) by NWS, National Water Level Observation Network (NWLON) managed by National Ocean Service, Advanced Surface Observation System (ASOS) by NWS, temporary and permanent gauging stations operated by USGS, and the network by Weatherflow Inc. The Weatherflow network consists of more than 100 stations near coastal urban areas, which are specifically designed to withstand the conditions of a landfalling hurricane with a less than 1% failure rate in surviving and recording winds up to 121 knots. The locations of these stations are shown in Fig. 2a.
To validate storm surge simulations, we collected water level series from NDBC, NWLON and USGS stations, and high-water marks (HMWs) from USGS (Fig. 2b). These data sets were remapped to NAVD 88 in order to validate the Delft3D-FM simulated water levels. It is worth pointing out that a more extensive set of HMWs are available along the NC coasts, but only those along the lower reach of the Neuse River were used for validation because the Delft3D model grid mesh is sufficiently fine in this region (~ 100 m). The list of collected data from different sources is presented in Table 1.

3
Conventional metrics were employed for judging model performance, including percentage bias (PB), Pearson's correlation (R), and root-mean-squared error (RMSE). These metrics are defined as follows: where Q obs,i and Q sim,i are observed and simulated datasets (i.e., wind speed, barometric pressure, or water level), respectively; and n is the number of records in the time series.

Validation of wind speed and pressure
The H10 and HWRF wind fields for Florence are first examined through the radial wind profiles for a 48-h window surrounding the landfall that starts from 0z 14 September. Figure 3 displays the NHC best track over the time window with the evolution of the radius of maximum wind (RMW) highlighted. It is evident that RMW increased slightly through the landfall (i.e., 37 km). Between 12 and 18z 15 September, RMW expanded greatly (i.e., 278 km), corresponding to the degeneration of the storm into a tropical depression.
Our examination will focus on the distribution of horizontal wind speed for each 6-h increment along the Southwest (225° azimuth) to Northeast (45° azimuth) transect (Fig. 3). Several authors (e.g., Hu et al. 2012) chose to analyze the distribution of wind speed along both southwest-northeast and northwest-southeast transects. We focus our attention to the former as it roughly aligns with the direction of the Carolina coastline; in particular, the wind intensity in the northeast quadrant (NEQ), which was directly towards the shore during the landfall, most likely had disproportionate impacts on the magnitude of storm surge. Figure 4 displays the wind profiles from HWRF and H10 for each 6-h increment within the 48-h window. Superimposed in each panel are surface wind speed observations at stations that fall in the northeast and southwest (SWQ) quadrants (i.e., the quadrants that intersect with the transect). Prominent observations are summarized and briefly discussed below.
Prior to landfall, the wind speeds from both models exhibit symmetry with comparable Vmax in the NEQ and SWQ ( Fig. 4a and b). The maximum wind speeds represented by H10 and HWRF are close-though HWRF appears to produce slightly higher Vmax values than those supplied by the best track at 0z (Fig. 4a). In addition, H10 wind profile declines at a faster rate with radial distance. At both 0 and 6z on 14 September (Figs. 4a and b), HWRF wind speed shows close agreement with surface observations for the NEQ, whereas the opposite is true for the SWQ (Fig. 4a and b). These suggest that the HWRF performs reasonably well in reproducing the easterly wind that blew towards the land, but  Fig. 2, circles represent the radius of maximum wind for selected synoptic hours during the landfall. The location of New Bern is highlighted in purple, and each star on the track of Florence represents 6-h increments. The southwest to northeast transect is shown as a dotted line, and the plots for this direction are shown in Fig. 4 somehow exaggerates the land-bound westerly. By contrast, H10 appears to underestimate the wind speed in NEQ and this bias is more pronounced further away from the center of the storm. Fig. 4 Validation of radial wind profiles along the azimuth angles of 45° (northeast) and 225° (southwest) produced by HWRF and H10 against surface observations (i.e., black dots) over the time window surrounding the landfall of Florence (from 0z on 14 September to 18z on 15 September). Left to right represents the southwest to northeast direction as shown in Fig. 3 as a dotted line. Note the surface observations are sampled from stations with the quadrants in which each radial direction is embedded. Here the purple circles represent V max of best track data, and green circle represents wind speed at New Bern On 06z September 14 (Fig. 4b), a secondary maximum emerges in the NEQ at about 120 km from the center. Right at landfall (12z on 14 September), HWRF wind maxima in both quadrants broadly exceed the Vmax from NHC best track. In the NEQ, HWRF wind profile is in close agreement with observations at further distance (> 120 km) from the center of the storm, whereas it is mostly above observations in the SWQ, echoing the observations in the preceding time step. H10 underestimates the wind in both quadrants compared to observations. Florence made landfall on 12z September 14 (Fig. 4c). Perhaps the most notable feature during this time is that HWRF grossly overestimated the Vmax over both quadrants (i.e., NEQ and SWQ). At this hour, the HWRF Vmax values over the two quadrants approach 80 knots, whereas the NHC estimates are about 60 knots. Again, further away from the center of the storm HWRF wind profile appears to be accurate in the NEQ, but is biased high in the SWQ. By contrast the H10 wind remains negatively biased across the transect.
After landfall, HWRF wind speed weakens rapidly over NEQ, with Vmax declining from 75 to less than 40 knots between 12 and 18z on 14 September. Meanwhile, HWRFbased Vmax for the SWQ shows a relatively minor reduction (from 80 to 65 knots), and closely follows the observations. This results in a sharp asymmetry in the HWRF wind profile. At 18z on 14 September, HWRF wind speed generally agrees with the observations along the transect. Whereas, H10 wind speed shows clear, negative bias outside the RMW that tends to be increasingly severer at further distance to the center of the storm (i.e., > 200 km).
The relative performance of H10 and HWRF for three subsequent snapshots (0, 6 and 12z on September 15) broadly resembles that at 18z on September 14, except that the HWRF profile appears to be mostly below the observed in both quadrants (Figs. 4e, f and g). H10 wind speed appears more consistent with observations within the RMW, but it declines sharply to near zero after 100 km, leading to a conspicuous underrepresentation of wind speed beyond RMW. This decline can be partially attributed to the lack of representation of background wind in the H10 model. Between 12 and 18z on September 15, RMW expanded abruptly as a result of dissipating storm intensity. On 18z, H10 wind speed is broadly higher than the observations across the entire transect, and it exhibits a curious slow rate of decline with radius. By contrast, HWRF wind speed more closely matches the observations for this hour, though it appears to be consistently lower than the latter across distance.
To diagnose factors underlying the differential accuracy of wind profiles as represented by HWRF and H10, we compare the surface pressure profiles from the two models. Figures 5a-h show the pressure profile along the NW-SE transect for the same time instants used in the wind analysis shown in Fig. 4. Broadly speaking, HWRF surface pressure matches closely with the observations throughout the 48-h window. By contrast, H10 pressure tends to be biased low further away from the pressure center, presumably because its specification of ambient pressure is much lower than the observed. In addition, H10 produces a conspicuously lower center pressure at the center relative to HWRF from 6z to 12z on 15 September. As H10 uses the center pressure from the best track directly, it is clear that HWRF was unable to fully resolve the pressure drop. Another notable observation is that the HWRF pressure profile at the last instant (18z on 15 September) is consistently higher than the observations (Fig. 5h). This feature corresponds to, and is mostly likely causatively related to the systematically lower wind speed simulated by HWRF shown in Fig. 4h. It appears that Florence's intensity declined at a somewhat slower rate on the 15th than what was predicted by HWRF. The performance statistics are presented as tabular form in the appendix.

3
To further assess the accuracy of the two sets of wind fields from the two models as the storm progressed, we compare the two wind speed time series against in situ observations at 6 stations close to the track (Fig. 3). Among these stations, B002 and B024 are offshore NDBC stations: B002 is further away from the shore whereas B024 is located near shore and close to the track. XFED and XOCR are Weatherflow stations situated along the shore; and the former is located to the site of landfall. KNKT and KNCA are ASOS stations in Cherry Point and Jacksonville, NC, respectively, and both  Table 2.
For B002, the most striking feature is that H10 produces a sizable secondary peak on 18z of 15 September 15 that is not observed by the majority of the stations (Fig. 6a). This is apparently a result of H10 not being able to accurately represent the wind structure during and after the decay of the storm into a depression. As shown in Fig. 4h, when the storm reduces to a depression, wind speed produced by H10 declines at a much slower rate with radius that in previous hours, and this results in a sharp, artificial expansion of the region with high wind speed that translates to a secondary wind peak over the periphery of the Fig. 6 Validation of wind speed time series produced by H10 and HWRF model. Shown in the first, second and third row are time series at two NDBC buoy stations offshore (a and b), two Weather Flow (WF) stations located on land along the coast; and two ASOS stations situated further inland storm. By contrast, HWRF wind track the observed time series closely, though it exhibits a persistent negative bias that is the most pronounced around the time of peak. At B024, the station near the track, this secondary peak is absent in the time series of H10 wind speed, and H10 clearly underrepresents the peak wind speed (Fig. 6b). HWRF accurately reproduced the peak, but it exhibits a negative bias during both the rising and falling limbs of the series.
Over the two Weatherflow stations (XFED and XOCR, Figs. 6c and d), the bogus secondary maximum is again evident in the H10 wind series. In addition, at both sites H10 tends to underrepresent the maximum wind speed during landfall. HWRF outperforms H10 at both sites, though there is a slight, negative bias at XOCR (Fig. 6d). For the two ASOS stations in the north (KNCA and KNKT; Figs. 6e and f), H10 conspicuously overrepresents the peak wind speed at KNCA, whereas at KNKT its accurately captures the peak. Once again, at both sites the bogus secondary peak is present and the outperformance of HWRF is evident. Figure 7 shows the correlation coefficients of H10 and HWRF wind fields against surface observations for each individual station. Note that the HWRF reanalysis was on a 6-h resolution and was first interpolated into hourly intervals prior to computing the correlation. For H10, the correlation contrasts sharply between stations situated to the north and south of the tracks. For the stations in the north, correlation is broadly poorer (< 0.4) with only over a few stations close to the track exhibiting values higher than 0.5. By comparison, a cluster of stations to the south of the track exhibit good correlation-almost all of these stations are located in the vicinity of the track. Further away to the south, correlation declines sharply. For HWRF, the correlation is generally good for a majority of stations, though it tends to be relatively low for a cluster of stations further north along the Chesapeake Bay.
The bias in the peak wind speed as represented by the two wind data sets at each station is shown in Fig. 8. On both sides of the track, the peak wind speed from H10 is positively biased for a majority of stations further away from the track; whereas close to the track, it is negatively biased or neutral for a slight majority of stations. The negative bias near the track is attributable to the fast decline in the H10 wind speed away from the center, which makes it unable to reproduce the peak wind speed at the stations with moderate distance to the track. By contrast, the positive bias for the far-away stations is related to the artificial increase in wind speeds as featured in H10 after the storm reduced to a depression. To elaborate, the observed peak wind speeds during Florence's landfall were weak at these stations; so were the coincidental wind speeds from H10. For H10, the peak wind speeds  actually arrived later in the series during the depression phase of the storms, and these values were broadly higher than the observed, leading to the positive bias. The bias in the HWRF peak wind also varies widely among stations. Relative to H10, HWRF product is nearly bias-neutral for a significant number of stations to the north of the track. For the remaining stations the bias is mixed: there are a few stations along the NC coast where the bias is clearly positive. To the south of the track, the bias becomes overall negative further away from the track.

Comparison of simulated water level and inundation extent
The Delft3D storm surge simulations using wind and pressure fields from H10 and HWRF reanalysis are first validated against in situ observations. Two validation stations located to the north of the track are selected for this purpose (Fig. 9). The first one is a COOPS station (ID 8658163) situated offshore of Wrightsville Beach, NC, and the second is a temporary USGS gauge placed in the Neuse River near New Bern (Table 4). Note that there are several other COOPS stations with water level records during the event, but these are not used as their locations are either too far away from the storm center, or over areas shielded from storm surge.  Tables 3 and 4. For reference, the wind time series are shown alongside the water level series in each plot. At the Wrightsville Beach station, as indicated earlier the HWRF wind series features two sharp drop-offs that are not consistent with the observations (Fig. 10a). Florence Fig. 10 Wind and water level times series at COOPS station 8658163 (Wrightsville Beach, NC) that is located close to the track: a surface wind speed from H10, HWRF and surface station, and b water levels produced by Delft-3D simulations driven by H10 and HWRF along with the observations 1 3 Fig. 11 As in Fig. 10, except at the USGS station near New Bern (USGS 02092576) that is located in the Neuse River at around 95 km from the shoreline produced a small surge around 18UTC on the 14th, with maximum water level reaching 2 m (Fig. 10b). Delft3D simulation driven by HWRF yields a surge on the same day but 12 h earlier, and the maximum water level in the next tidal cycle declined to normal level (~ 1 m). By contrast, using H10 as the forcing results in slightly more accurate water level series on the 14th than that produced by HWRF: the simulated surge level near 6 UTC remains positively biased but is slightly lower than that by HWRF, and the peak surge level near 18 UTC is much closer to the observed. Another noticeable difference is that using H10 leads to depressed ebb levels throughout the time window, whereas HWRF-driven simulation largely resolves the tidal cycles. The inability of HWRF-driven simulation to reproduce the peak water level, as we surmise, is a consequence of the fast decline in HWRF wind on the 14th which was associated with the passage of the eye (Fig. 9). At the USGS station near New Bern, the contrast between H10 and HWRF wind series is stark, broadly echoing that shown earlier in Fig. 6. H10 produces an earlier rise in wind speed, overpredicts the peak wind, and features a sharp drop-off to 15 September. As has been shown earlier in Fig. 6. H10 wind series features a bogus secondary peak on the 15th that is related to the rapid expansion in the RMV during the weakening of Florence to a depression. On the other hand, HWRF wind series closely track the observations, though it features a secondary peak about 12 h following the primary one. Simulated water level series driven by H10 lags slightly behind that of observed, and, somewhat paradoxically, the peak level is visibly lower than the observed despite the higher wind speeds during the landfall over the location as featured by H10. By comparison, those forced by HWRF are closely correlated with the observations, and the peak surge is nearly perfectly captured by the model. In addition, in the H10-driven simulations, the bogus secondary wind peak on 15 September translates to a distinctive bogus spike in water level.
The summary validation statistics including bias, correlation and RMSE all point to broad outperformance of HWRF wind fields over both locations, but the impacts on water level varies (Tables 3 and 4). For the Wrightsville Beach station, H10 performs better in terms of correlation, but worse as indicated by RMSE (Table 3). This difference is reflective of the lower ebb levels seen in Fig. 10. For the New Bern station, HWRF-related wind and water series both outperform by a wide margin (Table 4). Figure 12 compares the maps of maximum inundation extents computed from the H10 and HWRF-driven simulations, and Fig. 12(c) shows the difference field. Broadly speaking, H10 produces higher surge inland and upstream of the Neuse and Pamlico River, and along the NC coast extending from the Wilmington (near the landfall) to Morehead City. HWRF, by contrast, produces higher surge over an area stretching from the coast of Pamlico Sound to the lower reaches of the two major rivers. This contrast is clearly a product of the differing radial wind profiles as represented by the two models demonstrated earlier (see Fig. 4). Prior to and during landfall, H10 underestimates wind speed away from storm center (> 100 km), and this reduces the momentum for propelling the storm surge over the Northern portion of the domain. On the other hand, after landfall, H10 produces artificially high wind speed at closer range, apparently an indication that its current structural framework is unable to realistically account for the reduction in wind speeds as a result of increased friction on land. H10's artificially high peak wind speeds close to the eye naturally translate into higher storm surge near the track and upstream of the estuaries.

Fig. 12
Maps showing maximum flood extent, generated from 3-hourly output water level of Sep. 13 and 14, driven by wind and pressure fields produced by a H10, and b HWRF. c shows difference of simulated maximum water level (m) by Delft-3D forced by wind and pressure fields of HWRF and H10. The maximum water level is computed for each cell from 3-hourly series over 13-14 September (surrounding the landfall). Positive/negative values indicate higher maximum surge produced using HWRF/H10 product 1 3 Fig. 13 Scatter plots of simulated maximum inundation depths versus HWM observations along the Neuse River; a Maximum surge with H10, and b Maximum surge with HWRF Figure 13 shows the validation of simulated maximum surge produced by Delft3D driven by H10 and HWRF-driven against HWMs collected downstream of Neuse River (see Fig. 2 for locations). Both simulations are closely correlated with observations, but that driven by HWRF features lower RMSE and higher correlation. In terms of bias, the HWRF-driven maximum surges show a slightly positive overall bias, whereas those based on H10 are negatively biased. These results are consistent with the observations made in comparing the time series at New Bern (Fig. 11), and suggest that the HWRF reanalysis is likely a superior source forcing for the surge simulations at least along the Neuse River and over the adjacent offshore locations.

Discussions
Accurate meteorological forcing is a key ingredient in the analysis and prediction of coastal flooding caused by storm surges. For decades, parametric wind fields have been a major source of forcing input for storm surge prediction and analysis. Thus far, very few studies have touched upon the relative strengths of wind fields as derived from parametric models versus those based on NWP model, or the impacts of differential accuracy of wind fields on the fidelity of storm surge simulations. The present study addresses this gap by offering a detailed assessment of Hurricane Florence wind fields as represented by a parametric model (H10) and the HWRF reanalysis, and the storm surge simulations driven by the respective data set.
The comparisons underscore several key shortcomings of H10 model. These include its overly sharp decline with distance, its tendency to overpredict wind speeds on land, and its overall inability to resolve wind fields after the storm weakened into a depression. These shortcomings are apparent reflections of structural limitations of H10 model and the premises on which the model was formulated. Note that H10 improves upon the original H80 in major respects, such as relating pressure drops directly to wind at the surface, rather than at the gradient level, relaxation of the cyclostrophic balance assumption for the inner core, and the use of new regression relations and parameter b s . Nonetheless, it is evident that these improvements are inadequate for the model to reproduce the wind fields over the periphery of the storm, or to capture the complex, rapid evolution of vortex structure after landfall. Earlier studies, notably Hu et al. (2012), found positive bias in H80 wind fields and attempted to remedy this bias by incorporating canopy-based adjustment factors. Apparently, this bias remains an issue in H10 wind fields at least during the landfall of Florence despite the aforementioned enhancements. Note that the bias is not entirely a result of structural inadequacy of H10: our analysis reveals suspiciously high Vmax in the best track data on which the model relies on to derive radial profile, and this can be a major contributor to the positive bias. In addition, H10's lack of representation of realistic ambient pressure, and more precisely its inability to reproduce the sharp transition in temperature gradient, may have played an important role in rendering the negative bias in the farther range. Further, gradient wind or cyclostrophic assumptions may become increasingly poor approximations of the wind-pressure relationship as the storm weakened after landfall. Kepert (2001), for example, postulated the existence of a jet-like feature in the boundary layer of tropical storms, which later was verified by empirical observations (Hirth et al. 2012). Moreover, in situ observations hint the presence of a secondary wind maximum at farther range of Florence during and after landfall. Unfortunately, though H10 offers a mechanism for explicitly representing the secondary maximum, the implementation of this scheme is deemed impractical as it requires observational data to identify and define the secondary maximum-such data are hardly available a priori. How to leverage remotely sensed product such as brightness temperature from satellite imagers or sounders for this purpose will be a topic of future research.
While our analyses confirm the realism of HWRF wind and pressure reanalysis, these also uncover a few shortcomings of the product. Perhaps the most glaring is the sharp progression of the bias through landfall. Prior to, and even during the landfall, the radial profile of HWRF reanalysis exhibits a slightly positive bias over both NEQ and SWQ. After landfall, however, the bias became progressively negative. It should be noted that the bias calculations were not exactly rigorous, in that the in situ stations used for the comparisons over each quadrant are not situated precisely along the azimuth angle for which the HWRF wind profile was drawn. This caveat aside, this transition in bias after landfall is conspicuous enough to warrant close attention-it may well be reflective of potential mechanistic deficiencies in HWRF that inhibit its ability to accurately resolve the dissipation phase of the Hurricane. Possible mechanisms include inadequate representations of land surface conditions, over-smoothing of the wind profile of the 2-D vortex in creating the bogus vortex that results in discrepancies between model simulated and observed maximum wind speed, and the lack of assimilation of surface observations after landfall. In particular, high soil moisture and inundation as a result of heavy rain are known to help sustain the intensity of tropical storms after landfall through the so-called brown ocean effects (Nair et al. 2019;Yoo et al. 2020). Florence produced torrential rain during its landfall, and this may have helped slow down the dissipation of the storm. The degree to which HWRF's coupling scheme resolves the interplays between land surface and atmosphere is a topic that requires additional scrutiny.
The mixed results from the comparisons of Delft3D simulated storm water levels driven by the two sets of wind products are in fact illuminating. Using the H10 wind and pressure fields as input leads to higher storm surge in regions near the track, consistent with the observation that H10 tends to feature higher wind intensity close to the center of the storm. On the other hand, H10-driven model simulation underpredicts the surge at the New Bern even though the H10 features higher peak wind speed locally. This seeming contradiction points to the fact that the magnitude of storm surge is not determined exclusively, or even strongly, by local wind speed and direction, but by a complex aggregate of wind/pressure offshore as well as over land, geometries of coastline and estuaries, and the interplays between these factors. Comparisons of maximum surge level between H10 and HWRF clearly illustrate that the latter is able to induce higher surge over much of the lower Pamlico Sound stretching from the Barrier Island to New Bern, apparently a reflection of the ability of HWRF to resolve wind fields over further distance prior to the landfall. These findings collectively underscore the challenges in quantifying errors in storm surge simulations that arise from errors in simulated wind fields.

Conclusions
The wind validation was performed using data from 95 observing stations from both public and private sources. These include National Ocean Service (NOS) stations deployed nearshore, United States Geological Survey (USGS) temporary sensors along major rivers, Advanced Surface Observation System (ASOS) stations in airports, and sensors operated by Weather Flow Inc along the coast. In this study, a hydrodynamic model Delft3D was configured to simulate storm surge along the southeast using the Holland (2010) known as H10 and the Hurricane Weather Research and Forecasting (HWRF) wind/pressure fields. In order to minimize the complicating effects of model calibration, the model incorporates simple, spatially uniform roughness coefficients which underwent only light calibration. Key findings are summarized as follows: 1. The HWRF model captures the wind and pressure fields more accurately compared to the H10 model before the landfall of Hurricane Florence (2018). The latter features lower wind speeds away from the storm center, possibly an outcome of lacking representation of ambient wind. 2. During the landfall, H10 performs slightly better for the stations located within 100 km from the storm center, whereas HWRF tends to overpredict the peak wind. Yet, H10 wind speed drops sharply further away from the storm center, resulting in large negative biases at those validation stations. 3. After landfall, Florence's strength declines rapidly. H10 is unable to reproduce the decline over the inner range (close to the storm center), resulting in large positive biases across stations. Over the outer range, however, H10's rapid drop-off in wind speed leads to broad negative biases. Note that the positive bias near storm center is partly a result of overly high maximum wind speed (Vmax) values in the best track data. 4. HWRF more accurately depicts the evolution of Florence wind fields in time after landfall, yet it produces overly suppressed wind speeds across the southwest-northeast transect and this suppression is particularly pronounced on the 15th. 5. After the weakening of Florence into a depression, the Radius of Maximum Wind Speed (RMV) expands drastically, and H10 is unable to produce realistic wind fields. 6. Storm surge simulations driven by HWRF and H10 yield mixed outcomes. HWRFdriven simulation accurately reproduces the surge near New Bern, NC, whereas H10driven simulation features a slightly lower, and delayed peak. In an offshore station near the storm track, using H10 as forcing leads to a slightly better depiction of the magnitude and timing of the surge. 7. Inundation depth produced by HWRF-driven simulation is conspicuously higher than that from the H10-forced simulation over the downstream portions of Neuse and Pamlico Rivers. By contrast, it is broadly lower offshore, along the upper reaches of the two rivers, and over areas close to the track.
In broad terms, the study illustrates the strength of HWRF model in reproducing the radial wind profile of Hurricane Florence, in depicting the evolution of the wind fields during and after the landfall, and in capturing the spatial patterns of wind through during the weakening phase of the storm on the 15th. When compared against HWRF, the shortcomings of H10 model are evident. These include its overly sharp decline with distance, its tendency to overpredict wind speeds on land, and its overall inability to resolve wind fields after the storm weakened into a depression. These shortcomings are apparent reflections of structural limitations of H10 model and the premises on which the model was formulated. It is clear that improvements in structure and parameter estimation scheme introduced in H10, including the relaxation of the cyclostrophic balance assumption for the inner core, and the use of new regression relations and parameter b s , are inadequate for the model to reproduce the wind profile on land. Hu et al. (2012) noted a similar positive bias of H80 wind field on land, and attempted to remedy this bias by incorporating canopy-based adjustment factors. It is, however, worth drawing the distinctions between H10 and H80, as the latter explicitly relies on the cyclostrophic balance and has a rather imprecise definition of elevation associated with its wind profile whereas the former does neither. One major contributor to H10's positive bias in the inner range is the suspiciously high Vmax in the best track data on which the model relies on to derive radial profile. While in situ data used in this study was insufficiently dense near the RMV to directly appraise the validity of the Vmax, they point to a distinct possibility that the Vmax is biased high in the later part of the window. In addition, H10's lack of representation of realistic ambient pressure may have played an important role-it leads to artificially suppressed pressure gradient which translates into lower wind speed. Further, the structure of H10 wind profile is perhaps broadly unsuitable for modeling the wind field of tropical storms after landfall as these are strongly modulated by interactions of the storm with terrain features. It is worth pointing out in situ observations hint the presence of a secondary wind maximum at farther range of Florence during and after landfall. H10 does offer a mechanism for representing this maximum, but it was not implemented due to a lack of a priori information to establish the magnitude and location of this maximum, and the concern of barotropic instability.
Similar to many studies conducted earlier on retrospective analysis of storm surge, our investigation was constrained by data availability and computational demand. Owing to a shortage of in situ observation and the focus on a single storm, we were unable to perform detailed, spatially distributed validation of storm surge simulations along the Pamlico sound and adjacent land, or assess the ability of model calibration to compensate for errors in forcings. As a result, a number of questions concerning the fidelity of the two sets of surge simulations remain unanswered. In addition, for computational tractability, some of the mechanisms that impact the surge were omitted. In this study, the impacts of wave setup were assumed negligible owing to the consideration of the unique geographic of the study region, where such impacts were likely subdued due to the presence the Barrier Islands. Extending the comparisons to cover additional landfalling tropical storms, and surge cases over different geographic domains, will offer additional insights into the differential strengths of H10 and HWRF, their structural underpinnings, and manifestations of mechanistic deficiencies in wind fields in surge simulations and predictions. In particular, it is of great interest to assess the effects of errors in the parametric wind models associated with the lack of explicit representation of the ambident wind fields on surge in situations where wave setup likely plays prominent a role through coupled wave-surge simulations, and to investigate potential mechanisms for alleviating such effects, e.g., inclusion of additional velocity-range pairs and superposition of parametric wind fields on ambient winds from numerical weather model simulations.