Next Article in Journal
A Study on Darcy versus Forchheimer Models for Flow through Heterogeneous Landfills including Macropores
Next Article in Special Issue
Numerical Investigation of Sediment Flushing and Morphological Changes in Tamsui River Estuary through Monsoons and Typhoons
Previous Article in Journal
Runoff Probability Prediction Model Based on Natural Gradient Boosting with Tree-Structured Parzen Estimator Optimization
Previous Article in Special Issue
Solution of Shallow-Water Equations by a Layer-Integrated Hydrostatic Least-Squares Finite-Element Method
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Graduate Institute of Hydrological and Oceanic Sciences, National Central University, Chungli 320317, Taiwan
2
Institute of Physics, Academia Sinica, Taipei 115201, Taiwan
3
Research Center for Environmental Changes, Academia Sinica, Taipei 115201, Taiwan
*
Author to whom correspondence should be addressed.
Water 2022, 14(4), 547; https://doi.org/10.3390/w14040547
Submission received: 31 December 2021 / Revised: 27 January 2022 / Accepted: 8 February 2022 / Published: 12 February 2022
(This article belongs to the Special Issue Hydrodynamics in Ocean Environment: Experiment and Simulation)

Abstract

:
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.

1. 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 grid-nesting 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 grid-nesting 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/regional-scale 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 storm-induced 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 [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. Section 5 will explore the nonlinear advection term with the fixed/moving shoreline in simulating coastal storm surges of the 2013 Typhoon Haiyan. Finally, Section 6 will explore the nonlinear advection term with the fixed/moving shoreline in simulating coastal storm surges of the 2013 Typhoon Haiyan, we will conclude this study and illustrate some work for the future.

2. 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.

2.1. Governing Equation of the Storm Surge Model

The mass and momentum equations of COMCOT-SURGE are presented as follows:
η t + P x + Q y = 0 ,
P t + x P 2 H + y P Q H f Q = g H η x H ρ w P a x + τ s x ρ w τ b x ρ w + A h 2 P x 2 + 2 P y 2 ,
Q t + x P Q H + y Q 2 H + f P = g H η y H ρ w P a y + τ s y ρ w τ b y ρ w + A h 2 Q x 2 + 2 Q y 2 ,
where η is the free surface elevation (unit: m), t is the time (unit: s), ( x , y ) are the spatial notations of the x- and y-directions (unit: m), ( P , Q ) are the volume-flux components (i.e., depth-integrated volume fluxes per unit length) of the x- and y-directions (unit: m2/s), H is the total water depth ( H = h + η ; unit: m), h is the still water depth (unit: m), g is the gravitational acceleration (=9.81 m/s2), f is the Coriolis parameter ( f = 2 ω s i n φ ; unit: 1/s), ω is the Earth angular velocity (=7.2921 × 10−5 rad/s), P a is the sea-level air pressure (unit: N/m2), ( τ s x , τ s y ) are the wind shear stresses (unit: N/m2), ( τ b x , τ b y ) are the bottom frictional shear stresses (unit: N/m2), ρ w is the water density (unit: kg/m3), and A h is the horizontal eddy diffusion coefficient (unit: m2/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:
P t f Q = g H η x H ρ w P a x + τ s x ρ w τ b x ρ w + A h 2 P x 2 + 2 P y 2 ,
Q t + f P = g H η y H ρ w P a y + τ s y ρ w τ b y ρ w + A h 2 Q x 2 + 2 Q y 2 .
The Manning’s formula, from the conception of the open-channel flow, is used to model the bottom friction:
τ b x = ρ w g N 2 H 7 / 3 P P 2 + Q 2 ,
τ b y = ρ w g N 2 H 7 / 3 Q P 2 + Q 2 ,
where N is the Manning’s roughness coefficient (unit: s/m1/3).
The quadric law is used to model the wind shear stresses on the water surface:
τ s x = ρ a C d u 10 u 10 x ,
τ s y = ρ a C d u 10 u 10 y ,
in which C d is the wind-drag coefficient between the water surface and the air, u 10 is the 10-m wind speed (unit: m/s), ( u 10 x ,   u 10 y ) 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]:
C d = 1.2875 × 10 3 ,       u 10 < 7.5 , 0.8 + 0.065 × u 10 × 10 3 ,       7.5 u 10 < 25 , 0.8 + 0.065 × 25 × 10 3 ,       u 10 25
Here we note that C d is a dimensionless coefficient.

2.2. 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:
η i , j n + 1 / 2 = η i , j n 1 / 2 + Δ t Δ x P i + 1 / 2 , j n P i 1 / 2 , j n + Δ t Δ y Q i , j + 1 / 2 n Q i , j 1 / 2 n ,
where Δ t is the time step (unit: s) and ( Δ x , Δ y ) 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:
P i + 1 / 2 , j n + 1 = 1 v x Δ t 1 + v x Δ t P i + 1 / 2 , j n Δ t 1 + v x Δ t 1 Δ x λ 11 P i + 3 / 2 , j n 2 H i + 3 / 2 , j n + λ 12 P i + 1 / 2 , j n 2 H i + 1 / 2 , j n + λ 13 P i 1 / 2 , j n 2 H i 1 / 2 , j n Δ t 1 + v x Δ t 1 Δ y λ 21 P i + 1 / 2 , j + 1 n · Q i + 1 / 2 , j + 1 n H i + 1 / 2 , j + 1 n + λ 22 P i + 1 / 2 , j n · Q i + 1 / 2 , j n H i + 1 / 2 , j n + λ 23 P i + 1 / 2 , j 1 n · Q i + 1 / 2 , j 1 n H i + 1 / 2 , j 1 n Δ t 1 + v x Ɗ t g H i + 1 / 2 , j n + 1 / 2 η i + 1 , j n + 1 / 2 η i , j n + 1 / 2 Δ x + 1 ρ w H i + 1 / 2 , j n + 1 / 2 P a i + 1 , j n + 1 P a i , j n + 1 Δ x + Δ t 1 + v x Δ t ρ a ρ w C d i + 1 / 2 , j n + 1 U 10 i + 1 / 2 , j n + 1 U 10 x i + 1 / 2 , j n + 1 + f i + 1 / 2 , j · Q i + 1 / 2 , j n + Δ t 1 + v x Δ t A h P i + 3 / 2 , j n 2 P i + 1 / 2 , j n + P i 1 / 2 j n Δ x 2 + P i + 1 / 2 , j + 1 n 2 P i + 1 / 2 , j n + P i + 1 / 2 , j 1 n Δ y 2 ,
and,
Q i , j + 1 / 2 n + 1 = 1 v y Δ t 1 + v y Δ t Q i , j + 1 / 2 n Δ t 1 + v y Δ t 1 Δ x λ 31 P i + 1 , j + 1 / 2 n · Q i + 1 , j + 1 / 2 n H i + 1 , j + 1 / 2 n + λ 32 P i , j + 1 / 2 n · Q i , j + 1 / 2 n H i , j + 1 / 2 n + λ 33 P i 1 , j + 1 / 2 n · Q i 1 , j + 1 / 2 n H i 1 , j + 1 / 2 n Δ t 1 + v y Δ t 1 Δ y λ 41 Q i , j + 3 / 2 n 2 H i , j + 3 / 2 n + λ 42 Q i , j + 1 / 2 n 2 H i , j + 1 / 2 n + λ 43 Q i , j 1 / 2 n 2 H i , j 1 / 2 n Δ t 1 + v y Δ t g H i , j + 1 / 2 n + 1 / 2 η i , j + 1 n + 1 / 2 η i , j n + 1 / 2 Δ y + 1 ρ w H i , j + 1 / 2 n + 1 / 2 P a i , j + 1 n + 1 / 2 P a i , j n + 1 / 2 Δ y + Δ t 1 + v y Δ t ρ a ρ w C d i , j + 1 / 2 n + 1 U 10 i , j + 1 / 2 n + 1 U 10 y i , j + 1 / 2 n + 1 f i , j + 1 / 2 · P i , j + 1 / 2 n + Δ t 1 + v y Δ t A h Q i + 1 , j + 1 / 2 n 2 Q i , j + 1 / 2 n + Q i 1 , j + 1 / 2 n Δ x 2 + Q i , j + 3 / 2 n 2 Q i , j + 1 / 2 n + Q i , j 1 / 2 n Δ y 2 .
The parameters v x and v y are
v x = 1 2 g N 2 H i + 1 / 2 , j n 7 / 3 P i + 1 / 2 , j n 2 + Q i + 1 / 2 , j n 2 1 / 2 ,
v y = 1 2 g N 2 H i , j + 1 / 2 n 7 / 3 P i , j + 1 / 2 n 2 + Q i , j + 1 / 2 n 2 1 / 2 .
The coefficients for the upwind advective terms are:
λ 11 = 0 ,     λ 12 = 1 ,     λ 13 = 1 ,     i f   P i + 1 / 2 , j n 0 , λ 11 = 1 ,     λ 12 = 1 ,     λ 13 = 0 ,     i f   P i + 1 / 2 , j n < 0 ,
λ 21 = 0 ,     λ 22 = 1 ,     λ 23 = 1 ,     i f   Q i + 1 / 2 , j n 0 , λ 21 = 1 ,     λ 22 = 1 ,     λ 23 = 0 ,     i f   Q i + 1 / 2 , j n < 0 ,
λ 31 = 0 ,     λ 32 = 1 ,     λ 33 = 1 ,     i f   P i , j + 1 / 2 n 0 , λ 31 = 1 ,     λ 32 = 1 ,     λ 33 = 0 ,     i f   P i , j + 1 / 2 n < 0 ,
λ 41 = 0 ,     λ 42 = 1 ,     λ 43 = 1 ,     i f   Q i , j + 1 / 2 n 0 , λ 41 = 1 ,     λ 42 = 1 ,     λ 43 = 0 ,     i f   Q i , j + 1 / 2 n < 0 .
The total water depths at the cell centers are evaluated by
H i , j n + 1 / 2 = η i , j n + 1 / 2 + h i , j .
Furthermore, the total water depths at discharge points (i.e., locations of volume-flux components) are calculated as
H i + 1 / 2 , j n + 1 / 2 = 0.5 · H i , j n + 1 / 2 + H i + 1 , j n + 1 / 2 ,
H i , j + 1 / 2 n + 1 / 2 = 0.5 · H i , j n + 1 / 2 + H i , j + 1 n + 1 / 2 .
The total water depths of discharge points at t = n + 1 / 2 are not calculated in the discretized mass equations; thus, they are evaluated from both time and space averages:
H i + 1 / 2 , j n = 0.25 · H i , j n 1 / 2 + H i , j n + 1 / 2 + H i + 1 , j n 1 / 2 + H i + 1 , j n + 1 / 2 ,
H i , j + 1 / 2 n = 0.25 · H i , j n 1 / 2 + H i , j n + 1 / 2 + H i , j + 1 n 1 / 2 + H i , j + 1 n + 1 / 2 .
In a similar manner, the volume-flux components P i , j + 1 / 2 n and Q i + 1 / 2 , j n are
P i , j + 1 / 2 n = 0.25 · P i 1 / 2 , j n + P i 1 / 2 , j + 1 n + P i + 1 / 2 , j n + P i + 1 / 2 , j + 1 n ,
Q i + 1 / 2 , j n = 0.25 · Q i , j 1 / 2 n + Q i + 1 , j 1 / 2 n + Q i , j + 1 / 2 n + Q i + 1 , j + 1 / 2 n .
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 upwind-discretized advection terms, the total water depth threshold is also adopted.

2.3. 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:
Δ t C r Δ s 2 g h m a x ,
in which C r is the Courant Number, h m a x is the maximum still water depth (unit: m), and Δ s is the diagonal distance of two neighboring cell grids ( Δ s = 2 Δ x if in a square grid; unit: m). As discussed in (1998) [26], the Courant Number in the depth-integrated 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 t = n Δ t o u t , 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 t = n + 1 / 4 Δ t i n are solved by the mass equation (see Step 3 of Figure 2), and the volume-flux components at t = n + 1 / 2 Δ t i n are evaluated by the momentum equations (see Step 4 of Figure 2). To solve the free surface elevations of the finer grid at t = n + 3 / 4 Δ t i n , the volume-flux components along the connected boundaries at t = n + 1 / 2 Δ t i n are required. To get the volume-flux components along the connected edges at t = n + 1 / 2 Δ t i n , 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 t = n Δ t o u t and t = n + 1 Δ t o u t . It is noted here that the volume-flux components of the outer grid at t = n + 1 Δ t o u t 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 t = n + 1 / 2 Δ t o u t . By having the volume-flux boundary conditions at t = n + 1 / 2 Δ t i n , the free surface elevations of the inner grid at t = n + 3 / 4 Δ t i n can be smoothly solved (see Step 6 of Figure 2). Subsequently, the volume-flux components of the inner grid at t = n + 1 Δ t i n 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 t = n + 1 / 4 Δ t i n and t = n + 3 / 4 Δ t i n 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 t = n + 1 Δ t o u t are solved.

2.4. 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 volume-flux 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 H f 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, H f , is 0.5 × H i . 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 P and Q :
P i + 1 / 2 , j n + 1 = 1 v x Δ t 1 + v x Δ t P i + 1 / 2 , j n Δ t 1 + v x Δ t 1 Δ x λ 11 P i + 3 / 2 , j n 2 H i + 3 / 2 , j n + λ 12 P i + 1 / 2 , j n 2 H i + 1 / 2 , j n + λ 13 P i 1 / 2 , j n 2 H i 1 / 2 , j n Δ t 1 + v x Δ t 1 Δ y λ 21 P i + 1 / 2 , j + 1 n · Q i + 1 / 2 , j + 1 n H i + 1 / 2 , j + 1 n + λ 22 P i + 1 / 2 , j n · Q i + 1 / 2 , j n H i + 1 / 2 , j n + λ 23 P i + 1 / 2 , j 1 n · Q i + 1 / 2 , j 1 n H i + 1 / 2 , j 1 n Δ t 1 + v x Δ t g H f i + 1 / 2 , j n + 1 / 2 η i + 1 , j n + 1 / 2 η i , j n + 1 / 2 Δ x ,
and,
Q i , j + 1 / 2 n + 1 = 1 v y Δ t 1 + v y Δ t Q i , j + 1 / 2 n Δ t 1 + v y Δ t 1 Δ x λ 31 P i + 1 , j + 1 / 2 n · Q i + 1 , j + 1 / 2 n H i + 1 , j + 1 / 2 n + λ 32 P i , j + 1 / 2 n · Q i , j + 1 / 2 n H i , j + 1 / 2 n + λ 33 P i 1 , j + 1 / 2 n · Q i 1 , j + 1 / 2 n H i 1 , j + 1 / 2 n Δ t 1 + v y Δ t 1 Δ y λ 41 Q i , j + 3 / 2 n 2 H i , j + 3 / 2 n + λ 42 Q i , j + 1 / 2 n 2 H i , j + 1 / 2 n + λ 43 Q i , j 1 / 2 n 2 H i , j 1 / 2 n Δ t 1 + v y Δ t g H f i , j + 1 / 2 n + 1 / 2 η i , j + 1 n + 1 / 2 η i , j n + 1 / 2 Δ y .
H f indicates the inland flood depth (unit: m), which replaces the total water depth, H , 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.

2.5. OpenMP Parallel Computing

COMCOT-SURGE supports using the OpenMP (Open Multi-Processing) parallel-computing 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.

3. Model Validation—Solitary Wave Runup on a Circular Island

3.1. 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.

3.2. Computational Setting

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
η x , t = A s e c h 2 κ x c t ,
and,
κ = 3 A 4 h 3 ,
c = g h + A
where A is the incident solitary wave height (unit: m), κ is the effective wavenumber (unit: 1/m), and c 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 4a, the domain of Grid 01 is (0.0–30.0) m and (0.0–25.0) m, and the domain of Grid 02 is (10.0–20.0) and (9.0–17.0) m. The frictionless linear shallow water equations are used in Grid 01; the nonlinear shallow water equations with the bottom friction (Manning’s coefficient is 0.013 s/m1/3) are adopted in Grid 02 for the regions near the circular island. The grid size of Grid 01 is 0.1 m, and the grid size ratio between Grid 01 and Grid 02 is 3:1. The time step of Grid 01 is 0.01 s, and the time step ratio between Grid 01 and Grid 02 is 2:1. The boundary conditions for the left and right sides of Grid 01 are the wall-boundary condition, and the downstream side of Grid 01 is the radiation boundary condition [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/m3 (i.e., the reference density of pure water).
Ref. [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).

3.3. Computed Free Surface Elevations

The computed free surface elevations on the frontal and lee sides of the circular island are presented in Figure 5 and Figure 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–h).

3.4. 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 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/Laboratory_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 [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).

3.5. 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.

4. Case Study of Storm Surges—2013 Super Typhoon Haiyan

4.1. 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.

4.2. 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 grid-size 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].
Figure 9. Computational domains of the three-layer nested grids in Leyte Gulf and San Pedro Bay. (a) Domain D01 with Domain D02. The color shading indicates the water depths in m. The dashed black line indicates Domain D02. The dashed white line indicates the 50-m water-depth contour. (b) Domain D02 with Domain D03. The color shading indicates the water depths in m. The solid black line indicates Domain D03. The dashed blue line indicates the 10-m water-depth contour. The green circles show the numerical gauge locations (see Table 3).
Figure 9. Computational domains of the three-layer nested grids in Leyte Gulf and San Pedro Bay. (a) Domain D01 with Domain D02. The color shading indicates the water depths in m. The dashed black line indicates Domain D02. The dashed white line indicates the 50-m water-depth contour. (b) Domain D02 with Domain D03. The color shading indicates the water depths in m. The solid black line indicates Domain D03. The dashed blue line indicates the 10-m water-depth contour. The green circles show the numerical gauge locations (see Table 3).
Water 14 00547 g009
Table 2. Nested-grid computational domains, grid sizes, and time steps for calculating storm surges induced by the 2013 Super Typhoon Haiyan: D01-coarse grid, D02-medium grid, and D03-fine grid.
Table 2. Nested-grid computational domains, grid sizes, and time steps for calculating storm surges induced by the 2013 Super Typhoon Haiyan: D01-coarse grid, D02-medium grid, and D03-fine grid.
DomainLongitude
(Unit: °E)
Latitude
(Unit: °N)
Grid   Size   ( Δ x ,   Δ y )
(Unit: m)
Time   Step   ( Δ t )
(Unit: s)
D01124.9–126.010.4–11.4(437.24, 445.28)0.4
D02124.98–125.1510.85–11.30(145.66, 148.42)0.2
D03124.99–125.0411.18–11.26(48.52, 49.48)0.1
Table 3. Numerical gauge locations and their domain belongings.
Table 3. Numerical gauge locations and their domain belongings.
Station NameLongitude
(Unit: °E)
Latitude
(Unit: °N)
Domain
Basey125.052311.2723D02
Tacloban125.000411.2538D03
Palo125.011811.1561D02
Tanauan125.026411.1016D02
Dulag125.041310.9465D02
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 m2/s [50]; the water density is the reference density of seawater (=1025 kg/m3) [51]; the Manning’s coefficient is 0.025 s/m1/3 [12].

4.3. 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).

4.4. 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 [28]. Furthermore, the maximum storm surges amplified from Leyte Gulf to San Pedro Bay can be also found in the discussion of [12,46,52]. 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.

4.5. 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 [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 [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).

5. Numerical Experiments

Section 6 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.

5.1. 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).

5.2. 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 presents the clock time corresponding to the usage of 1, 2, 4, 8, 12, 16, 20, and 24 threads.

6. Conclusions and Future Work

6.1. 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.

6.2. 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.

Author Contributions

Y.-L.T. developed the storm surge model, performed the model validation/storm surge simulations, and wrote the manuscript. T.-R.W. organized this study, coordinated with all co-authors, and found the research funding. E.Y. and S.C.L. provided the conceptualization of this study and found the research funding. C.-Y.L. provided the knowledge of meteorological storm fields. All authors have read and agreed to the published version of the manuscript.

Funding

T.-R.W. would like to appreciate the support received from the Central Weather Bureau (Taiwan) under the Grant No. of MOTC-CWB-110-O-03, and E.Y. and S.C.L. want to appreciate the grants funded by UND TEIN*CC under the Grant No. of Asi@Connect-18-066.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The benchmark data for a solitary wave runup on a circular island case can be found at the NOAA Pacific Marine Environmental Laboratory (PMEL) website: 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).

Acknowledgments

The authors want to thank P.L.-F. Liu (National University of Singapore) for consulting on the moving boundary scheme used in COMCOT-SURGE, GEBCO Compilation Group (2021) for the GEBCO 2021 Grid (doi:10.5285/c6612cbe-50b3-0cff-e053-6c86abc09f8f), and Tajima et al., (2014)’s paper for providing valuable field survey data.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Fritz, H.M.; Blount, C.; Sokoloski, R.; Singleton, J.; Fuggle, A.; McAdoo, B.G.; Moore, A.; Grass, C.; Tate, B. Hurricane Katrina storm surge distribution and field observations on the Mississippi Barrier Islands. Estuarine Coast. Shelf Sci. 2007, 74, 12–20. [Google Scholar] [CrossRef]
  2. Emanuel, K. Increasing destructiveness of tropical cyclones over the past 30 years. Nature 2005, 436, 686–688. [Google Scholar] [CrossRef] [PubMed]
  3. Sun, J.; Oey, L.; Xu, F.-H.; Lin, Y.-C. Sea level rise, surface warming, and the weakened buffering ability of South China Sea to strong typhoons in recent decades. Sci. Rep. 2017, 7, 7418. [Google Scholar] [CrossRef]
  4. Zhang, W.-Z.; Shi, F.; Hong, H.-S.; Shang, S.-P.; Kirby, J.T. Tide-surge interaction intensified by the Taiwan Strait. J. Geophys. Res. Earth Surf. 2010, 115. [Google Scholar] [CrossRef] [Green Version]
  5. Tang, Y.M.; Sanderson, B.; Holland, G.; Grimshaw, R. A numerical study of storm surges and tides, with application to the North Queensland coast. J. Phys. Oceanogr. 1996, 26, 2700–2711. [Google Scholar] [CrossRef] [Green Version]
  6. Mastenbroek, C.; Burgers, G.; Janssen, P.A.E.M. The dynamical coupling of a wave model and a storm surge model through the atmospheric boundary layer. J. Phys. Oceanogr. 1993, 23, 1856–1866. [Google Scholar] [CrossRef] [Green Version]
  7. Bunya, S.; Dietrich, J.C.; Westerink, J.J.; Ebersole, B.A.; Smith, J.M.; Atkinson, J.H.; Jensen, R.; Resio, D.T.; Luettich, R.A.; Dawson, C.; et al. A high-resolution coupled riverine flow, tide, wind, wind wave, and storm surge model for southern Louisiana and Mississippi. Part I: Model development and validation. Mon. Weather Rev. 2010, 138, 345–377. [Google Scholar] [CrossRef] [Green Version]
  8. Dietrich, J.C.; Bunya, S.; Westerink, J.J.; Ebersole, B.A.; Smith, J.M.; Atkinson, J.H.; Jensen, R. A high-resolution coupled riverine flow, tide, wind, wind wave, and storm surge model for southern Louisiana and Mississippi. Part II: Synoptic description and analysis of Hurricanes Katrina and Rita. Mon. Weather Rev. 2010, 138, 378–404. [Google Scholar] [CrossRef] [Green Version]
  9. Kennedy, A.; Gravois, U.; Zachry, B.C.; Westerink, J.J.; Hope, M.E.; Dietrich, J.; Powell, M.; Cox, A.T.; Luettich, R.A.; Dean, R.G. Origin of the Hurricane Ike forerunner surge. Geophys. Res. Lett. 2011, 38. [Google Scholar] [CrossRef]
  10. Jelesnianski, C.P. SLOSH: Sea, Lake and Overland Surges from Hurricanes; U.S. Department of Commerce, National Oceanic and Atmospheric Administration, National Weather Service: Silver Spring, MD, USA, 1992; Volume 48.
  11. Zhang, K.; Li, Y.; Liu, H.; Rhome, J.; Forbes, C. Transition of the Coastal and Estuarine Storm Tide Model to an Operational Storm Surge Forecast Model: A Case Study of the Florida Coast. Weather Forecast. 2013, 28, 1019–1037. [Google Scholar] [CrossRef]
  12. Kim, S.; Mori, N.; Mase, H.; Yasuda, T. The role of sea surface drag in a coupled surge and wave model for Typhoon Haiyan 2013. Ocean Model. 2015, 96, 65–84. [Google Scholar] [CrossRef]
  13. Sheng, Y.P.; Alymov, V.; Paramygin, V.A. Simulation of storm surge, wave, currents, and inundation in the Outer Banks and Chesapeake Bay during Hurricane Isabel in 2003: The importance of waves. J. Geophys. Res. Earth Surf. 2010, 115. [Google Scholar] [CrossRef]
  14. Warner, J.C.; Armstrong, B.; He, R.; Zambon, J.B. Development of a coupled ocean–atmosphere–wave–sediment transport (COAWST) modeling system. Ocean Model. 2010, 35, 230–244. [Google Scholar] [CrossRef] [Green Version]
  15. Weisberg, R.H.; Zheng, L. Hurricane storm surge simulations comparing three-dimensional with two-dimensional formulations based on an Ivan-like storm over the Tampa Bay, Florida region. J. Geophys. Res. Ocean. 2008, 113. [Google Scholar] [CrossRef] [Green Version]
  16. Zhang, Y.J.; Ye, F.; Stanev, E.V.; Grashorn, S. Seamless cross-scale modeling with SCHISM. Ocean Model. 2016, 102, 64–81. [Google Scholar] [CrossRef] [Green Version]
  17. Dietrich, J.; Zijlema, M.; Westerink, J.; Holthuijsen, L.; Dawson, C.; Luettich, R.; Jensen, R.; Smith, J.; Stelling, G.; Stone, G. Modeling hurricane waves and storm surge using integrally-coupled, scalable computations. Coast. Eng. 2011, 58, 45–65. [Google Scholar] [CrossRef]
  18. Cheung, K.; Phadke, A.; Wei, Y.; Rojas, R.; Douyere, Y.-M.; Martino, C.; Houston, S.; Liu, P.; Lynett, P.; Dodd, N.; et al. Modeling of storm-induced coastal flooding for emergency management. Ocean Eng. 2003, 30, 1353–1386. [Google Scholar] [CrossRef]
  19. Flather, R.A. Existing operational oceanography. Coast. Eng. 2000, 41, 13–40. [Google Scholar] [CrossRef]
  20. Hasegawa, H.; Kohno, N.; Itoh, M. Development of Storm Surge Model in Japan Meteorological Agency. In Proceedings of the 2nd JCOMM Scientific and Technical Symposium, Key West, FL, USA, 8–13 November 2015. [Google Scholar]
  21. Yu, Y.-C.; Chen, H.; Shih, H.-J.; Chang, C.-H.; Hsiao, S.-C.; Chen, W.-B.; Chen, Y.-M.; Su, W.-R.; Lin, L.-Y. Assessing the potential highest storm tide hazard in Taiwan based on 40-year historical typhoon surge hindcasting. Atmosphere 2019, 10, 346. [Google Scholar] [CrossRef] [Green Version]
  22. Chen, W.-B.; Liu, W.-C. Assessment of storm surge inundation and potential hazard maps for the southern coast of Taiwan. Nat. Hazards 2016, 82, 591–616. [Google Scholar] [CrossRef]
  23. Li, N.; Yamazaki, Y.; Roeber, V.; Cheung, K.F.; Chock, G. Probabilistic mapping of storm-induced coastal inundation for climate change adaptation. Coast. Eng. 2018, 133, 126–141. [Google Scholar] [CrossRef]
  24. Medeiros, S.C.; Hagen, S.C. Review of wetting and drying algorithms for numerical tidal flow models. Int. J. Numer. Methods Fluids 2012, 71, 473–487. [Google Scholar] [CrossRef]
  25. Forbes, C.; Rhome, J.; Mattocks, C.; Taylor, A. Predicting the storm surge threat of Hurricane Sandy with the National Weather Service SLOSH Model. J. Mar. Sci. Eng. 2014, 2, 437–476. [Google Scholar] [CrossRef] [Green Version]
  26. Liu, P.L.-F.; Woo, S.-B.; Cho, Y.-S. Computer Programs for Tsunami Propagation and Inundation; Cornell University: Ithaca, NY, USA, 1998. [Google Scholar]
  27. Liu, P.L.F.; Cho, Y.-S.; Briggs, M.J.; Kanoglu, U.; Synolakis, C.E. Runup of solitary waves on a circular Island. J. Fluid Mech. 1995, 302, 259–285. [Google Scholar] [CrossRef]
  28. Tsai, Y.-L.; Wu, T.-R.; Lin, C.-Y.; Lin, S.C.; Yen, E.; Lin, C.-W. Discrepancies on storm surge predictions by parametric wind model and numerical weather prediction model in a semi-enclosed bay: Case study of typhoon Haiyan. Water 2020, 12, 3326. [Google Scholar] [CrossRef]
  29. Lin, Y.-H.; Fang, M.-C.; Hwung, H.-H. Transport reversal due to Typhoon Krosa in the Taiwan Strait. Open Ocean Eng. Journa 2010, 3, 143–157. [Google Scholar] [CrossRef]
  30. Lin, Y.-H.; Hwung, H.-H.; Fang, M.-C. The Numerical Simulation of Storm-Surge and Coastal Flooding in Western Taiwan: A Case Study of 2007 Typhoon SEPAT. J. Shipp. Ocean. Eng. 2011, 1. [Google Scholar] [CrossRef]
  31. Cho, Y.-S. Numerical Simulations of Tsunami Propagation and Run-up. Ph.D. Thesis, Cornell University, Ithaca, NY, USA, 1995. [Google Scholar]
  32. Lin, S.C.; Wu, T.-R.; Yen, E.; Chen, H.-Y.; Hsu, J.; Tsai, Y.-L.; Lee, C.-J.; Philip, L.-F.L. Development of a tsunami early warning system for the South China Sea. Ocean Eng. 2015, 100, 1–18. [Google Scholar] [CrossRef]
  33. Yen, E.; Lin, S.C.; Wu, T.-R.; Tsai, Y.-L.; Chung, M.-J. Knowledge-Building Approach for Tsunami Impact Analysis Aided by Citizen Science. Front. Earth Sci. 2020, 8, 315. [Google Scholar] [CrossRef]
  34. Wu, J. Wind-stress coefficients over sea surface from breeze to hurricane. J. Geophys. Res. Ocean. 1982, 87, 9704–9706. [Google Scholar] [CrossRef]
  35. WAMDI. The WAM model—A third generation ocean wave prediction model. J. Phys. Oceanogr. 1988, 18, 1775–1810. [Google Scholar] [CrossRef] [Green Version]
  36. Wang, X.; Power, W. COMCOT: A Tsunami Generation, Propagation and Run-Up Model; GNS Science: Lower Hutt, New Zealand, 2011. [Google Scholar]
  37. Briggs, M.J.; Synolakis, C.E.; Harkins, G.S.; Green, D.R. Laboratory experiments of tsunami runup on a circular island. Pure Appl. Geophys. 1995, 144, 569–593. [Google Scholar] [CrossRef]
  38. Titov, V.; Synolakis, C.E. Numerical modeling of tidal wave runup. J. Waterw. Port Coast. Ocean. Eng. 1998, 124, 157–171. [Google Scholar] [CrossRef]
  39. Lynett, P.; Wu, T.-R.; Liu, P. Modeling wave runup with depth-integrated equations. Coast. Eng. 2002, 46, 89–107. [Google Scholar] [CrossRef]
  40. Dean, R.G.; Dalrymple, R.A. Water Wave Mechanics for Engineers and Scientists; World Scientific Publishing Company: Singapore, 1991. [Google Scholar] [CrossRef]
  41. NDRRMC. Effects of Typhoon “YOLANDA” (HAIYAN); Technical Report; National Disaster Risk Reduction and Management Council: Quezon City, Philippines, 2014.
  42. Schiermeier, Q. Did Climate Change Cause Typhoon Haiyan? Nature 2013, 11. [Google Scholar] [CrossRef]
  43. Takagi, H.; Esteban, M.; Shibayama, T.; Mikami, T.; Matsumaru, R.; De Leon, M.; Thao, N.; Oyama, T.; Nakamura, R. Track analysis, simulation, and field survey of the 2013 Typhoon Haiyan storm surge. J. Flood Risk Manag. 2014, 10, 42–52. [Google Scholar] [CrossRef]
  44. Tajima, Y.; Yasuda, T.; Pacheco, B.M.; Cruz, E.C.; Kawasaki, K.; Nobuoka, H.; Miyamoto, M.; Asano, Y.; Arikawa, T.; Ortigas, N.M.; et al. Initial report of JSCE-PICE joint survey on the storm surge disaster caused by Typhoon Haiyan. Coast. Eng. J. 2014, 56, 1450006. [Google Scholar] [CrossRef]
  45. Mas, E.; Bricker, J.; Kure, S.; Adriano, B.; Yi, C.; Suppasri, A.; Koshimura, S. Field survey report and satellite image interpretation of the 2013 Super Typhoon Haiyan in the Philippines. Nat. Hazards Earth Syst. Sci. 2015, 15, 805–816. [Google Scholar] [CrossRef] [Green Version]
  46. Soria, J.L.A.; Switzer, A.D.; Villanoy, C.L.; Fritz, H.M.; Bilgera, P.H.T.; Cabrera, O.C.; Siringan, F.P.; Maria, Y.Y.-S.; Ramos, R.D.; Fernandez, I.Q. Repeat storm surge disasters of Typhoon Haiyan and its 1897 predecessor in the Philippines. Bull. Am. Meteorol. Soc. 2016, 97, 31–48. [Google Scholar] [CrossRef]
  47. Mikami, T.; Shibayama, T.; Takagi, H.; Matsumaru, R.; Esteban, M.; Thao, N.D.; Kumagaim, K. Storm surge heights and damage caused by the 2013 Typhoon Haiyan along the Leyte Gulf coast. Coast. Eng. J. 2016, 58, 1640005. [Google Scholar] [CrossRef]
  48. Weatherall, P.; Marks, K.M.; Jakobsson, M.; Schmitt, T.; Tani, S.; Arndt, J.E.; Rovere, M.; Chayes, D.; Ferrini, V.; Wigley, R. A new digital bathymetric model of the world’s oceans. Earth Space Sci. 2015, 2, 331–345. [Google Scholar] [CrossRef]
  49. Holland, G.J. An analytic model of the wind and pressure profiles in hurricanes. Mon. Weather Rev. 1980, 108, 1212–1218. [Google Scholar] [CrossRef]
  50. Zu, T.; Gan, J.; Erofeeva, S.Y. Numerical study of the tide and tidal dynamics in the South China Sea. Deep Sea Res. Part I Oceanogr. Res. Pap. 2008, 55, 137–154. [Google Scholar] [CrossRef]
  51. Jan, S.; Yang, Y.-J.; Wang, J.; Mensah, V.; Kuo, T.-H.; Chiou, M.-D.; Chern, C.-S.; Chang, M.-H.; Chien, H. Large variability of the Kuroshio at 23.75°N east of Taiwan. J. Geophys. Res. Oceans 2015, 120, 1825–1840. [Google Scholar] [CrossRef]
  52. Mori, N.; Kato, M.; Kim, S.; Mase, H.; Shibutani, Y.; Takemi, T.; Tsuboki, K.; Yasuda, T. Local amplification of storm surge by Super Typhoon Haiyan in Leyte Gulf. Geophys. Res. Lett. 2014, 41, 5106–5113. [Google Scholar] [CrossRef] [Green Version]
  53. Sepúlveda, I.; Tozer, B.; Haase, J.S.; Liu, P.L.; Grigoriu, M. Modeling uncertainties of bathymetry predicted with satellite altimetry data and application to tsunami hazard assessments. J. Geophys. Res. Solid Earth 2020, 125, e2020JB019735. [Google Scholar] [CrossRef]
  54. Kowalik, Z.; Murty, T.S. Chapter III Two-Dimensional Numerical Models. In Numerical Modeling of Ocean Dynamics; World Scientific Publishing Co. Pte. Ltd.: Singapore, 1993; pp. 105–215. [Google Scholar] [CrossRef]
Figure 1. Sketch of a nested-grid domain in the grid size ratio of 3:1 (a) at the upper left corner of a fine grid; (b) at the lower right corner of a fine grid. Circles indicate the free surface elevations (black–the coarse grid; blue–the fine grid). Arrows show the volume-flux components of the x- and y-directions (black–the coarse grid; blue–the fine grid).
Figure 1. Sketch of a nested-grid domain in the grid size ratio of 3:1 (a) at the upper left corner of a fine grid; (b) at the lower right corner of a fine grid. Circles indicate the free surface elevations (black–the coarse grid; blue–the fine grid). Arrows show the volume-flux components of the x- and y-directions (black–the coarse grid; blue–the fine grid).
Water 14 00547 g001
Figure 2. Grid-nesting procedure between the outer and inner grid in Δ t o u t / Δ t i n = 2 . We note that Δ t o u t and Δ t i n indicate the time step of the outer and inner grids, respectively. η is the free surface elevation, ( P , Q ) are the volume-flux components, and the superscript n denotes the time label. This figure is replotted from [26].
Figure 2. Grid-nesting procedure between the outer and inner grid in Δ t o u t / Δ t i n = 2 . We note that Δ t o u t and Δ t i n indicate the time step of the outer and inner grids, respectively. η is the free surface elevation, ( P , Q ) are the volume-flux components, and the superscript n denotes the time label. This figure is replotted from [26].
Water 14 00547 g002
Figure 3. One-dimensional illustration for the moving boundary scheme: (a) the shoreline stops between the grid cells i and i+1; (b) the shoreline stops between the grid cells i+1 and i+2. Circles indicate the cell centers (blue–wet cells; black–dry cells); rectangles imply the cell edges for volume-flux components (green–nonzero flux; red–zero flux); the blue-shaded object identifies the water body.
Figure 3. One-dimensional illustration for the moving boundary scheme: (a) the shoreline stops between the grid cells i and i+1; (b) the shoreline stops between the grid cells i+1 and i+2. Circles indicate the cell centers (blue–wet cells; black–dry cells); rectangles imply the cell edges for volume-flux components (green–nonzero flux; red–zero flux); the blue-shaded object identifies the water body.
Water 14 00547 g003
Figure 4. (a) Top view of the wave basin and the circular island. The blue asterisks indicate the locations of the wave gauges. The black arrow implies the incident wave direction. The dashed red line shows the domain of Grid 02. (b) Side view of the circular island on the x-z plane along the centerline of the circular island. The blue line indicates the still water surface.
Figure 4. (a) Top view of the wave basin and the circular island. The blue asterisks indicate the locations of the wave gauges. The black arrow implies the incident wave direction. The dashed red line shows the domain of Grid 02. (b) Side view of the circular island on the x-z plane along the centerline of the circular island. The blue line indicates the still water surface.
Water 14 00547 g004
Figure 5. Snapshots for the computed free surface elevations of Grid 02 on the frontal side of the circular island (A/h = 0.091) at t = (a) 8.5, (b) 9.0, (c) 9.5, (d) 10.0, (e) 10.5, (f) 11.0, (g) 11.5, and (h) 12.0 sec. The incident solitary wave propagates along the +y-direction. The color bar indicates the free surface elevations in the unit of m. The gray-shaded object shows the circular island.
Figure 5. Snapshots for the computed free surface elevations of Grid 02 on the frontal side of the circular island (A/h = 0.091) at t = (a) 8.5, (b) 9.0, (c) 9.5, (d) 10.0, (e) 10.5, (f) 11.0, (g) 11.5, and (h) 12.0 sec. The incident solitary wave propagates along the +y-direction. The color bar indicates the free surface elevations in the unit of m. The gray-shaded object shows the circular island.
Water 14 00547 g005aWater 14 00547 g005b
Figure 6. Snapshots for the computed free surface elevations of Grid 02 on the lee side of the circular island (A/h = 0.091) at t = (a) 11.0, (b) 11.5, (c) 12.5, (d) 13.5, (e) 14.0, (f) 14.5, (g) 15.0, and (h) 16.0 sec. The incident solitary wave propagates along the +y-direction. The color bar indicates the free surface elevations in the unit of m. The gray-shaded object shows the circular island.
Figure 6. Snapshots for the computed free surface elevations of Grid 02 on the lee side of the circular island (A/h = 0.091) at t = (a) 11.0, (b) 11.5, (c) 12.5, (d) 13.5, (e) 14.0, (f) 14.5, (g) 15.0, and (h) 16.0 sec. The incident solitary wave propagates along the +y-direction. The color bar indicates the free surface elevations in the unit of m. The gray-shaded object shows the circular island.
Water 14 00547 g006aWater 14 00547 g006b
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.
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.
Water 14 00547 g007
Figure 8. Comparison between the maximum computed and measured runup heights (A/h = 0.091). Both numerical results and measured data are projected onto the circular island. Note that the incident solitary wave propagates along the +y-direction from the bottom boundary. The red line shows the computed maximum runup heights; the blue asterisks indicate the measured runup heights.
Figure 8. Comparison between the maximum computed and measured runup heights (A/h = 0.091). Both numerical results and measured data are projected onto the circular island. Note that the incident solitary wave propagates along the +y-direction from the bottom boundary. The red line shows the computed maximum runup heights; the blue asterisks indicate the measured runup heights.
Water 14 00547 g008
Figure 10. Computed storm surges of Domain D01: (a) at 21:00 UTC on 7 November 2013; (b) at 23:30 UTC on 7 November 2013; (c) at 00:00 UTC on 8 November 2013; (d) at 00:30 UTC on 8 November 2013. The color shading indicates the free surface elevations in m, and the black arrows show the 10-m wind directions.
Figure 10. Computed storm surges of Domain D01: (a) at 21:00 UTC on 7 November 2013; (b) at 23:30 UTC on 7 November 2013; (c) at 00:00 UTC on 8 November 2013; (d) at 00:30 UTC on 8 November 2013. The color shading indicates the free surface elevations in m, and the black arrows show the 10-m wind directions.
Water 14 00547 g010
Figure 11. Computed storm surges of Domain D03: (a) at 23:00 UTC on 7 November 2013; (b) at 23:30 UTC on 7 November 2013; (c) at 00:00 UTC on 8 November 2013; (d) at 00:30 UTC on 8 November 2013. The color shading indicates the free surface elevations in m. The blue circle shows the location of the Tacloban Airport.
Figure 11. Computed storm surges of Domain D03: (a) at 23:00 UTC on 7 November 2013; (b) at 23:30 UTC on 7 November 2013; (c) at 00:00 UTC on 8 November 2013; (d) at 00:30 UTC on 8 November 2013. The color shading indicates the free surface elevations in m. The blue circle shows the location of the Tacloban Airport.
Water 14 00547 g011
Figure 12. Computed storm-induced currents of Domain D03: (a) at 23:00 UTC on 7 November 2013; (b) at 23:30 UTC on 7 November 2013; (c) at 00:00 UTC on 8 November 2013; (d) at 00:30 UTC on 8 November 2013. The color shading indicates the flow speed in m/s. The black arrows show the flow velocities in m/s. The red circle shows the location of the Tacloban Airport.
Figure 12. Computed storm-induced currents of Domain D03: (a) at 23:00 UTC on 7 November 2013; (b) at 23:30 UTC on 7 November 2013; (c) at 00:00 UTC on 8 November 2013; (d) at 00:30 UTC on 8 November 2013. The color shading indicates the flow speed in m/s. The black arrows show the flow velocities in m/s. The red circle shows the location of the Tacloban Airport.
Water 14 00547 g012
Figure 13. (a) Maximum storm surges of Domain D02. The color shading indicates the maximum storm surges in m. The solid black line indicates the coastline in the numerical model. (b) Maximum flood depth of Domain D03. The color shading indicates the maximum flood depth in m. The solid black line indicates the coastline in the numerical model. The blue circle shows the location of the Tacloban DZR Airport. The green circles imply the locations of measured flood depths (see Table 5). The yellow rectangles with the dashed black line mark inundation areas of [44].
Figure 13. (a) Maximum storm surges of Domain D02. The color shading indicates the maximum storm surges in m. The solid black line indicates the coastline in the numerical model. (b) Maximum flood depth of Domain D03. The color shading indicates the maximum flood depth in m. The solid black line indicates the coastline in the numerical model. The blue circle shows the location of the Tacloban DZR Airport. The green circles imply the locations of measured flood depths (see Table 5). The yellow rectangles with the dashed black line mark inundation areas of [44].
Water 14 00547 g013
Figure 14. Time series of computed storm surges at specified numerical gauges: Basey, Tacloban, Palo, Tanauan, and Dulag. The y-axis indicates storm surge height in m, and the x-axis shows the time from 12:00 UTC on 7 November 2013 to 12:00 UTC on 8 November 2013. The locations of these numerical gauges can be found in Table 3.
Figure 14. Time series of computed storm surges at specified numerical gauges: Basey, Tacloban, Palo, Tanauan, and Dulag. The y-axis indicates storm surge height in m, and the x-axis shows the time from 12:00 UTC on 7 November 2013 to 12:00 UTC on 8 November 2013. The locations of these numerical gauges can be found in Table 3.
Water 14 00547 g014
Figure 15. Maximum storm surges of the (a) nonlinear equation model with the moving shoreline; (b) nonlinear equation model with the fixed shoreline; and (c) linear equation with the fixed shoreline. The color shading indicates the maximum storm surges in m. The solid black line shows the shoreline of COMCOT-SURGE.
Figure 15. Maximum storm surges of the (a) nonlinear equation model with the moving shoreline; (b) nonlinear equation model with the fixed shoreline; and (c) linear equation with the fixed shoreline. The color shading indicates the maximum storm surges in m. The solid black line shows the shoreline of COMCOT-SURGE.
Water 14 00547 g015
Figure 16. Storm-induced current fields by the nonlinear equation model with the fixed shoreline at (a) 00:00 UTC on 8 November 2013 and (b) 00:30 UTC on 8 November 2013. Storm-induced current fields by the linear equation model with the fixed shoreline at (c) 00:00 UTC on 8 November 2013 and (d) 00:30 UTC on 8 November 2013. The color shading indicates the flow speed in m/s. The black arrows show the flow velocities in m/s.
Figure 16. Storm-induced current fields by the nonlinear equation model with the fixed shoreline at (a) 00:00 UTC on 8 November 2013 and (b) 00:30 UTC on 8 November 2013. Storm-induced current fields by the linear equation model with the fixed shoreline at (c) 00:00 UTC on 8 November 2013 and (d) 00:30 UTC on 8 November 2013. The color shading indicates the flow speed in m/s. The black arrows show the flow velocities in m/s.
Water 14 00547 g016
Figure 17. Clock time versus thread numbers of the workstation. It is noted that 2 threads are equal to 1 CPU. The blue rectangles indicate the clock time using different numbers of threads. The red line is the curve-fitted result: y = 42.16 + 196.3 × exp   0.365 x , and x and y here imply the x and y-axis of this figure.
Figure 17. Clock time versus thread numbers of the workstation. It is noted that 2 threads are equal to 1 CPU. The blue rectangles indicate the clock time using different numbers of threads. The red line is the curve-fitted result: y = 42.16 + 196.3 × exp   0.365 x , and x and y here imply the x and y-axis of this figure.
Water 14 00547 g017
Table 1. Wave gauge locations around the circular island.
Table 1. Wave gauge locations around the circular island.
Gauge NumberX Location
(Unit: m)
Y Location
(Unit: m)
G615.009.40
G915.0010.40
G1617.5813.00
G2215.0015.60
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.
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.
DomainSea-Level Pressure/Wind Shear StressAdvection TermCoriolis Force TermBottom FrictionHorizontal Eddy Diffusion TermMoving Boundary Scheme
D01OXOOOX
D02OOOOOX
D03OOXOXO
Table 5. Measured and predicated flood depths at P1, P2, and P3. The measured flood depths are from [44].
Table 5. Measured and predicated flood depths at P1, P2, and P3. The measured flood depths are from [44].
Field Survey PointsLongitude
(Unit: °E)
Latitude
(Unit: °N)
Measured Flood Depth
(Unit: m)
Predicted
Flood Depth
(Unit: m)
P1125.024711.20055.93.095
P2125.022411.22713.53.783
P3125.000411.24573.50.891
Table 6. Clock time with the corresponding usage of the thread number in the workstation. It is noted here that 1 CPU has 2 threads.
Table 6. Clock time with the corresponding usage of the thread number in the workstation. It is noted here that 1 CPU has 2 threads.
Thread
Number
124812162024
Clock Time
(Unit: min)
185.47123.5890.8766.1737.2243.9539.8340.11
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Tsai, Y.-L.; Wu, T.-R.; Yen, E.; Lin, C.-Y.; Lin, S.C. Parallel-Computing Two-Way Grid-Nested Storm Surge Model with a Moving Boundary Scheme and Case Study of the 2013 Super Typhoon Haiyan. Water 2022, 14, 547. https://doi.org/10.3390/w14040547

AMA Style

Tsai Y-L, Wu T-R, Yen E, Lin C-Y, Lin SC. Parallel-Computing Two-Way Grid-Nested Storm Surge Model with a Moving Boundary Scheme and Case Study of the 2013 Super Typhoon Haiyan. Water. 2022; 14(4):547. https://doi.org/10.3390/w14040547

Chicago/Turabian Style

Tsai, Yu-Lin, Tso-Ren Wu, Eric Yen, Chuan-Yao Lin, and Simon C. Lin. 2022. "Parallel-Computing Two-Way Grid-Nested Storm Surge Model with a Moving Boundary Scheme and Case Study of the 2013 Super Typhoon Haiyan" Water 14, no. 4: 547. https://doi.org/10.3390/w14040547

APA Style

Tsai, Y. -L., Wu, T. -R., Yen, E., Lin, C. -Y., & Lin, S. C. (2022). Parallel-Computing Two-Way Grid-Nested Storm Surge Model with a Moving Boundary Scheme and Case Study of the 2013 Super Typhoon Haiyan. Water, 14(4), 547. https://doi.org/10.3390/w14040547

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop