Discrepancies on Storm Surge Predictions by Parametric Wind Model and Numerical Weather Prediction Model in a Semi-Enclosed Bay: Case Study of Typhoon Haiyan

: This study explores the discrepancies of storm surge predictions driven by the parametric wind model and the numerical weather prediction model. Serving as a leading-order storm wind predictive tool, the parametric Holland wind model provides the frictional-free, steady-state, and geostrophic-balancing solutions. On the other hand, WRF-ARW (Weather Research and Forecasting-Advanced Research WRF) provides the results solving the 3D time-integrated, compressible, and non-hydrostatic Euler equations, but time-consuming. To shed light on their discrepancies for storm surge predictions, the storm surges of 2013 Typhoon Haiyan in the Leyte Gulf and the San Pedro Bay are selected. The Holland wind model predicts strong southeastern winds in the San Pedro Bay after Haiyan makes landfall at the Leyte Island than WRF-ARW 3 km and WRF-ARW 1 km. The storm surge simulation driven by the Holland wind model ﬁnds that the water piles up in the San Pedro Bay and its maximum computed storm surges are almost twice than those driven by WRF-ARW. This study also ﬁnds that the storm surge prediction in the San Pedro Bay is sensitive to winds, which can be a ﬀ ected by the landfall location, the storm intensity, and the storm forward speed. The numerical experiment points out that the maximum storm surges can be ampliﬁed by more 5–6% inside the San Pedro Bay if Haiyan’s forward speed is increased by 10%. This study has explored the discrepancies of the storm surge predictions driven by a parametric wind (PW) model and a numerical weather prediction (NWP) model. By shedding light upon the discrepancies, the storm surges of 2013 Typhoon Haiyan in the Leyte Gulf and the San Pedro Bay have been studied. The frictional-free, steady-state, and geostrophic-balancing solutions on the storm winds and sea-level pressures given by the Holland wind model (HWM) have been first verified against the observed data at the Guiuan weather station. The 10-m winds predicted by WRF-ARW (Weather Research and Forecasting-Advanced Research WRF), which solves the 3D time-integrated, compressible, and non-hydrostatic Euler equations, have been compared with the ones by HWM. In the results of computed storm surges driven by HWM, the storm surges of more than 4.5 m have been found in the San Pedro Bay due to the water piling up and the pattern is agreed well with the field survey data. However, in the results by both WRF-ARW 3 km and WRF-ARW 1 km, the maximum storm surge in the San Pedro Bay is almost half than that driven by HWM. The significant differences among HWM, WRF-ARW 3 km, and WRF-ARW 1 km are caused by the predicted landfall location of Haiyan, the storm intensity, and the storm forward speed. The numerical experiment focused on the storm forward speed has pointed out that the extreme values of storm surges inside the San Pedro Bay can be intensified by 4.5–7.0% when the storm forward speed is increased by 10%. This study hopes to pave the way for future research in understanding the discrepancies in predicting storm surges driven by either of the PW or NWP models and thus benefiting storm surge predictions in a semi-enclosed bay.


Introduction
Storm surges caused by tropical storms can have a significant impact on coastal facilities [1]. In some cases, storm surges can interact with tides [2,3] and wind waves [4][5][6]. When a storm approaches coastal regions, winds can amplify the height of storm surges locally at least two times more than that in the open ocean [5,6]. Hence, to perform an accurate storm surge modeling, having a better understanding of storm-wind fields is critical.
In the Western Pacific Ocean, more intense typhoons potentially occur in the future and the low-lying regions become more vulnerable to storm surges [7,8]. To minimize the disaster of storm surges, the accuracy of predicting storm surges needs to be improved. The uncertainties of storm surge predictions can attribute to the input data (e.g., bathymetry and winds), physics consideration (coupling with waves, tides, river flows, etc.), and numerical settings (e.g., grid size, time step, and an

Storm Surge Model
In this study, the depth-integrated storm surge model solving the set of linear shallow water equations with the Coriolis effect and the bottom friction is adopted. This storm surge model is based on the COMCOT (COrnell Multi-grid COupled Tsunami) tsunami model. COMCOT has been validated by the solitary wave run-up benchmark problem [44] and been used to study tsunami events such as 1960 Chilean Tsunami [45], 2004 Indian Ocean Tsunami [46,47] and 2011 Japan Tohoku Tsunami [48]. Furthermore, COMCOT has been developed into a cloud computing platform [48]. The external forcing terms of the sea-level pressures and the wind shear stresses are added to the momentum equations of COMCOT to simulate the storm surge propagation.
The mass and momentum equations of the linear shallow water equations in the Cartesian coordinate are presented as follows: ∂P ∂Q ∂t Water 2020, 12, 3326 4 of 27 where η is the water surface elevation, t is the time, x,y are spatial notations, (P,Q) are the depth-integrated volume fluxes per unit length, H is the total water depth (H = h + η), h is the still water depth, g is the gravitational acceleration (=9.81 m/s 2 ), f is the Coriolis parameter ( f = 2ωsinϕ), ω is the Earth angular velocity (=7.2921 × 10 −5 rad/s), (τ sx , τ sy ) are the wind shear stresses, (τ bx , τ by ) are the bottom frictional shear stresses, ρ w is the water density (=1000 kg/m 3 ), and P a is the sea-level air pressure. The Manning's formula, originally from the conception of the open-channel flow, is used to model the bottom friction: τ by = ρ w gn 2 where n is the Manning's roughness coefficient. By introducing appropriate Manning's coefficient, the roughness of the material can be determined in this formula. As shown in Equations (4) and (5), the bottom friction is proportional to the squared of volume fluxes per unit length and the Manning's coefficient and inversely proportional to the water depth to the power of 7/3. Hence, the bottom friction using this formula is less important in deep waters, more significant in shallow waters, and crucial in coastal regions (including inland areas). Since the inundation calculation is not involved in this study, the Manning's coefficient is 0.025 over the whole computational domain, which is also adopted by other studies [18,19]. In fact, the Manning's coefficient shall be determined carefully, especially the inundation calculation involved. The use of the Manning's roughness coefficient to storm surge modeling has been investigated [49][50][51].
The quadric law is used to model the wind shear stresses on the water surface: τ sy = ρ a C d u 10 u 10y , (7) in which C d is the wind drag coefficient between the water surface and the air, u 10 is the 10-m wind speed, (u 10x , u 10y ) are the components of the 10-m wind speed in the xand y-directions. The wind drag coefficient of Wu (1982) is used [52]: C d = 1.2875 × 10 −3 , u 10 < 7.5, (0.8 + 0.065) × 10 −3 , u 10 ≥ 7.5.
The boundary conditions along the computational domain are given by the free surface elevations caused by the pressure drop (i.e., the difference between the sea-level pressure and the ambient pressure) in wet cells: This storm surge model is discretized by the explicit finite difference method in time and space. The leap-frog scheme is used to discretize the free surface elevation at (n + 1/2)∆t and the volume fluxes per unit length at (n + 1)∆t [53,54]. The forward finite difference method is used to solve the mass equation and the momentum equations. The staggered grid system, Arakawa C grid, is used [55]. Hence, the volume fluxes are defined at the boundaries of cells and the free surface elevations are in the center of cells. The boundary conditions between wet and dry cells are non-flux boundary conditions. It implies the volume fluxes of water cannot go across the coastline (i.e., the inundation calculation is not considered and the coastline is fixed in the simulation). Hence, the mass conservation is solved only in wet cells.
The tide effect and its interaction with storm surges are not considered in this study because the predicted tidal level was 0.15 m (about 2-3% of the peak surge) when the peak surge occurred at Water 2020, 12, 3326 5 of 27 00:00 UTC 8 November 2013 [13,18] and the maximum tidal range was within 1.0 m in the San Pedro Bay [35]. The wave radiation stresses are not considered since they are insignificant in generating storm surges along the coasts of the Leyte Gulf [13,56] and the maximum wave-enhanced storm surge is about 0.3 m (approximately 6% of the peak storm surges) at Tacloban.

Parametric Holland Wind Model
In this study, the well-known parametric Holland wind model (hereafter, HWM), is adopted to reconstruct the meteorological fields of Typhoon Haiyan. HWM has been widely used in storm surge modeling and wind-wave modeling [30,57]. The exponential distribution of the atmospheric pressure proposed by Schloemer (1954) [58] is extended by Holland (1980) [24]: where P c is the central storm pressure, ∆P is the pressure drop (∆P = P n − P c ), P n is the ambient pressure, R max is the radius of maximum wind, r is the distance from the specified grid to the storm center, and B is the scaling parameter. Note that P n can vary in the regions from about 1010 to 1013.25 hPa [59]. For the sake of simplicity, P n = 1013.25 hPa is used in this study. By using the gradient wind balance equation, the solution of the radial wind profile with the Coriolis effect is yielded [24,60]: V g is the computed gradient Holland wind speed. Note that this solution is steady axisymmetric and geostrophic. The maximum gradient wind speed exists where r = R max where the Coriolis force is relatively insignificant compared to the pressure gradient and the centrifugal force, As shown in Equation (12), the maximum gradient wind speed is proportional to the square root of the scaling parameter B, and the pressure drop ∆P.
The 1999 Harper and Holland's scaling parameter B for the computed pressure and wind profiles is shown below [60]: where the scaling parameter B is limited from 1.0 to 2.5. Note that if the scaling parameter B = 1, the pressure distribution proposed by Schloemer (1954) [58] is assumed. The radial distribution of the wind speed from the storm center is already given by Equations (10)- (13). Now, the wind directions are needed to define. As the assumption of a steady axisymmetric storm, the wind direction is cyclonic toward the storm center. However, while the storm is moving forward, the wind flows are crossing the isobars [60,61], which indicates an inward angle shall exist. In this study, the inward angle of 25 degrees is used [62]. By adopting the assumption of the cyclonic winds toward a storm center and the inward angle deviating from the isobars, the vector form of the gradient winds V g is yielded from the function form of Equation (11).
The correction factor is introduced here to convert the computed gradient Holland wind speed to the wind speed on the standard level, the 10-m height above the mean sea-level: Water 2020, 12, 3326 where K m is based on the logarithmic law, in which z 0 is the roughness length and z g is the height of the gradient wind. However, z 0 and z g are hard to define during intense tropical storms. Alternatively, K m can be approximately used as 0.7 for the 10-m wind conversion of the Holland wind model [60]. Further discussions of K m on storm surge and wind wave modeling can be found in the studies such as Chueng et al. (2003) [30] and   [57]. While a storm is moving, the wind pattern becomes asymmetric because winds are amplified on the right and weaken on the left in a storm forward direction. The correction term for an asymmetric wind pattern is proposed [23]: U cor is the correction term and V f is the forward velocity of a storm. The maximum corrected term exists where r = R max and is equivalent to the half of the forward storm velocity.
Finally, by assuming the exponential pressure distribution with the scaling parameter, using the gradient wind balance equation, and considering the effect of the storm forward motion, the 10-m winds are expressed by Note again that the 10-m winds come from the parametric wind model used in this study are frictional-free, steady-state, and geostrophic-balancing solutions.

WRF-ARW Model
In this study, the 10-m wind fields and sea-level pressures performed by the WRF-ARW (Weather Research and Forecasting-Advanced Research WRF) model are involved. WRF-ARW solves the 3D time-integrated, compressible, and non-hydrostatic Euler equations [10]. The model results of WRF-ARW version 3.8.1 used in this study are Cases 3 km F2 and 1 km F2 of Kueh et al. (2019) [15]. They will be denoted as WRF-ARW 3 km and WRF-ARW 1 km in the following discussion. The horizontal domain of WRF-ARW is from (117.88 • E-147.12 • E) and (0.0 • N-19.62 • N), and 45 vertical levels are employed. More modeling details and numerical settings of WRF-ARW can be found in Kueh et al. (2019) [15].

Storm Track and Parameters of Haiyan
The storm track of Haiyan is illustrated in Figure 1a, and the best-track parameters from the JMA (Japan Meteorological Agency) database [63] are listed in Table 1. The JMA best-track database provides the best-track parameters of a storm every 6 h, such as the central pressures, the maximum sustained wind speeds, and the locations of a storm center. As shown in Figure 1a, Typhoon Haiyan passed over the central regions of the Philippines and made landfall at the Samar and Leyte Islands between 18:00 UTC 7 November 2013 and 00:00 UTC 8 November 2013. In Figure 1b, Haiyan passed over the Leyte Gulf and the San Pedro Bay. From the best-track parameters listed in Table 1, Haiyan has the lowest central pressure of 895 hPa and sustained wind speed of about 64.30 m/s at 12:00 UTC and 18:00 UTC 7 November 2013. It implies that Haiyan made landfall in the Philippines while it has a peak storm intensity.   The parameters required in HWM have been described in Table 1, except for the radius of the maximum wind (RMW, R max ). RMW is an important parameter which decides the scale of the storm winds. The effect of R max to storm surge modeling has been investigated by Kim et al. (2015) [18] and R max of 50 km with the leveling-off drag coefficient fulfills the survey data. However, smaller R max (accurately, R max = 23 and 30 km) is used to perform the storm surge simulation of Haiyan by Kumagai et al. (2016) [19]. However, the pressure distributions used in Kim et al. (2015) and Kumagai et al. (2016) are different from the one in this study. Therefore, RMW shall be examined again. In the following section, the effect of RMW will be discussed and verified against the observed weather data at the Guiuan weather station.

Computational Setup of Storm Surge Modeling
The computational domain of the storm surge model is shown in Figure 2a. The domain is designed to study the storm surges induced by Haiyan in the Leyte Gulf and the San Pedro Bay where had been impacted dramatically. The model domain is from (124.9 • E-126.1 • E) and (10.0 • N-11.4 • N). As shown in Figure 2a, the coastline along the Samar Island to the Leyte Island is similar to the funnel shape, which makes this region more susceptible to storm surges [35]. In the San Pedro Bay, the water depths are less than 30 m, which depicts that the wind shear stresses can play an important role in intensifying storm surges locally. The source of the bathymetry is from GEBCO (the GEneral Bathymetric Chart of the Oceans) 2014 grid [64]. The grid size of the computational domain is 277.7 m over the whole computational domain. The model time step is 0.5 s to fulfill the Courant-Friedrichs-Lewy (CFL) condition in the explicit finite difference model [65]. The storm surge simulations are performed from 00:00 UTC 7 November 2013 to 12:00 UTC 8 November 2013, a total of 129,600 s. In our study, because of insufficient information of bathymetry and land data for the San Juanico Strait, our simulations end in the frontal part of the strait (see Figure 2a) and the boundary condition used is non-flux. The storm surges induced by Haiyan around the San Juanico Strait have been discussed by Bricker et al. (2014) [66]. They indicated that Haiyan's storm surges can extend up from the San Pedro Bay through the San Juanico Strait to Bogulibas and rapidly dissipate where the strait narrows. Note that the results reached by Bricker et al. (2014) [66] were based on the simulations using rough bathymetry from GEBCO and topography from SRTM (Shuttle Radar Topography Mission) despite some local detailed data such as the nautical charts were adopted. The results are expected to be better if more accurate bathymetry and land data can be acquired in future works.  Table 2.

Validation of Parametric Holland Wind Model
Before performing the storm surge simulation, the results of HWM shall be verified. Hence, in this section, the 10-m winds and the sea-level pressures of HWM are compared to the observed data at the Guiuan weather station, East Samar. In the Guiuan weather station, the wind speeds are measured on the vertical level of 60 m and are converted to the standard level of 10 m using the power-law formula [67]. The observation data of the Guiuan weather station are digitalized from Kim et al. (2015) [18]. This study chooses the RMW of 30 km following Kumagai et al. (2016) [19] and 50 km from Kim et al. (2015) [18]. The study of Kumagai et al. (2016) [18] used Rmax = 23 and 30 km to reconstruct the storm winds of Haiyan. This study fixes smaller RMW in the storm wind modeling to have a fair comparison and analyze the differences in the use of RMW.
In Figure 3, the results of HWM using the pressure distribution scaled by the parameters of Harper and Holland (1999) [60] are compared to the weather data measured at the Guiuan weather station. In Figure 3a, the observed sea-level pressures start to decrease gradually from about 04:00 UTC 7 November 2013, drop dramatically from 12:00 UTC 7 November 2013, and then reach the  Table 2. The time histories of storm surges are provided at the specified numerical gauges in the storm surge model. Figure 2b shows the map of these numerical gauges, and Table 2 lists their locations. Note again that the boundary conditions between wet and dry cells are non-flux (i.e., the volume fluxes of water are not allowed to cross the coastline and set to zero along the coastline). Hence, the locations of the numerical gauges are moved two to three grids away from the coastline to avoid the land effect caused by the boundaries.

Validation of Parametric Holland Wind Model
Before performing the storm surge simulation, the results of HWM shall be verified. Hence, in this section, the 10-m winds and the sea-level pressures of HWM are compared to the observed data at the Guiuan weather station, East Samar. In the Guiuan weather station, the wind speeds are measured on the vertical level of 60 m and are converted to the standard level of 10 m using the power-law formula [67]. The observation data of the Guiuan weather station are digitalized from Kim et al. (2015) [18]. This study chooses the RMW of 30 km following Kumagai et al. (2016) [19] and 50 km from Kim et al. (2015) [18]. The study of Kumagai et al. (2016) [18] used R max = 23 and 30 km to reconstruct the storm winds of Haiyan. This study fixes smaller RMW in the storm wind modeling to have a fair comparison and analyze the differences in the use of RMW.
In Figure 3, the results of HWM using the pressure distribution scaled by the parameters of Harper and Holland (1999) [60] are compared to the weather data measured at the Guiuan weather station. In Figure 3a, the observed sea-level pressures start to decrease gradually from about 04:00 UTC 7 November 2013, drop dramatically from 12:00 UTC 7 November 2013, and then reach the minimum values of 910 hPa at 22:00 UTC 7 November 2013. Unfortunately, the measurement was interrupted because the weather station was damaged by Haiyan [68]. Using 1999 Harper and Holland's pressure distribution with R max = 50 km, the peak value of predicted sea-level pressures is close to the observation. For the results using R max = 30 km, they decrease at 15:00 UTC 7 November 2013, which changes are slower than the ones using R max = 50 km. Both of them slightly over-predicted the sea-level pressures from 06:00 UTC to 15:00 UTC on 7 November 2013 and give earlier arrival time of the peak sea-level pressure, 1 h before the observation. In Figure 3b, using R max = 50 km over-predicts wind speeds. Compared to the use of R max = 50 km, the predicted winds using R max = 30 km match closer to the converted 10-m winds.
interrupted because the weather station was damaged by Haiyan [68]. Using 1999 Harper and Holland's pressure distribution with Rmax = 50 km, the peak value of predicted sea-level pressures is close to the observation. For the results using Rmax = 30 km, they decrease at 15:00 UTC 7 November 2013, which changes are slower than the ones using Rmax = 50 km. Both of them slightly over-predicted the sea-level pressures from 06:00 UTC to 15:00 UTC on 7 November 2013 and give earlier arrival time of the peak sea-level pressure, 1 h before the observation. In Figure 3b, using Rmax = 50 km overpredicts wind speeds. Compared to the use of Rmax = 50 km, the predicted winds using Rmax = 30 km match closer to the converted 10-m winds.
Without using the 1999 Harper and Holland's scaling parameter, the pressure distribution of Schloemer (1954) [58] is assumed. As shown in Figure 4a, the results using Rmax = 30 km match the trend of the observed sea-level pressures closer than the ones using Rmax = 50 km. However, both of them underestimate the peak value of the sea-level pressures measured at the Guiuan weather station. The comparison against the observed winds is illustrated in Figure 4b. In Figure 4b, both results do not match the observed wind pattern, which increases quickly from 14:00 UTC 7 November 2013. The predicted winds increase gently from 00:00 UTC 7 November 2013 and reach the peak values at about 21:00 UTC 8 November 2013.
Despite using the 1954 Schloemer's scaling parameter with Rmax = 30 km seems to have a better comparison with the observation data in terms of the sea-level pressures, the predicted winds cannot match the dramatic increase from 15:00 UTC 7 November 2013. By verifying against the observation data at the Guiuan weather station, the use of the pressure distribution scaled by the 1999 Harper and Holland's parameters with Rmax = 30 km has the best agreement in the wind comparison. Hence, the results of HWM presented in the following sections will use this configuration.  Without using the 1999 Harper and Holland's scaling parameter, the pressure distribution of Schloemer (1954) [58] is assumed. As shown in Figure 4a, the results using R max = 30 km match the trend of the observed sea-level pressures closer than the ones using R max = 50 km. However, both of them underestimate the peak value of the sea-level pressures measured at the Guiuan weather station. The comparison against the observed winds is illustrated in Figure 4b

Haiyan's Storm Surge Simulation by HWM
After examining the 10-m winds and sea-level pressures against the observation data at Guiuan weather station, the storm surge modeling driven by HWM are explored in this section.
In storm surge modeling, winds can amplify the surge heights in coastal regions when wind directions are more or less perpendicular to coasts. Besides, the intensity of winds is relative to the wind shear stresses over a water surface. Hence, before performing storm surge predictions, the wind fields on the standard level of 10 m shall be examined. In Figure 5, the snapshots of the 10-m winds predicted by HWM are presented. In Figure 5a, the storm center of Haiyan just enters the Leyte Gulf and passes over the small island in the south of the Samar Island. In a while, the value of 10-m winds blowing to the Samar Island is higher than 60 m/s, and the strong offshore winds occur in the San Pedro Bay. As shown in Figure 5b, Haiyan moves to the central regions of the Leyte Gulf, and the HWM-predicted winds hit the southern regions of the Samar Island, such as Giporlos and Balangiga. When Haiyan makes landfall at the Leyte Island, the severe onshore winds towards the San Pedro Bay can be found in Figure 5c. The winds of more than 60 m/s hit the regions in particular Tacloban, Basey, and Palo. After Haiyan makes landfall at Leyte Island, the southeastern winds keep blowing from the Leyte Gulf to the San Pedro Bay (see Figure 5d). As shown in Figure 5d, after the landfall, the cyclonic winds still maintain the wind speeds of more than 45 m/s. Despite using the 1954 Schloemer's scaling parameter with R max = 30 km seems to have a better comparison with the observation data in terms of the sea-level pressures, the predicted winds cannot match the dramatic increase from 15:00 UTC 7 November 2013. By verifying against the observation data at the Guiuan weather station, the use of the pressure distribution scaled by the 1999 Harper and Holland's parameters with R max = 30 km has the best agreement in the wind comparison. Hence, the results of HWM presented in the following sections will use this configuration.

Haiyan's Storm Surge Simulation by HWM
After examining the 10-m winds and sea-level pressures against the observation data at Guiuan weather station, the storm surge modeling driven by HWM are explored in this section.
In storm surge modeling, winds can amplify the surge heights in coastal regions when wind directions are more or less perpendicular to coasts. Besides, the intensity of winds is relative to the wind shear stresses over a water surface. Hence, before performing storm surge predictions, the wind fields on the standard level of 10 m shall be examined. In Figure 5, the snapshots of the 10-m winds predicted by HWM are presented. In Figure 5a, the storm center of Haiyan just enters the Leyte Gulf and passes over the small island in the south of the Samar Island. In a while, the value of 10-m winds blowing to the Samar Island is higher than 60 m/s, and the strong offshore winds occur in the San Pedro Bay. As shown in Figure 5b, Haiyan moves to the central regions of the Leyte Gulf, and the HWM-predicted winds hit the southern regions of the Samar Island, such as Giporlos and Balangiga. When Haiyan makes landfall at the Leyte Island, the severe onshore winds towards the San Pedro Bay can be found in Figure 5c. The winds of more than 60 m/s hit the regions in particular Tacloban, Basey, and Palo. After Haiyan makes landfall at Leyte Island, the southeastern winds keep blowing from the Leyte Gulf to the San Pedro Bay (see Figure 5d). As shown in Figure 5d, after the landfall, the cyclonic winds still maintain the wind speeds of more than 45 m/s. In Figure 6, the computed storm surges driven by HWM are illustrated. Due to the strong offshore winds in the San Pedro Bay and the Leyte Gulf, the negative storm surges are generated accordingly (see Figure 6a,b). As shown in Figure 6c, when Haiyan makes landfall at the Leyte Island, the cyclonic winds push up the water hitting the regions such as Tacloban, Palo, and Tanauan. Storm surges of more than 3.5 m are observed in the simulation. After Haiyan makes landfall at the Leyte Island, the winds in the San Pedro Bay become southeastern (see Figure 6d). Because of these winds, the waters are piled up inside the bay and generate the significant storm surges of more than 4.0 m, which hit Tacloban, Palo, and Basey. In Figure 6, the computed storm surges driven by HWM are illustrated. Due to the strong offshore winds in the San Pedro Bay and the Leyte Gulf, the negative storm surges are generated accordingly (see Figure 6a,b). As shown in Figure 6c, when Haiyan makes landfall at the Leyte Island, the cyclonic winds push up the water hitting the regions such as Tacloban, Palo, and Tanauan. Storm surges of more than 3.5 m are observed in the simulation. After Haiyan makes landfall at the Leyte Island, the winds in the San Pedro Bay become southeastern (see Figure 6d). Because of these winds, the waters are piled up inside the bay and generate the significant storm surges of more than 4.0 m, which hit Tacloban, Palo, and Basey. The maximum computed storm surges driven are shown in Figure 7, and the maximum computed surge heights along the coasts of the Leyte and Samar Islands are compared to the measured field survey data as well. In Figure 7a,b, the field survey data of Tajima et al. (2014) [33] with the Reliability A (clear mark and small error of survey) are presented and the tide effect has been removed. The profiles of maximum computed surge heights are acquired one grid away from the coasts to avoid the boundary effect. In Figure 7a,b, the simulation results agree well with the patterns of the measured water levels along the coasts of the Leyte Island and the Samar Island generally. Besides, relatively higher storm surges are found near the inner of the San Pedro Bay in both model results and the measured water levels. However, the measured water levels of more than 6 m can be found near the innermost part of the San Pedro Bay (see Figure 7a), where the wave overtopping and the up-wash of wind waves could be severe [19]. The model results may give an underestimation of water levels near the bay. Moreover, the comparisons are expected to be better if the inundation calculation can be performed with the detailed land elevation model and the shape of coastline in the Leyte Gulf and the San Pedro Bay [18]. As shown in Figure 7c, inside the San Pedro Bay, the maximum storm surges of more than 4.5 m are found in Basey, Tacloban, Palo, and Tanauan. These prominent maximum storm surges are attributed to the southeastern winds from the Leyte Gulf to the San Pedro Bay, making the water pile up inside the bay after Haiyan makes landfall at the Leyte Island (see Figure 6d). The maximum computed storm surges driven are shown in Figure 7, and the maximum computed surge heights along the coasts of the Leyte and Samar Islands are compared to the measured field survey data as well. In Figure 7a,b, the field survey data of Tajima et al. (2014) [33] with the Reliability A (clear mark and small error of survey) are presented and the tide effect has been removed. The profiles of maximum computed surge heights are acquired one grid away from the coasts to avoid the boundary effect. In Figure 7a,b, the simulation results agree well with the patterns of the measured water levels along the coasts of the Leyte Island and the Samar Island generally. Besides, relatively higher storm surges are found near the inner of the San Pedro Bay in both model results and the measured water levels. However, the measured water levels of more than 6 m can be found near the innermost part of the San Pedro Bay (see Figure 7a), where the wave overtopping and the up-wash of wind waves could be severe [19]. The model results may give an underestimation of water levels near the bay. Moreover, the comparisons are expected to be better if the inundation calculation can be performed with the detailed land elevation model and the shape of coastline in the Leyte Gulf and the San Pedro Bay [18]. As shown in Figure 7c, inside the San Pedro Bay, the maximum storm surges of more than 4.5 m are found in Basey, Tacloban, Palo, and Tanauan. These prominent maximum storm surges are attributed to the southeastern winds from the Leyte Gulf to the San Pedro Bay, making the water pile up inside the bay after Haiyan makes landfall at the Leyte Island (see Figure 6d). In Figure 8, the time series of computed storm surges at numerical gauges are presented. In all stations, the records show the declines of the water-level first and then following quick rises within 3 h. This phenomenon corresponds to the witness reports about the "tsunami-like" wavefront hit Tacloban and Basey after the water-level receded [32,35]. Maximum storm surges of more than 4.0 m are found in Tacloban, Basey, Tanauan, and Palo (see Figure 8a,d,e,h), which corresponds to the maximum computed surges and field survey shown in Figure 7. There are two water-level measurements at Tacloban and Guiuan tidal stations during the passage of Haiyan [19]. The observation of the Tacloban tidal station was interrupted before the peak storm surge occurred; hence, this observation didn't record the complete time series of water-level. In the Guiuan tide station which is located at the southeastern part of the Samar Island not the San Pedro Bay, the wind wave and storm surge interactions could be significant. Hence, we only use the field survey data of Tajima   In Figure 8, the time series of computed storm surges at numerical gauges are presented. In all stations, the records show the declines of the water-level first and then following quick rises within 3 h. This phenomenon corresponds to the witness reports about the "tsunami-like" wavefront hit Tacloban and Basey after the water-level receded [32,35]. Maximum storm surges of more than 4.0 m are found in Tacloban, Basey, Tanauan, and Palo (see Figure 8a,d,e,h), which corresponds to the maximum computed surges and field survey shown in Figure 7. There are two water-level measurements at Tacloban and Guiuan tidal stations during the passage of Haiyan [19]. The observation of the Tacloban tidal station was interrupted before the peak storm surge occurred; hence, this observation didn't record the complete time series of water-level. In the Guiuan tide station which is located at the southeastern part of the Samar Island not the San Pedro Bay, the wind wave and storm surge interactions could be significant. Hence, we only use the field survey data of Tajima  The information of numerical gauges can be found in Figure 2b and Table 2. The x-axis indicates the time from 12:00 UTC 7 November 2013 to 12:00 UTC 8 November 2013.

Haiyan's Storm Surge Simulation by WRF-ARW
In this section, the results of storm surge modeling driven by WRF-ARW are explored. The 10m winds and sea-level pressure from WRF-ARW 3 km and WRF-ARW 1 km with the 1-h time interval are used to drive the storm surge model. To reach the storm intensity of Haiyan in the WRF-ARW modeling, the nudging scheme, four-dimensional data assimilation (FDDA), is used from 00:00 UTC 4 November 2013 to 00:00 UTC 5 November 2013 (24 h). Afterward, the simulation is carried out by WRF-ARW only without the further use of data assimilation. Hence, the translation speed of Haiyan predicted by WRF-ARW is slower than the JMA best-track data (see Figure 9). As illustrated in Figure  9, the storm landfall locations between the JMA best-track and WRF-ARW 3 km are more or less the same and WRF-ARW 1 km predicted more south landfall location at the Leyte Island. The landfall times predicted by WRF-ARW 3 km and WRF-ARW 1 km at the Leyte Island are 4-5 h late than the best-track. The storm centers of WRF-ARW 3 km and WRF-ARW 1 km and their distance from the best-track are tabulated in Table 3. More discussions about the WRF-ARW modeling can be found in Kueh et al. (2019) [15]. By the deviations of landfall locations and the translation speed of Haiyan from WRF-ARW, the results of storm surge prediction driven by WRF-ARW are anticipated to show different patterns from those driven by HWM. The information of numerical gauges can be found in Figure 2b and Table 2. The x-axis indicates the time from 12:00 UTC 7 November 2013 to 12:00 UTC 8 November 2013.

Haiyan's Storm Surge Simulation by WRF-ARW
In this section, the results of storm surge modeling driven by WRF-ARW are explored. The 10-m winds and sea-level pressure from WRF-ARW 3 km and WRF-ARW 1 km with the 1-h time interval are used to drive the storm surge model. To reach the storm intensity of Haiyan in the WRF-ARW modeling, the nudging scheme, four-dimensional data assimilation (FDDA), is used from 00:00 UTC 4 November 2013 to 00:00 UTC 5 November 2013 (24 h). Afterward, the simulation is carried out by WRF-ARW only without the further use of data assimilation. Hence, the translation speed of Haiyan predicted by WRF-ARW is slower than the JMA best-track data (see Figure 9). As illustrated in Figure 9, the storm landfall locations between the JMA best-track and WRF-ARW 3 km are more or less the same and WRF-ARW 1 km predicted more south landfall location at the Leyte Island. The landfall times predicted by WRF-ARW 3 km and WRF-ARW 1 km at the Leyte Island are 4-5 h late than the best-track. The storm centers of WRF-ARW 3 km and WRF-ARW 1 km and their distance from the best-track are tabulated in Table 3. More discussions about the WRF-ARW modeling can be found in Kueh et al. (2019) [15]. By the deviations of landfall locations and the translation speed of Haiyan from WRF-ARW, the results of storm surge prediction driven by WRF-ARW are anticipated to show different patterns from those driven by HWM.  In Figures 10 and 11, the snapshots of 10-m winds of WRF-ARW 3 km and WRF-ARW 1 km are presented respectively. As shown in Figure 10a,b, the storm center of Haiyan predicted by WRF-ARW 3 km just arrives at the small island in the south of East Samar and then passes over it. The locations of the storm center are similar to the JMA's best track (see Figure 6a). At these moments, the winds hit southeastern regions of the Samar Island and offshore winds occur in the San Pedro Bay. In Figure 10c, the storm center predicted by WRF-ARW 3 km tends to land near Mayorga, which is relatively south than the one predicted by HWM (see Figure 6c). After making landfall at the Leyte Island, the southeastern winds occur in the San Pedro Bay (see Figure 6d); however, the wind intensity is not as strong as HWM.
In Figure 11a, the storm center of Haiyan predicted by WRF-ARW 1 km passes over Dinagat islands, which is further south than the ones predicted by WRF-ARW 3 km and the JMA best track. Compared to wind speeds of WRF-ARW 3 km, WRF-ARW 1 km's results are much more potent, about 20% stronger than that. As shown in Figure 11b, the winds of WRF-ARW 1 km are more parallel to the southern regions of the Samar Island. The wind speeds of more than 55 m/s are found in the east and northeast part to the storm center. WRF-ARW 1 km gives a further south landfall location than others; hence, the strong onshore winds hit the regions such as Macarthur, Dulag, and Tolosa (see Figure 11c). After making landfall at the Leyte Island, the southeastern winds blowing from the Leyte Gulf to the San Pedro Bay occur (see Figure 11d). Moreover, the stronger winds of WRF-ARW 1 km hit the regions such as Tolosa and Dulag and the maximum 10-m wind speed of more than 45 m/s is still found near the east coasts of the Leyte Island.  In Figures 10 and 11, the snapshots of 10-m winds of WRF-ARW 3 km and WRF-ARW 1 km are presented respectively. As shown in Figure 10a,b, the storm center of Haiyan predicted by WRF-ARW 3 km just arrives at the small island in the south of East Samar and then passes over it. The locations of the storm center are similar to the JMA's best track (see Figure 6a). At these moments, the winds hit southeastern regions of the Samar Island and offshore winds occur in the San Pedro Bay. In Figure 10c, the storm center predicted by WRF-ARW 3 km tends to land near Mayorga, which is relatively south than the one predicted by HWM (see Figure 6c). After making landfall at the Leyte Island, the southeastern winds occur in the San Pedro Bay (see Figure 6d); however, the wind intensity is not as strong as HWM.
Water 2020, 12, x FOR PEER REVIEW 17 of 27 In the comparison between WRF-ARW 3 km and WRF-ARW 1 km, the significant discrepancy is the storm track, which affects the predicted landfall location at the Leyte Island and the southeastern winds blowing from the Leyte Gulf to the San Pedro Bay.    In the comparison between WRF-ARW 3 km and WRF-ARW 1 km, the significant discrepancy is the storm track, which affects the predicted landfall location at the Leyte Island and the southeastern winds blowing from the Leyte Gulf to the San Pedro Bay.   In Figure 11a, the storm center of Haiyan predicted by WRF-ARW 1 km passes over Dinagat islands, which is further south than the ones predicted by WRF-ARW 3 km and the JMA best track. Compared to wind speeds of WRF-ARW 3 km, WRF-ARW 1 km's results are much more potent, about 20% stronger than that. As shown in Figure 11b, the winds of WRF-ARW 1 km are more parallel to the southern regions of the Samar Island. The wind speeds of more than 55 m/s are found in the east and northeast part to the storm center. WRF-ARW 1 km gives a further south landfall location than others; hence, the strong onshore winds hit the regions such as Macarthur, Dulag, and Tolosa (see Figure 11c). After making landfall at the Leyte Island, the southeastern winds blowing from the Leyte Gulf to the San Pedro Bay occur (see Figure 11d). Moreover, the stronger winds of WRF-ARW 1 km hit the regions such as Tolosa and Dulag and the maximum 10-m wind speed of more than 45 m/s is still found near the east coasts of the Leyte Island.
In the comparison between WRF-ARW 3 km and WRF-ARW 1 km, the significant discrepancy is the storm track, which affects the predicted landfall location at the Leyte Island and the southeastern winds blowing from the Leyte Gulf to the San Pedro Bay.
The computed storm surges driven by WRF-ARW 3 km and WRF-ARW 1 km are illustrated in Figures 12 and 13, respectively. As shown in Figure 12a,b, due to the cyclonic winds, the negative storm surges occur in the Leyte Gulf and the San Pedro Bay. When the storm center of Haiyan approaches the Leyte Island, the positive storm surges of more than 1 m are generated (see Figure 12c). In Figure 12d, after Haiyan makes landfall at the Leyte Island, the storm surges of more than 2.0 m pile up in the San Pedro Bay and hit Tacloban, Palo, and Tanauan. The pattern of computed storm surges is similar to that predicted by HWM (see Figure 6d) but lower. In the results driven by WRF-ARW 1 km, because of the more south predicted storm track, the storm surges of more than 2.0 m first hit Dulag and Macarthur (see Figure 13b,c). After Haiyan makes landfall at the Leyte Island, the storm surges of more than 2.0 m occur in the San Pedro Bay (see Figure 13d). The computed storm surges driven by WRF-ARW 3 km and WRF-ARW 1 km are illustrated in Figures 12 and 13, respectively. As shown in Figure 12a,b, due to the cyclonic winds, the negative storm surges occur in the Leyte Gulf and the San Pedro Bay. When the storm center of Haiyan approaches the Leyte Island, the positive storm surges of more than 1 m are generated (see Figure  12c). In Figure 12d, after Haiyan makes landfall at the Leyte Island, the storm surges of more than 2.0 m pile up in the San Pedro Bay and hit Tacloban, Palo, and Tanauan. The pattern of computed storm surges is similar to that predicted by HWM (see Figure 6d) but lower. In the results driven by WRF-ARW 1 km, because of the more south predicted storm track, the storm surges of more than 2.0 m first hit Dulag and Macarthur (see Figure 13b,c). After Haiyan makes landfall at the Leyte Island, the storm surges of more than 2.0 m occur in the San Pedro Bay (see Figure 13d).  In Figures 14 and 15, the maximum computed storm surges driven by WRF-ARW 3 km and WRF-ARW 1 km are illustrated and the maximum computed surge heights along the coasts of the Leyte Island and the Samar Island are compared to the measured field survey data [33]. In Figure 14c, the maximum storm surges of more than 2 m driven by WRF-ARW 3 km are found in Tacloban, Palo, Tanauan, and Basey, which of them are inside the San Pedro Bay. The pattern is similar to the one driven by HWM; however, the amplitude of maximum computed storm surges inside the San Pedro Bay is almost twice lower. Since lower maximum computed storm surges are predicted, the results cannot match perfectly with the field survey data (see Figure 14a,b). In the result driven by WRF-ARW 1 km, the prominent computed maximum storm surges of more than 2 m are also found on the east coast of Leyte Island such as Dulag, Macarthur, and Abuyog (see Figure 15b). As expected, the results cannot match well with the field survey data as well (see Figure 15a,c). In Figures 14 and 15, the maximum computed storm surges driven by WRF-ARW 3 km and WRF-ARW 1 km are illustrated and the maximum computed surge heights along the coasts of the Leyte Island and the Samar Island are compared to the measured field survey data [33]. In Figure 14c, the maximum storm surges of more than 2 m driven by WRF-ARW 3 km are found in Tacloban, Palo, Tanauan, and Basey, which of them are inside the San Pedro Bay. The pattern is similar to the one driven by HWM; however, the amplitude of maximum computed storm surges inside the San Pedro Bay is almost twice lower. Since lower maximum computed storm surges are predicted, the results cannot match perfectly with the field survey data (see Figure 14a,b). In the result driven by WRF-ARW 1 km, the prominent computed maximum storm surges of more than 2 m are also found on the east coast of Leyte Island such as Dulag, Macarthur, and Abuyog (see Figure 15b). As expected, the results cannot match well with the field survey data as well (see Figure 15a,c).
In Figure 16, the time series of computed storm surges in the numerical gauges are shown. Since both of WRF-ARW 3 km and WRF-ARW 1 km predict slower storm forward motion than the best track, the arrival time of peak storm surges occur behind those in the HWM's results. In the stations in the inner of the San Pedro Bay such as Tacloban and Basey, the results driven by WRF-ARW 3 km and WRF-ARW 1 km have the similar patterns. In these stations, the peak predicted storm surges driven by WRF-ARW 3 km are slightly higher than those by WRF-ARW 1 km; however, both of them are still lower than those driven by HWM. For stations such as Abuyog, Macarthur, and Dulag, because WRF-ARW 1 km predicts the further south landfall location at the Leyte Island, the winds blow toward these regions and cause storm surges about 0.8-1.0 m higher than those in WRF-ARW 3 km. Moreover, in these stations, the differences of maximum storm surges driven between HWM and WRF-ARW are not as prominent as those located in the San Pedro Bay.     driven by WRF-ARW 3 km are slightly higher than those by WRF-ARW 1 km; however, both of them are still lower than those driven by HWM. For stations such as Abuyog, Macarthur, and Dulag, because WRF-ARW 1 km predicts the further south landfall location at the Leyte Island, the winds blow toward these regions and cause storm surges about 0.8-1.0 m higher than those in WRF-ARW 3 km. Moreover, in these stations, the differences of maximum storm surges driven between HWM and WRF-ARW are not as prominent as those located in the San Pedro Bay.  Figure 2b and Storm surge modeling driven by both WRF-ARW 3 km and WRF-ARW 1 km cannot reproduce obvious storm surges in the San Pedro Bay. For WRF-ARW 3 km, the predicted storm track is close to the JMA's best track; however, WRF-ARW 3 km underestimates the wind intensity. Despite WRF-ARW 1 km gives comparable higher wind intensity, it predicts relatively south storm track of Haiyan, which makes the underestimation of storm surges in the San Pedro Bay. Moreover, both of WRF-ARW 3 km and WRF-ARW 1 km predict slower storm forward motion, which also affects the storm surge calculation. Next, the effects of the storm track and the storm forward speed on generating storm surges in the San Pedro Bay and the Leyte Gulf will be explored. Storm surge modeling driven by both WRF-ARW 3 km and WRF-ARW 1 km cannot reproduce obvious storm surges in the San Pedro Bay. For WRF-ARW 3 km, the predicted storm track is close to the JMA's best track; however, WRF-ARW 3 km underestimates the wind intensity. Despite WRF-ARW 1 km gives comparable higher wind intensity, it predicts relatively south storm track of Haiyan, which makes the underestimation of storm surges in the San Pedro Bay. Moreover, both of WRF-ARW 3 km and WRF-ARW 1 km predict slower storm forward motion, which also affects the storm surge calculation. Next, the effects of the storm track and the storm forward speed on generating storm surges in the San Pedro Bay and the Leyte Gulf will be explored.

Experiment for Storm Forward Speed
In previous sections, the storm surge modeling driven by HWM, WRF-ARW 3 km, and WRF-ARW 1 km presents the prominent differences due to the deviations of predicted storm tracks and the forward speed of storm. Tajima et al. (2016) [38] and Kumagai et al. (2016) [19] have conducted a series of sensitivity analysis on the variations of storm track and the storm forward speed to the storm surges in the San Pedro Bay. Tajima et al. (2016) [38] found that Haiyan's storm track seems to be the worse one for the regions near Palo and Tanauan. Besides, they pointed out that maximum storm surge at Tacloban will be reduced by about 45% if the forward speed of Haiyan is decreased by 50%. Kumagai et al. (2016) [19] indicated that the storm surges near Dulag will increase if Haiyan's storm track become more south, but storm surges in the San Pedro Bay will reduce. Both of these studies correspond to our findings in the storm surge modeling driven by HWM, WRF-ARW 3 km, and WRF-ARW 1 km.
In this section, the numerical experiment focused on a storm forward speed is carried out since only few stations have been explored [38] and a constant forward speed was assumed [19]. Because Haiyan's storm surge simulation affected by the storm-track variations [19,38] and the radius of maximum winds have been explored carefully [18], these discussions will not be included here. In this numerical experiment, besides the original storm forward speed, this study includes two comparable cases: forward speed increased by 10% and decreased by 10%. The 10-m winds and the sea-level pressures in this experiment are from HWM. To analyze the hydrodynamical response of storm surges to a storm forward speed, the wind asymmetry modified by a storm forward motion is considered by the original forward speed for the sake of simplicity.
The time series of storm surges in the numerical gauges are illustrated in Figure 17. In the stations located inside the San Pedro Bay such as Tacloban, Basey, Palo, and Tanauan, the maximum computed storm surges are intensified by a faster storm forward speed (see Figure 17a,d,e,h). On the other hand, by considering the slower storm forward speed, the maximum computed storm surges are reduced. However, in the stations outside the San Pedro Bay such as Balangiga and Abuyog, the maximum computed storm surges are reduced in the case with the faster storm forward speed (see Figure 17b,g). The statistical results of the maximum computed storm surges are tabulated in Table 4. From Table 4, considering the storm forward speed increased by 10%, the extreme values of storm surges in the stations inside the San Pedro Bay such as Tacloban and Tanauan are intensified by 4.5-7.0%; however, the reductions of 5.3-6.8% on the maximum values are shown when considering the slower forward speed.
Water 2020, 12, x FOR PEER REVIEW 23 of 27 Figure 17. (a-i) Time series of the computed storm surges in the numerical experiment on a storm forward motion: original speed (black lines), increased forward speed (red lines), and decreased forward speed (blue lines) in the numerical gauges. The detailed information of numerical gauges can be found in Figure 2b and Table 2. The x-axis indicates the time from 12:00 UTC 7 November 2013 to 12:00 UTC 8 November 2013.

Conclusions
This study has explored the discrepancies of the storm surge predictions driven by a parametric wind (PW) model and a numerical weather prediction (NWP) model. By shedding light upon the discrepancies, the storm surges of 2013 Typhoon Haiyan in the Leyte Gulf and the San Pedro Bay have been studied. The frictional-free, steady-state, and geostrophic-balancing solutions on the storm winds and sea-level pressures given by the Holland wind model (HWM) have been first verified Figure 17. (a-i) Time series of the computed storm surges in the numerical experiment on a storm forward motion: original speed (black lines), increased forward speed (red lines), and decreased forward speed (blue lines) in the numerical gauges. The detailed information of numerical gauges can be found in Figure 2b and Table 2. The x-axis indicates the time from 12:00 UTC 7 November 2013 to 12:00 UTC 8 November 2013.

Conclusions
This study has explored the discrepancies of the storm surge predictions driven by a parametric wind (PW) model and a numerical weather prediction (NWP) model. By shedding light upon the discrepancies, the storm surges of 2013 Typhoon Haiyan in the Leyte Gulf and the San Pedro Bay have been studied. The frictional-free, steady-state, and geostrophic-balancing solutions on the storm winds and sea-level pressures given by the Holland wind model (HWM) have been first verified against the observed data at the Guiuan weather station. The 10-m winds predicted by WRF-ARW (Weather Research and Forecasting-Advanced Research WRF), which solves the 3D time-integrated, compressible, and non-hydrostatic Euler equations, have been compared with the ones by HWM. In the results of computed storm surges driven by HWM, the storm surges of more than 4.5 m have been found in the San Pedro Bay due to the water piling up and the pattern is agreed well with the field survey data. However, in the results by both WRF-ARW 3 km and WRF-ARW 1 km, the maximum storm surge in the San Pedro Bay is almost half than that driven by HWM. The significant differences among HWM, WRF-ARW 3 km, and WRF-ARW 1 km are caused by the predicted landfall location of Haiyan, the storm intensity, and the storm forward speed. The numerical experiment focused on the storm forward speed has pointed out that the extreme values of storm surges inside the San Pedro Bay can be intensified by 4.5-7.0% when the storm forward speed is increased by 10%. This study hopes to pave the way for future research in understanding the discrepancies in predicting storm surges driven by either of the PW or NWP models and thus benefiting storm surge predictions in a semi-enclosed bay.