Parallel-Computing Two-Way Grid-Nested Storm Surge Model with a Moving Boundary Scheme and Case Study of the 2013 Super Typhoon Haiyan

This study presents a numerical tool for calculating storm surges from offshore, nearshore, and coastal regions using the finite-difference method, two-way grid-nesting function in time and space, and a moving boundary scheme without any numerical filter adopted. The validation of the solitary wave runup on a circular island showed the perfect matches between the model results and measurements for the free surface elevations and runup heights. After the benchmark problem validation, the 2013 Super Typhoon Haiyan event was selected to showcase the storm surge calculations with coastal inundation and flood depths in Tacloban. The catastrophic storm surges of about 8 m and wider, storm-induced inundation due to the Super Typhoon Haiyan were found in the Tacloban Airport, corresponding to the findings from the field survey. In addition, the anti-clockwise, storm-induced currents were explored inside of Cancabato Bay. Moreover, the effect of the nonlinear advection terms with the fixed and moving shoreline and the parallel efficiency were investigated. By presenting a storm surge model for calculating storm surges, inundation areas, and flood depths with the model validation and case study, this study hopes to provide a convenient and efficient numerical tool for forecasting and disaster assessment under a potential severe tropical storm with climate change.


Introduction
Storm surges due to a tropical cyclone cause a catastrophic impact on coastal communities [1]. Due to climate change resulting in the occurrence of more severe typhoons [2,3], having a better understanding of storm surges is required. However, storm surges involve multiple physical factors, thus complicating hydrodynamical modeling. First, storm surges interacting with tides and wind waves may intensify coastal surge heights [4][5][6][7][8]. Second, when a tropical cyclone closes coasts, surface winds amplify the amplitude of storm surges in shallow waters [7,8]. Lastly, the Coriolis force with a particular environment will generate unignored forerunner surges before main surges [9]. Hence, to counter the coastal storm surge impact with the multiple physical factors, an accurate numerical tool for evaluating the threat of storm surges is in high demand.
From a numerical point of view, storm surge models can be divided into two groups: (1) structured-grid model; and (2) unstructured-grid model. Structured-grid storm surge models are easily implemented and developed for operational purposes in the early stage.
An example of a structured-grid model is SLOSH (Sea, Lake, and Overland Surges from Hurricanes [10]). However, these kinds of models may not handle the change of the wavelength of storm surges perfectly because of limitations of a grid size [11], ignoring the advection/horizontal eddy diffusion terms on simulating coastal storm surges [10], or simulating storm surges without inundation areas and flood depths [10,11]. Thus, the gridnesting function becomes a good option for a structured-grid storm surge model to consider better physics and descriptions for calculating coastal storm surges. Using the gridnesting function, a structured-grid storm surge model can simulate storm surges from the ocean to coasts with the appropriate grid sizes for both deep and shallow waters. For example, some of these nested-grid models are NCTSM (Nested Coupled Tide-Surge Model [11]) and SuWAT (Surge, WAve, and Tide [12]). Besides using the grid-nesting function, structured models adopting a curvilinear grid improve the simulation of coastal storm surges around a complicated coastline, such as CH3D-SSMS (Curvilinear Hydrodynamic in 3D Storm Surge Modeling System [13]) and ROMS (Regional Ocean Modeling System [14]). Here, it is noted that a curvilinear-grid storm surge model sometimes adopts a relatively small domain; thus, boundary conditions need to be provided by a basin/regionalscale model [13]. In addition to structured-grid models, unstructured models have recently become more prevalent in simulating storm surges. Most unstructured-grid storm surge models are extended from ocean current models, allowing for the simulation of storm surges in only one computational grid with a more extensive range [7,8] or resolving the three-dimensional structure of storm-induced currents [15]. These models are, for example, SCHISM (Semi-implicit Cross-scale Hydroscience Integrated System Model [16]), ADCIRC (ADvanced CIRCulation [7,8]), and FVCOM (Finite Volume Community Ocean Model [15]). Although an unstructured-grid storm surge model allows the computational grid to gradually change from coarse to fine, it needs to shoulder the convergence and stability issue of the grid system [17], which usually relies on professional grid-maker software. This implies that the grid size of an unstructured grid is difficult to modify after the grid has been built, especially for grids near the coastline [17]. Hence, by comparing the pros and cons between the structured-grid and unstructured-grid models, the structured-grid models with the grid-nesting function should be a user-friendly choice implemented for users and developers to study storm surges from deep to shallow waters.
From a practical point of view, a storm surge model should be able to evaluate future storm surge hazards [18] or predict storm surges in a deterministic/ensemble manner before the landfall of a tropical cyclone [19,20]. First, since climate change has caused the sea-level to rise and these water levels have become a severe issue in the upcoming future [2,3], a hazard map for storm-induced inundations and floods is vital to coastal communities [21][22][23]. However, calculating storm-induced inundations and floods requires an accurate moving boundary scheme (i.e., moving shoreline scheme) with a stable nonlinear advection (convection) solver of a storm surge model, which usually results in some numerical issues during computations [24] and is usually ignored in the operational forecasting [10]. Thus, storm surge calculations with the moving boundary scheme tracing flood depth and inundation due to a tropical cyclone are challenging. Second, the prediction of coastal storm surge heights and warning of potential storm surges before a tropical cyclone landfall are essential for a global or local operational organization; such predictive simulations include Storm Surge Maximum of the Maximum (MOM) [11] and Maximum Envelope of High Water (MEOW) [11,25]. Hence, an operational storm surge model shall simulate storm surges from the basin/regional scale to the coastal scale. In addition, the model shall have enough efficiency to conduct ensemble simulations (such as MOM and MEOW) and update the warning messages appropriately. Moreover, predicting storminduced inundation areas and flood depths is more important than solely calculating coastal storm surge heights [10,11]. In summary, a storm surge model should have the ability to (1) evaluate potential storm surge hazard maps by simulating coastal storm surges with inundations and flood depths and (2) predict storm surges from offshore to nearshore with flood depths and inundation areas with enough operational efficiency.
After reviewing the structured/unstructured storm surge model development and practical uses of the storm surge predictions, a well-developed storm surge model should satisfy accuracy, stability, efficiency, flexibility, and convenience when performing storm surge computations from offshore to nearshore. Thus, this paper aims to present a numerical tool for calculating storm surges satisfying the aforementioned requirements. The storm surge model developed in this study is extended from the well-known tsunami model, COMCOT (COrnell Multi-grid Coupled Tsunami Model) [26,27] and allows the two-way grid-nesting function to increase the spatial resolution of nearshore regions. In addition, the moving boundary scheme with the nonlinear advection-term solver is expected to handle the climbing of coastal storm surges. Moreover, the performance and efficiency are enhanced by the parallel-computing technique. Here, we note that the linear equation model presented in Tsai et al. (2020) [28] is also based on COMCOT. However, that model ignores the nonlinear advection terms, horizontal eddy diffusions, moving boundary scheme, and grid-nesting function. Thus, the model presented in this study can be considered as an extension of the 2020's model from both physical and numerical points of view. Furthermore, some other studies have also conducted storm surge modeling based on COMCOT [29,30]. Yet, they did not involve any parallel-computing function in the simulations.
This paper is organized as follows. In Section 2, we will first present the storm surge model used in this study from the governing equations, discretization, grid-nesting/moving boundary scheme, to the parallel-computing technique. Next, in Section 3, this study will demonstrate a benchmark problem for validating the moving boundary scheme and the advection-term solver. Afterward, in Section 4, we will show the 2013 Super Typhoon Haiyan event with storm surges, storm-induced currents, and inundations around Leyte Gulf and San Pedro Bay. Finally, in Section 5, we will conclude this study and illustrate some work for the future.

Methodology
This study carries out storm surge computation using a depth-integrated equation model with the Coriolis effect, bottom friction, and horizontal eddy diffusion [11,12]. This storm surge model is extended from the well-developed COMCOT tsunami model to have better simulation accuracy. The reasons for developing COMCOT from simulating tsunamis to storm surges are as follows: (1) COMCOT has been validated by several benchmark problems [27,31] and real tsunami events, such as the 2011 Japan Tohoku Tsunami [32] and 2019 Indonesia Palu Tsunami [33]; thus, the accuracy of the numerical algorithm such as the upwind advection solver has been proved. (2) Both tsunamis and storm surges are in the regime of long waves; the physical equation can easily handle the computation change from tsunami to storm surge by adding some physical terms (e.g., wind shear stress terms) with slight modifications. (3) COMCOT has a two-way grid-nesting function in time and space, allowing seamless storm surge simulations from offshore to nearshore with flexible grid sizes/time steps and more accurate model results. (4) The parallel-computing function has been added to COMCOT for a fast-calculation purpose; hence, this parallel-computing function should be easily included into storm surge calculations [32]. From the reasons mentioned above, the storm surge model developed in this study will be nicknamed COMCOT-SURGE (COrnell Multi-grid COupled Tsunami-Storm SURGE). The governing equations, discretization, two-way grid-nesting, moving boundary scheme, and parallel-computing function of COMCOT-SURGE will be illustrated in the following subsections.

Governing Equation of the Storm Surge Model
The mass and momentum equations of COMCOT-SURGE are presented as follows: where is the free surface elevation (unit: m), is the time (unit: s), ( , ) are the spatial notations of the x-and y-directions (unit: m), ( , ) are the volume-flux components (i.e., depth-integrated volume fluxes per unit length) of the x-and y-directions (unit: m 2 /s), is the total water depth ( = ℎ + ; unit: m), ℎ is the still water depth (unit: m), g is the gravitational acceleration (= 9.81 m/s 2 ), is the Coriolis parameter ( = 2 ; unit: 1/s), is the Earth angular velocity (= 7.2921 × 10 −5 rad/s), is the sea-level air pressure (unit: N/m 2 ), ( , ) are the wind shear stresses (unit: N/m 2 ), ( , ) are the bottom frictional shear stresses (unit: N/m 2 ), is the water density (unit: kg/m 3 ), and ℎ is the horizontal eddy diffusion coefficient (unit: m 2 /s). Here we note that Equations (2) and (3) ignore the advection terms where the nonlinear effect becomes insignificant in deep waters. Thus, the linear momentum equations are shown below: The Manning's formula, from the conception of the open-channel flow, is used to model the bottom friction: where is the Manning's roughness coefficient (unit: s/m 1/3 ). The quadric law is used to model the wind shear stresses on the water surface: in which is the wind-drag coefficient between the water surface and the air, 10 is the 10-m wind speed (unit: m/s), ( 10 , 10 ) are the components of the 10-m wind speed in the x-and y-directions. The wind-drag coefficient proposed by Wu (1982) [34] with an imposed lower bound limit of WAMDI (1988) [35] is expressed; in addition, the cap of the wind-drag coefficient is determined [12,15] Here we note that is a dimensionless coefficient.

Discretization
COMCOT-SURGE adopts the finite-difference method to discretize the mass and momentum equations. Due to the leap-frog scheme, the mass and momentum equations are solved at a separate time step. Some studies have used a similar procedure to discretize their governing equations on storm surge modeling (e.g., Kim et al., 2015 [12]). The discretized mass equation is shown below: where ∆ is the time step (unit: s) and (∆ , ∆ ) are the grid sizes in the x-and y-directions (unit: m). The momentum equations are discretized using the upwind scheme for advection terms, implicit form for bottom friction terms, and central difference form for horizontal eddy diffusions. Thus, the consequence of discretized momentum equations for the x-and y-directions are as follows: and, , +1 2 ⁄ +1 The parameters and are The coefficients for the upwind advective terms are: The total water depths at the cell centers are evaluated by Furthermore, the total water depths at discharge points (i.e., locations of volume-flux components) are calculated as , +1 2 ⁄ +1 2 The total water depths of discharge points at = + 1 2 ⁄ are not calculated in the discretized mass equations; thus, they are evaluated from both time and space averages: , +1/2 = 0.25 • � , In a similar manner, the volume-flux components , +1/2 and +1 2 ⁄ , are We note that extremely large bottom frictions will occur when the total water depth is close to zero. Thus, the minimum total water depth threshold is required in simulations, which is 10 −5 m in this study. Following a similar procedure in calculating the upwinddiscretized advection terms, the total water depth threshold is also adopted.

Grid Nesting in Time and Space
COMCOT-SURGE adopts the two-way grid-nesting scheme to increase the resolutions in time and space for calculating storm surges. When adopting a smaller grid size, a corresponding time step following the CFL (Courant-Friedrichs-Lewy) condition is required to maintain the numerical stability in time and space: in which is the Courant Number, ℎ is the maximum still water depth (unit: m), and ∆ is the diagonal distance of two neighboring cell grids (∆ = √2∆ if in a square grid; unit: m). As discussed in Liu et al. (1998) [26], the Courant Number in the depthintegrated model using the leap-frog finite-difference method in the linear regime is approximately 0.70. For more practical wave modelling, the Courant Number is suggested as 0.65 for the linear model and 0.35 for the nonlinear model (see Wang and Power, 2011 [36]). It is noted that COMCOT-SURGE always needs to satisfy this condition when using the explicit leap-frog finite-difference method. Figure 1 shows the example of the grid-nesting function for the spatial domain in the grid-size ratio of 1:3 between a coarse and fine grid at the upper left and lower right corners. As shown in Figure 1, the cells of the inner grid (i.e., the fine grid) fully occupy the cells of the outer grid (i.e., coarse grid). The illustration here only explores the upper left and lower right corners between the coarse and fine grids, but the rest of the overlapped areas for the nested-grid domain are in the same manner. For waves propagating using a nested-grid domain, using a grid-size ratio from 3 to 5 to have a smooth transition across the connected boundaries between multiple grids is recommended [36]. It is also noted here that adopting a larger grid-size ratio needs to shoulder a smaller time step in the inner grid, restricted by the CFL condition. The grid-nesting function adopted in this study allows the outer and inner grids to use different time steps, but only in the time-step ratio of 2. Figure 2 illustrates the procedure of calculating the mass and momentum equations in the outer and inner grids having different time steps. As the volume-flux components have been computed at = Δ , the free surface elevations in the outer grid are evaluated by the mass equation (see Step 1 of Figure 2). After interpolating the volume-flux components from the outer grid to the inner grid (see Step 2 of Figure 2), the free surface elevations of the finer grid at = ( + 1/4)Δ are solved by the mass equation (see Step 3 of Figure 2), and the volumeflux components at = ( + 1/2)Δ are evaluated by the momentum equations (see Step 4 of Figure 2). To solve the free surface elevations of the finer grid at = ( + 3/4)Δ , the volume-flux components along the connected boundaries at = ( + 1/2)Δ are required. To get the volume-flux components along the connected edges at = ( + 1/2)Δ , they are interpolated in space and averaged in time (see Step 5 of Figure 2) from the volume-flux components of the outer grid at = Δ and = ( + 1)Δ . It is noted here that the volume-flux components of the outer grid at = ( + 1)Δ have not been addressed at that time; hence, they are "predicted" by solving the momentum equations with the information of the free surface elevations at = ( + 1/2)Δ . By having the volume-flux boundary conditions at = ( + 1/2)Δ , the free surface elevations of the inner grid at = ( + 3/4)Δ can be smoothly solved (see Step 6 of Figure 2). Subsequently, the volume-flux components of the inner grid at = ( + 1)Δ are evaluated by the discretized momentum equations (see Step 7 of Figure  2). If the two-way grid nesting is activated, the time-averaged free surface elevations of the finer grid at = ( + 1/4)Δ and = ( + 3/4)Δ are extrapolated back to the outer grid (see Step 8 of Figure 2). At the final step (Step 9 of Figure 2), the free surface elevations of the outer grid at = ( + 1)Δ are solved.

Moving Boundary Scheme
The moving boundary scheme adopted in this study allows the model to not only trace the moving shoreline with the free surface elevations, but also calculate the volumeflux components with the nonlinear shallow water equation model. Figure 3 shows the one-dimensional illustration for the adopted moving boundary scheme. As shown in Figure 3a, the shoreline stops between the grid cells i and i+1 because the free surface elevation at the grid cell i is not greater than the land elevation at i+1/2; thus, the volume-flux component at i+1/2 will not be calculated. As shown in Figure 3b, for another scenario that shows the free surface elevation at the cell i is greater than the land elevation at i+1/2, the flood depth exists and the volume-flux component at i+1/2 is calculated. Afterward, the free surface elevation at i+1 will be evaluated as the non-zero volume-flux component at i+1/2 exists. In this particular illustration, the flood depth, , is 0.5 × . We note here that the time-label of this one-dimensional illustration is ignored.  The moving boundary scheme focuses on dealing with the calculations of the volume-flux component with the movement of the shoreline. The nonlinear shallow water equations are adopted for calculating these physical properties for volume-flux components and : indicates the inland flood depth (unit: m), which replaces the total water depth, , when the flow crosses from dry cells to wet cells. If the volume-flux components exist between two wet cells, the momentum calculations will return to the standard procedures.

OpenMP Parallel Computing
COMCOT-SURGE supports using the OpenMP (Open Multi-Processing) parallelcomputing technique in a workstation, cluster, or personal computer, which has been used in iCOMCOT (i.e., cloud computing platform of COMCOT; see Lin et al., 2015 [32]). By using OpenMP, the paralleled version of COMCOT is about ten times faster than the serial version [32]. Thus, this OpenMP computing technique is extended from COMCOT to COMCOT-SURGE for calculating storm surges. The OpenMP parallel-computing technique is applied in the do-loops when solving the mass, momentum equations, and forcing terms (i.e., sea-level pressure gradient terms, wind shear stress terms, Coriolis terms, and horizontal eddy diffusion terms) during storm surge calculations.

Introduction
The experiment of the solitary wave runup on a circular island, designed to study the unexpected tsunami impacts on the lee side of a small island [27,37], provides a good data set to validate the model. The experimental data for the solitary wave runup on a circular island has been widely used to validate hydrodynamical models such as a shallow-water equation model [38] and Boussinesq-type model [39]. Thus, this study adopts the benchmark problem to examine: (1) the evolution of the free surface elevations near the coastline during the runup and rundown; (2) the model performance and numerical stability when tracing the moving shoreline with higher wave nonlinearity. Figure 4 shows the computational domain from top and side views. As illustrated in Figure 4a, the center of a circular island is located at x = 15 m and y = 13 m in a wave basin with the dimensions of 30 m in width and 25 m in length. The slope from the base of the island to the top is 1:4 (see Figure 4b). The top and base radius are 1.1 m and 3.6 m, respectively (see Figure 4b). The still water depth is 0.32 m, and the height from the still water surface to the island's top is 0.305 m. The four wave gauges recording the free surface elevations are available for the model validation, and the locations of the wave gauges can be found in Table 1. The incident solitary wave, generated by the wavemaker, propagates along the +y direction from the bottom boundary of the wave basin (the incident wave direction is indicated in Figure 4a). The formula of the incident solitary wave [40] is

Computational Setting
and, where is the incident solitary wave height (unit: m), is the effective wavenumber (unit: 1/m), and is the long-wave celerity (unit: m/s).  The two-layer nested computational domains (Grids 01 and 02) are adopted for the solitary wave runup on the circular island. As shown in Figure [36]. Grid 02 will accept the volume-flux components of Grid 01 as its boundary conditions, and it is noted that the two-way grid-nesting function is activated in time and space between Grids 01 and 02. The water density used in this benchmark simulation is 1000.0 kg/m 3 (i.e., the reference density of pure water). Briggs et al. (1995) [37] presented the experiments for the three different wave nonlinearities (A/h = 0.045, 0.091, and 0.181), but this study will only showcase the medium case. The reasons are: (1) the weakest wave-nonlinearity case (i.e., A/h = 0.045) may not be able to highlight the nonlinear effect with runup and rundown; (2) the largest wave-nonlinearity case (i.e., A/h = 0.181) has been found to have significant wave breaking phenomena around the circular island, which is not the main point of this study. Thus, this paper decides to shed light on the medium case (i.e., A/h = 0.091).

Computed Free Surface Elevations
The computed free surface elevations on the frontal and lee sides of the circular island are presented in Figures 5 and 6, respectively. At t = 8.5 and 9.0 s, the amplitude of the solitary wave becomes higher by the shoaling effect (see Figure 5a,b). At t = 9.5 and 10.0 sec, the amplified solitary wave inundates the frontal side of the circular island and generates a significant runup (see Figure 5c,d). In addition, the trapped waves propagate along the coastline, and the cylindrical wave patterns occur in front of the circular island after the wave runup generates on the frontal side of the island (see Figure 5). From t = 11.0, 11.5, 12.5, 13.5, and 14.0 s, the trapped waves propagate along the coastline (see Figure 6a-c), collide again (see Figure 6d), and generate another backward runup behind the island (see Figure 6e). It is noted that the inundated area on the frontal side is wider than those on the lee side. At t = 14.5, 15.0, and 16.0, the trapped waves propagate along the coastline again (see Figure 6f-

Time History of Free Surface Elevations
The measurement for free surface elevations at wave gauges on the frontal, lateral, and lee sides of the circular island (G6, G9, G16, and G22) from the Briggs et al. (1995)'s laboratory experiment are used to validate the model results. The measured water-level data in the time interval of 0.04 sec can be downloaded from the NOAA Center for Tsunami Research. (Website: https://nctr.pmel.noaa.gov/benchmark/Laboratory/Labora-tory_ConicalIsland/index.html, accessed on 10 January 2022). Figure 7 presents the comparisons between model results and measurements. In general, the computed free surface elevations agree well with the measured water levels in terms of the wave heights, arrival times, and wave shapes for the leading waves at G6, G9, G16, and G22 (see Figure 7). The correlation coefficients between the model results and measurements at G6, G9, G16, and G22 are 93.23%, 92.97%, 95.10%, and 92.87%, respectively. After the leading waves, the wave depressions are predicted less in the numerical model than in the measurements, which is consistent with the findings in the other depth-integrated model results [27,39]; this phenomenon has been discussed by Lynett et al. (2002) [39]. In addition, Titov and Synolakis (1998) [38] pointed out that wave breaking occurs on the lee side of the circular island; however, wave breaking occurs only on the lee side of the circular island and seems not to seriously affect the model prediction from the gauge comparison at G22 (see Figure 7).

Figure 7.
Comparisons between the computed free surface elevations (red lines) and measured water levels (black dots) (A/h = 0.091) at wave gauges G6, G9, G16, and G22, respectively. The wave gauge locations can be found in Table 1.

Runup Height and Inundation Area
The runup heights between the numerical results and measurements are projected onto the coordinates of the circular island in Figure 8. Here, the runup height is calculated as the maximum land elevation that waves arrive at from the original shoreline (z = 0 in the model simulations). As shown in Figure 8, the model results agree well with the measured runup heights on the frontal side of the island. However, on the lee side of the circular island where the wave breaking is found, the maximum runup heights predicted by the model are slightly lower than the measurements by about 0.5-1.0% (see Figure 8). Despite this, the leading order phenomena (i.e., significant runups on the lee side of the circular island) are well computed by the model. Thus, this implies that the maximum runup and inundation ranges are mainly contributed by the horizontal velocities rather than the vertical velocity changes due to the wave breaking in this particular runup case [38]. Hence, the model used here, without considering wave breaking, perfectly matches the measurements.

Introduction of 2013 Typhoon Haiyan
Typhoon Haiyan, also known as its local name Yolanda, was the strongest storm in 2013; it struck the Philippines with catastrophic storm surges, winds, and waves, and caused casualties of more than 6,300 [41]. The event of Typhoon Haiyan has three unique features: (1) a record-breaking wind speed of more than 310 km/hr [42]; (2) a fast forward motion of the storm at up to 41.0 km/hr [43]; (3) the notable induced storm surges and floods found in Leyte Gulf and San Pedro Bay [44][45][46][47]. Thus, these unique features make the 2013 Typhoon Haiyan a good case study for highlighting the model performance of predicting coastal storm surges and inundation areas.

Computational Setting
The three-layer nested-grid computational domains are adopted in this study to perform the storm surge simulation for the 2013 Typhoon Haiyan (see Figure 9). Table 2 tabulates each nested grid's computational domains, grid sizes, and corresponding time steps. By using the grid-nesting function, these three layers can simulate storm surge motions from offshore, nearshore, and coastal regions with appropriate grid sizes. The gridsize ratio between the outer grid and the inner grid is 3. In addition, we note here that all the settings of computational grids satisfy the CFL condition. The grid numbers of Domains D01, D02, and D03 are 69,276, 42,799, and 19,936, respectively; the total grid number is 132,011. The bathymetry data used in this study is from GEBCO (The General Bathymetric Chart of the Oceans) [48], which has a resolution of 15 arc-seconds in the latest GEBCO 2021 grid. The input 10-m winds and sea-level pressure fields are from the 1980 Holland Wind Model (Holland, 1980) [49] with the best-track data of JMA (Japan Meteorological Agency). The methodology for generating storm winds and surface air pressure fields, as well as the validation of these input meteorological fields with observations, is elaborated in the predecessor of this study [28].  Table 3). The switches for numerical calculations are modified accurately to simulate storm surges from offshore, nearshore, and coastal regions. This study assumes a storm surge propagating from the Western Pacific Ocean is under the linear regime; thus, the linear momentum equations are adopted to solve the offshore-scale storm surge motion. Afterward, when the storm surge propagates to nearshore and coastal regions, the nonlinear effect becomes more important; hence, the nonlinear momentum equations are conducted. Moreover, the inundation areas and flood depths shall be investigated near the Tacloban DZR (Daniel Z. Romualdez) Airport; thus, the moving boundary scheme (i.e., moving scheme) is turned on to trace shoreline and inland floods. In summary, Table 4 tabulates the switches for advection term, forcing terms (sea-level pressures, wind shear stresses, Coriolis force, bottom frictions, horizontal eddy diffusions), and moving boundary scheme for Domains D01, D02, and D03. In addition, some coefficients are required during the computations: the horizontal eddy diffusion coefficient is 100.0 m 2 /s [50]; the water density is the reference density of seawater (= 1025 kg/m 3 ) [51]; the Manning's coefficient is 0.025 s/m 1/3 [12]. Table 4. Switches for advection term, forcing terms, and moving boundary scheme. O indicates the switch is turned on; X indicates the switch is turned off.

Horizontal Eddy Diffusion Term
Moving Boundary Scheme

Storm Surges and Storm-Induced Currents
The Haiyan-induced storm surges propagate from the Western Pacific Ocean, to Leyte Gulf, and then to San Pedro Bay in the computational domains of our interests. Figure 10 shows the snapshots of computed storm surges in Domain D01. When the Super Typhoon Haiyan generates offshore winds, the negative storm surges occur accordingly in Leyte Gulf (see Figure 10a). When Typhoon Haiyan makes landfall, higher storm surges are generated near the coastline of Leyte Island and San Pedro Bay (see Figure 10b,c). With southeastern storm winds in San Pedro Bay, storm surges of more than 7 m occur (see Figure 10d), which causes dramatic impacts and damages to coastal communities in Tacloban (see, e.g., Mori et al., 2014 [52]; Soria et al., 2016 [46]). Figure 11 shows coastal storm surges and floods near Cancabato Bay and Tacloban DZR Airport, presented by the computed free surface elevations. The storm surges induced by Haiyan penetrates from Cancabato Bay to the coasts at 23:00 UTC on 7 November 2013 (see Figure 11a). Afterward, Haiyan-induced storm surges come from the east side of the Tacloban DZR Airport (see Figure 11b). After about 30 min to 1 h, larger coastal storm surges from San Pedro Bay enter Cancabato Bay and the east side of the Tacloban DZR Airport; thus, more inundation areas are found in the simulations (see Figure 11c,d). We note here that the elevation data used in the GEBCO 2021 grid may not be able to present accurate shapes and coastal infrastructures such as sea walls. In addition, the nearshore bathymetry is sensitive to hydrodynamical computations [53]. Hence, this particular case is used to showcase the model performance and ability of COMCOT-SURGE on the inundation calculation. Figure 12 presents the snapshots of the storm-induced current velocity fields. As shown in Figure 12a,b, an anti-clockwise vortex occurs inside Cancabato Bay. As found in the simulations, the maximum storm-induced current speed is about 3.5 m/s near the north tip of the Tacloban DZR Airport (see Figure 12b). Afterward, the storm-induced flows enter the Tacloban DZR Airport from both the east and north sides at 00:00 UTC on 8 November 2013 (Figure 12c). The storm-induced currents change to southward directions at 00:00 UTC on 8 November 2013, while the larger inundation areas are found in the simulation results (Figure 12d).

Maximum Storm Surges and Flood Depths
After exploring the storm surges and storm-induced currents around Cancabato Bay and the Tacloban Airport, this subsection further investigates the maximum storm surges and flood depths, which are essential to coastal communities. Figure 13a shows the maximum storm surges around the coasts of the Leyte Island. As shown in Figure 13a, the maximum computed surge heights are amplified from 3 m to 8 m in San Pedro Bay and the coasts of Leyte Island. The catastrophic storm surges contribute this phenomenon due to Typhoon Haiyan's southeastern winds propagating from Leyte Gulf to San Pedro Bay, as discussed in Figure 10. In addition, more detailed discussions about the storm surges related to storm winds can be found in Tsai et al. (2020) [28]. Furthermore, the maximum storm surges amplified from Leyte Gulf to San Pedro Bay can be also found in the discussion of Mori et al. (2014) [52], Kim et al. (2015) [12], and Soria et al. (2016) [46]. Figure 13b presents the maximum computed inland flood depths, which are the difference between the maximum computed free surface elevation and land elevation. As shown in Figure  13b, the maximum flood depth in the Tacloban DZR Airport shows the largest flood depths of about 6 m. Additionally, the coastal regions around Cancabato Bay suffered from significant floods of about 5 m. Although the land elevation from the GEBCO 2021 grid has not considered detailed infrastructures such as sea walls, the simulation results still expose dramatic flood depths around Cancabato Bay and the Tacloban DZR Airport corresponding to the field survey [46][47][48][49]. Basically, the storm-induced inundation areas agree with the measured extents by Tajima et al. (2014) [44] (see Figure 13b). Table 5 tabulates the flood depths between the measurements and model predictions of COMCOT-SURGE. Our storm surge model shows a good match with the flood depth at P2 but lower predictions at P1 and P2. On the one hand, we need to mention that the detailed digital elevation model (DEM) is not involved in our model; thus, the accuracy of inundation prediction can be improved once considering DEM or more accurate bathymetry [53]. On the other hand, the measured data may have some uncertainties, which may affect the validation of the inundation areas and flood depths. For example, P1, P2, and P3 show the reliabilities C, B, and A in the field survey, respectively (A: clear mark with small error; B: unclear mark but small error; C: unclear mark with large error [44]). Hence, the data (i.e., DEM, bathymetry, and field survey) are essential when discussing the storm-induced inundations.   Table  5). The yellow rectangles with the dashed black line mark inundation areas of Tajima et al. (2014) [44].

Time Series of Storm Surges
The time series of computed storm surges at specified numerical gauge stations (see Table 3) are illustrated in Figure 14. As shown in Figure 14, the negative storm surges occur from 18:00 UTC on 7 November 2013, and the water levels dramatically increase after 22:00 UTC on 7 November 2013 (see Basey, Tacloban, Palo, Tanauan, and Dulag). The maximum storm surge height of about 7.5 m is found at Tacloban (see Figure 14). This phenomenon is attributed to (1) the strong southeastern winds from Leyte Gulf to San Pedro Bay and (2) the changes of the wind directions due to the typhoon landfall. As discussed by Soria et al. (2016) [46], the water level from the trough to the crest within 1 to 2 h corresponds to the tsunami-like waves called by local residents, which is also discussed in Tsai et al. (2020) [28]. It is noted here that the maximum storm surges in the time-history data decay from north to south (see the Basey, Palo, Tanauan, and Dulag stations of Figure  14). In addition, after the largest amplitude, the following multiple crest-to-trough heights also diminish from north to south (see the Basey, Palo, and Tanauan stations of Figure 14).  Table 3.

Numerical Experiments
Section 5 has explored the 2013 Haiyan's storm surges and storm-induced flood/inundation areas using the three-layer nested-grid domain with the moving boundary scheme. Thus, this section will further discuss the effects of the nonlinear advection term and the fixed/moving shoreline, which are essential in simulating coastal storm surges, disaster assessments, and storm surge forecasting. In addition, the model efficiency boosted by OpenMP will be also investigated in this section.

Linear/Nolinear Equations with a Fixed or Moving Shoreline
When storm surges propagate to shallow waters, the surge amplitude increases by the wind shear stresses, wave radiation stresses, or bathymetry effect. Moreover, storm surges will penetrate from offshore to nearshore and inundate coastal communities. At these stages, the nonlinear effect may play an important role in simulating coastal storm surges; however, it is usually ignored in storm surge simulations or forecasting [10]. Additionally, the moving boundary (i.e., moving shoreline) with increased nonlinearity becomes challenging in storm surge simulations. Thus, this subsection will explore the nonlinear effect by conducting numerical experiments with the fixed and moving shorelines. Figure 15 explores the maximum storm surges of the nonlinear equation model with the moving/fixed shorelines and linear equation model with the fixed shoreline. The inland storm surges (i.e., free surface elevations or storm-induced floods) in Figure 15a are masked to compare each simulation. As shown in Figure 15, the maximum storm surges using the moving shoreline inside Cancabato Bay are about 0.5 lower than the fixed shoreline case. This corresponds to the arguments in Kowalik and Murty (1993) [54], that the water-level predictions of fixed shoreline are higher than the moving shoreline in a numerical model. However, the difference between the nonlinear and linear equation models under the fixed shoreline is relatively unclear, indicating a 0.1-0.2 m difference inside Cancabato Bay (see Figure 15b,c).  Figure 16 further investigates the storm-induced current fields between the nonlinear and linear equation model using the fixed shoreline. At 00:00 UTC on 8 November 2013, the nonlinear model shows an eddy inside Cancabato Bay, with its center located on the west coast of the Tacloban Airport (see Figure 16a). This eddy is also found in the model results using the nonlinear equation model with the moving boundary scheme (see Figure  12c). However, at 00:00 UTC on 8 November 2013, the simulation of the linear model with the fixed shoreline only shows anti-clockwise currents inside Cancabato Bay and no eddies are occurring (see Figure 16c). At 00:30 UTC on 8 November 2013, the storm-induced currents propagate along the east coasts of the Tacloban Airport in both nonlinear and linear models with the fixed shoreline (see Figure 16b,d). However, the results using the moving boundary scheme show the different flow patterns. Since the computed storm surges have inundated the Tacloban Airport, the flows pass over inland regions southeastward (see Figure 12d).

Parallel-Computing Efficiency
This section explores the model efficiency in the workstation with the CPU of AMD Ryzen 9 3900X (12 cores; in other words, 24 threads) and with a RAM of 64 GB under the operating system of the Linux-based CentOS 8. The clock time of each computation is calculated by the elapsed time between the first and the last output files. Here, we note again that the OpenMP algorithm is applied in the do-loop calculation when calculating the discretized mass and momentum equations and the forcing terms. But the output and input procedure are not accelerated in the model. Figure 17 shows that the clock time corresponds to the usage of thread numbers. As shown in Figure 17, the efficiency is dramatically enhanced when using more CPUs from 2 to 12 threads, but the efficiency stops being boosted after 12 threads (i.e., 6 CPUs). This implies that the OpenMP technique helps to enhance the storm surge calculation from the serial version to parallel version. However, after 12 threads in this particular case, the threads are overused. As shown in the running test in the workstation, the model efficiency follows the exponential curve (see Figure 17). Table 6

Conclusions
This study has developed a numerical tool with the two-way grid-nesting function in time and space, the moving boundary scheme in tracing the shoreline, and calculating coastal storm surges and inundation areas without any numerical filter adopted, which has been nicknamed COMCOT-SURGE. By validating the model performance and numerical stability against the solitary wave runup on a circular island, the time series of the free surface elevations and the runup heights perfectly matched the measurements. In addition, the free surface evolution of the incident solitary wave has been explored in the three-dimensional presentation with the runup, rundown, and trapped wave propagation around the circular island. After the benchmark problem validation, the extreme storm surge event of the 2013 Super Typhoon Haiyan was selected to showcase the storm surge computations of the developed model in this study. The three-layer nested-grid domains have been adopted to perform the storm surge simulations from the offshore, nearshore, to coastal regions. The extreme storm surges of about 8 m with the storm surge-induced inundation were explored in the Tacloban Airport and Cancabato Bay. Additionally, an anti-clockwise eddy in the storm-induced current fields was found inside Cancabato Bay.
The numerical experiments showed that the linear/nonlinear momentum equation models with the fixed shoreline predict higher storm surge amplitudes than the nonlinear momentum equation model with the moving shoreline around Cancabato Bay. However, the linear momentum equation model cannot resolve the eddy generated in Cancabato Bay. Furthermore, the model parallel efficiency has increased dramatically in using multiple threads as compared to a single thread. However, the boost of the model efficiency followed an exponential curve, implying the overuse of threads may not benefit the computations. By presenting a finite-difference-method numerical storm surge model with the two-way grid-nesting function and moving boundary scheme solving the nonlinear momentum equations, this study hopes to provide a convenient, useful, and flexible numerical tool for future research in evaluating storm surges in operational forecasting or disaster assessments.

Future Work
The future work is suggested as two parts: (1) open-source model to the research community and (2) model coupling with the spectral wave model. First, although the finite-difference model has been developed for a few decades, not many models have been shared or opened to the research community as open-source codes. Thus, an official website or GitHub/Gitlab webpage for COMCOT-SURGE will be developed for opening the source codes. Second, the wave-enhanced radiation stresses may play important roles in amplifying coastal storm surges [6][7][8]13]; thus, the coupling between COMCOT-SURGE and a spectral wave model such as SWAN (Simulating WAves Nearshore) shall be explored to improve prediction accuracy.  https://nctr.pmel.noaa.gov/benchmark/Laboratory/Laboratory_ConicalIsland/index.html (accessed on 10 January 2022). The GEBCO 2021 grid can be found at the website: https://www.gebco.net/data_and_products/gridded_bathymetry_data/gebco_2021/ (accessed on 10 January 2022).