A typhoon-induced storm surge numerical model with GPU acceleration based on an unstructured spherical centroidal Voronoi tessellation grid

1. Introduction

Storm surge is an abnormal rise and fall of the sea surface caused by strong atmospheric disturbances, which generally include strong winds and sudden pressure changes (Muis et al., 2016). Storm surge is often the marine disaster that poses the greatest threat to life and property in coastal areas. In particular, the combination of storm surges from tropical cyclones and high astronomical tides often causes devastation to coastal areas. Globally, many storm surge disasters associated with tropical cyclones have caused heavy casualties. For example, in 1959, Typhoon Vera struck Japan and caused over 5 000 deaths, most of which were related to the 3.5-m typhoon-induced storm surge in the Isewan area (Higaki et al., 2009). In 1970, the Bhola cyclone killed more than 300 000 people in East Pakistan (now Bangladesh) (Frank and Husain, 1971; Dube et al., 1997). In July 1980, Typhoon Joe (the highest storm surge was close to 6 m) made landfall in Guangdong Province, China and eventually caused 455 deaths (Yu et al., 2020). In 2005, Hurricane Katrina made landfall in New Orleans and caused approximately 1 100 deaths, many of which were directly or indirectly caused by storm surges (Rappaport, 2014).

Accurate and timely issuance of storm surge warnings to take appropriate countermeasures is an important means to reduce storm surge-related losses. There are two types of storm surge forecasting methods: empirical statistical forecasting methods and dynamic numerical simulation methods. Dynamic numerical simulation methods are the mainstream methods used in countries worldwide for storm surge forecasting. Numerical storm surge forecasting uses the numerical weather forecast or empirical typhoon model to provide the wind field and pressure over the sea surface, and it numerically solves the seawater motion equations to predict the storm surge in a specific sea area. Early storm surge models mostly use structured grids, such as the Sea, Lake & Overland Surges from Hurricanes (SLOSH) model in the United States and the storm surge model in Japan (Higaki et al., 2009). However, there are two problems with structured grids. First, the nearshore boundary is roughly defined, which affects the accuracy of the nearshore calculation results. Second, the high-resolution global simulation has large computational cost. Therefore, high-resolution simulation often requires nested refinement. However, nesting of simulation areas with different resolutions leads to instability at the boundary. To overcome the shortcomings of structured grids, unstructured grids have been adopted by many numerical storm surge models in the past 20 years, and the finite element method or the finite volume method has been used to solved these models (Kohno et al., 2018). Unstructured grid models, such as the triangular grid used by finite-volume coastal ocean model (FVCOM; Chen et al., 2003) and advanced circulation model (ADCIRC; Dietrich et al., 2011), that can well depict complex coastlines and can use high resolutions for nearshore areas and low resolutions for offshore areas.

Among the unstructured grids, spherical centroidal Voronoi tessellation (SCVT) is gradually attracting the attention of scientists. The SCVT consists of Voronoi diagrams which can handle arbitrary boundaries and refinement well. The SCVT grid can use transition refinement to generate high-resolution grids in key areas without nesting and low-resolution grids in non-key areas to save computational resources, thus ensuring high-resolution simulation of key areas while reducing computation time. Compared with nested refinement, transition refinement avoids rigid transitions and has better global consistency. Storm surge disasters occur in coastal areas, so SCVT grids that can support the transition refinement is very needed to more accurately describe the coastline.

Many scientists and institutions have developed marine and atmospheric numerical forecast models based on the SCVT grid. Thuburn et al. (2009) and Ringler et al. (2010) provided a spatial discrete scheme based on arbitrarily structured C-grids, known as Thuburn–Ringler–Skamarock–Klemp (TRiSK). Thuburn et al. (2009) solved the expression of the Coriolis force in arbitrarily structured grids, described the Coriolis force in arbitrarily structured C-grids based on the simplified vorticity equation, and accurately represented the geostrophic modes while ensuring that the Coriolis force does not perform work. Ringler et al. (2010) ensures strict conservation of mass, absolute vorticity, and instantaneous energy. And TRiSK finally developed into the shallow-water dynamic core of the Model for Prediction Across Scales (MPAS) by the National Center for Atmospheric Research (NCAR) and Los Alamos National Laboratory (LANL; Skamarock et al., 2012). Zhou et al. (2020) further introduced square conservation into the TriSK framework, that enables the model to maintain energy conservation, and developed a new shallow-water dynamic framework called TRiSK-based Multiple-Conservation dynamical cORE (TMCORE) by adopting a temporal integration scheme with square conservation (the new dynamic algorithm can strictly ensure the global conservation of mass, energy, and absolute vorticity without affecting the convergence accuracy), and introduced some high-precision advection algorithms.

Based on the above research, we believe that the SCVT grid has better application prospects in storm surge models. So We developed a numerical storm surge model based on TriSK and TMCORE framework to take advantage of SCVT grid, upgrade the National Marine Environmental Forecasting Center (NMEFC) storm surge model from a structured grid to an unstructured grid, and further improve the level of storm surge forecasting. The model is based on shallow water equations in vector-invariant form, and is discretized by Arakawa C grid. We applied it to storm surge simulations in nearshore areas of China, selected typical typhoon storm surge processes for numerical experiments, and analyzed the numerical simulation results in order to investigate the applicability of the model.

The storm surge is closely related to the forecast results of typhoon. It is necessary to consider dozens to hundreds of scenarios such as different tracks and intensities of typhoon, carry out ensemble forecast of typhoon storm surge and provide probability forecast. In the daily operational forecast of NMEFC, it is necessary to complete the ensemble forecast of 145 scenarios within half an hour based on the latest typhoon forecast results, which can be provided to forecasters for reference in time. At present, open-source unstructured grid numerical models such as ADCICR and FVCOM are all based on multi-core CPU parallel acceleration, which cannot meet such computing timeliness requirements. Using graphics processing unit (GPU) to accelerate computing can solve this problem by utilizing the characteristics of GPU computing cores, so our model has developed a GPU version.

The remaining contents of this paper are structured as follows. The second section introduces the dynamic framework of the numerical model in detail, including the numerical method, open boundary conditions, and the typhoon-driven model. The third section introduces the discretization of the model. The fourth section introduces GPU-accelerated model computing through the use of OpenACC directives. The fifth section reports the hindcast results of typical storm surges. The last section is the conclusion and discussion.

2. Dynamic framework

2.1 Governing equations

Storm surges are mainly caused by the effects of wind setup due to strong onshore winds over the sea surface and the inverted barometer effect associated with pressure drops in low-pressure systems. The current mainstream storm surge model adopts the two-dimensional shallow water equations plus the corresponding storm surge forcing terms, including wind stress, pressure gradient force, and bottom friction terms. Based on the TriSK and TMCORE framework, we introduce the storm surge driving terms into the vector-invariant shallow water equations. The vector-invariant shallow water equations given by Ringler et al. (2010), and we change the forecast variable from the total water depth to the storm surge.

where $ h $ is the water level change calculated from the mean sea level, i.e., storm surge; $ H = h + {H_0} $ is the total water depth, and $ {H_0} $ is the water depth at the mean sea level; $ {\boldsymbol{u}}$ is the velocity, and ${{\boldsymbol{u}}^ \bot } = {\boldsymbol{k}} \times {\boldsymbol{u}}$ is the tangential velocity, ${\boldsymbol{k}}$ is the unit vector, which points in the local vertical direction. $ q $ is the potential vorticity; $K = \dfrac{{{{\left| { {\boldsymbol{u}}} \right|}^2}}}{2}$ is the kinetic energy; ${p_{\rm{a}}}$ is the atmospheric pressure at the sea level; ${{\boldsymbol{\tau}} _{\mathbf{s}}}$ is wind stress, ${{\boldsymbol{\tau}} _{\mathbf{b}}}$ is the bottom friction stress, ${{\boldsymbol{\tau}} _{\mathbf{s}}}$ and ${{\boldsymbol{\tau}} _{\mathbf{b}}}$ are calculated using quadratic parameterization (${{\boldsymbol{\tau}} _{\rm{s}}} = {C_{\rm{d}}}{{{\rho}}_{\rm{a}}}{ {\boldsymbol{u}}_{\mathbf{a}}}\left| {{{ {\boldsymbol{u}}}_{\mathbf{a}}}} \right|$, ${\tau _{\rm{b}}} = c{\rho _{\rm{w}}}{\boldsymbol{u}}\left| { {\boldsymbol{u}}} \right|$); ${ {\boldsymbol{u}}_{\mathbf{a}}}$ is the wind vector at 10 m above the sea surface; ${\rho _{\rm{w}}}$ and ${\rho _{\rm{a}}}$ are seawater density and air density; $ c $ is the bottom friction coefficient, which can be taken as $ c = 2.6 \times {10^{ - 3}} $; ${C_{\rm{d}}}$ is the wind drag coefficient, which can be taken as ${C_{\rm{d}}} = (0.8 + 0.065 \times \left| {{{ {\boldsymbol{u}}}_{\mathbf{a}}}} \right|) \times {10^{ - 3}}$; $ g $ is the gravitational acceleration.

2.2 Initial and open boundary condition

The initial condition is not important in the storm surge simulation. The error in the initial condition can be overcome by the viscosity and the forcing very quickly. So, $ h $ and ${\boldsymbol{ u}}$ are set to zero at the initial time.

Storm surge is generally simulated in a limited area, and an open boundary is intentionally introduced to the simulation area. Therefore, it is necessary to reasonably handle the impact of the open boundary. The introduced open boundary should not reflect the waves propagating outward from the simulation area. To satisfy this need, a sponge boundary condition is used for the open boundary to absorb the waves propagating outward. The specific implementation method is as follows: create a seven-layer transition relaxation region at the open boundary, with the outermost layer as the given boundary condition, the innermost layer as the free development, and the intermediate transition layers relaxed according to the relaxation coefficient (Fig. 1).

Figure  1.  Schematic diagram of the open boundary.

DownLoad: Full-Size Img PowerPoint

The relaxation coefficient is related to the distance between the intermediate transition layer and the open boundary. A simple empirical value is as follows:

where ${h_{{n}}}$ represents the water level on the n-th layer relax cell. $ {h_{{\rm{inner}}}} $ and $ {h_{{\rm{outer}}}} $ represent the water level on the innermost and outermost layers respectively. ${w_{{n}}}$ represents the relaxation coefficient on the n-th layer.

The water level ${h_{{\rm{outer}}}}$ is the hydrostatic water level. When the atmospheric pressure decreases by 1 hPa, the water level increases by approximately 1 cm.

where $ {p_{\rm{a}}} $ represents the atmospheric pressure at the sea level; $ {P_{\rm{e}}} $ is the pressure at infinity (ambient pressure; 1010 hPa for the northwestern Pacific Ocean)

The rigid boundary condition is taken as the land boundary condition, and the normal velocity is 0.

2.3 Meteorological forcing

A storm surge model requires fields of surface wind and atmospheric pressure as external forcing and the quality of the atmospheric driving field determines the quality of storm surge simulation results. This study adopts the Holland typhoon model widely used in storm surge models.

where $ r $ is the distance from the target point to the center of the typhoon; ${R_{\rm{max}}}$ is the radius of maximum wind of the typhoon, and ${R_{\rm{max}}}$ is generally calculated using the empirical formula; $ {P_{\rm{o}}} $ is the central pressure of the typhoon; $ {P_{\rm{e}}} $ is the ambient pressure; $ B $ is the fitting parameter of the Holland model, which can be calculated using the empirical formula ($ B = 1.5 + (980 - {P_{\rm{o}}})/120 $) given by Hubbert et al. (1991). The air pressure $ P(r) $ at any point $ r $ can be obtained according to this empirical formula. Then, the entire typhoon wind field is calculated according to the gradient wind formula.

3. Discrete system

3.1 Notation

The discrete system requires the definition of three location (Fig. 2). The subscript $ i $ represents that the variable is on the centroid cell $ {X_i} $, and subscript $ e $ represents that the variable is on the edge $ {X_{{e}}} $. The ocean current $ {\boldsymbol{u}} $ is located on the edge of each cell $ {X_{{e}}} $, the water depth $ h $ is located at the centroid of each cell ${X_{{i}}}$, and the vorticity $ q $ is located at the vertices of each cell $ {X_{{v}}} $.The cells are nominally located at the Voronoi generating points, which, for centroidal Voronoi tessellations, are the mass centoids of the Voronoi cells with respect to a density function. The edges are nominally located at the midpoints of edges and the vertices are located at the corner of cell. The mass points are in the center of each cell, the velocity points are on the edges of each cell, and vorticity points are on the vertices of each cell.

Figure  2.  Variable locations on the SCVT grid.

DownLoad: Full-Size Img PowerPoint

In Fig. 2, ${F_{{e}}} = {\hat H_{{e}}} {\boldsymbol{u}}$ and ${F_{{e}}}^ \bot = {\hat H_{{e}}}{ {\boldsymbol{u}}^ \bot }$ represent the normal and tangential components of the water level flux per unit length on the edge $ {X_{{e}}} $, respectively. $ {\hat H_{{e}}} $ represents the water depth $ H $ reconstructed on the edge from its original position(cell centroid), and similarly, $ {\hat q_{{e}}} $ represents the potential vorticity $ q $ reconstructed on the edge from its original position (vertex). ${E_{{i}}} = g{h_{{i}}} + {K_{{i}}} + {\left(\dfrac{{{p_{\rm{a}}}}}{{{\rho _{\rm{w}}}}}\right)_{{i}}}$ is the sum of the water level, kinetic energy, and air pressure on the cell centroid ${X_{{i}}}$.

3.2 Spatial discretization

The finite volume method is used for spatial discretization in the horizontal direction. The variable configuration of the Arakawa C grid is used for discretization on the SCVT grid.

With the definition in Section 3.1, Eq. (1) and Eq. (2) are rewritten in a discrete form as follows:

In the governing Eq. (6) and Eq. (7), divergence, gradient and vorticity need to be calculated, so corresponding spatial discrete operators are required. The divergence is defined at the centroid ${X_{{i}}}$, and the divergence is the sum of the product of the flux per unit length on each edge $ {F_{{e}}} $ and the length of each edge $ {l_{{e}}} $ divided by the area of cell ${A_{{i}}}$.

where ${{e}} \in {{EC}}({{i}})$ represents the set of edges that on the cell. where ${n_{{{e, }}i}}$ is a indicator function that defines the direction of flux. ${n_{{{e, }}i}} =1$ represents the positive normal direction, namely, that the ocean current perpendicular to the edge flows toward the outside of cell, while when ${n_{{{e, }}i}} = - 1$ represents the positive normal direction, namely, that the ocean current perpendicular to the edge flows toward the inner side of cell.

The gradient is defined at the edge $ {X_{{e}}} $ as follows:

where ${{i}} \in {\rm{CE}}({{e}})$ represents the set of cells that share the edge. ${E_{{e}}}$ represents the scalar reconstructed on the edge from its original position ${E_{{i}}}$ (cell centroid).

The vorticity is defined at the vertex ${X_{\rm{v}}}$ as follows:

where ${{e}} \in {\rm{EV}}({{v}})$ represents the set of edges that share the vertex. The indicator function ${t_{{{e, v}}}}$ always points to the left side of ${n_{{{e, i}}}}$. If ${\boldsymbol{k}} \times {n_{{{e, i}}}}$ is directed toward vertex, then ${t_{{{e, v}}}} = 1$; otherwise, ${t_{{{e, v}}}} = - 1$. ${\boldsymbol{k}}$ is the unit vector, which points in the local vertical direction.

The spatial discretized shallow water Eq. (6) and Eq. (7) are written as follows:

3.3 Temporal integration scheme

The control equation is discretized in time with the explicit integration scheme, we use the third-order Runge‒Kutta (RK3) scheme.

In order to discuss the temporal integration scheme, the spatial discretization in Eq. (11) and Eq. (12) is simplified as $ L(h,{\boldsymbol{u}}) $ and $ F(h,{\boldsymbol{u}}) $ , respectively, and the following equations are obtained:

The temporal integration scheme is RK3. Take continuity Eq. (13) as an example to illustrate the discrete form of the integration scheme.

where $ h_i^t $ represents the water level at the previous time step, and $ h_i^{t + 1} $ represents the water level with forward integration of $ \Delta t $.

4. GPU accelerated computing

When forecasting typhoon storm surge, different tracks and intensities of typhoon shall be considered, and dozens of ensemble members shall be predicted. Therefore, parallel acceleration is required to meet the timeliness requirements of forecasting. In general, atmospheric or oceanic numerical models are mainly based on MPI technology and use multi-core CPUs for parallel acceleration. In recent years, with the maturity of GPU hardware, the technology of GPU accelerated computing is also gradually developing. For numerical computing tasks of the same scale, the power consumption of GPU is only 1/30 of that of traditional high-performance computer (HPC), and the purchase cost of GPU hardware is lower. The time and spatial discretization of our model adopts an explicit calculation scheme, which is a computationally intensive task. GPU devices generally have tens of thousands of computing cores, which are very suitable for computing intensive tasks.

In order to make the numerical mode written in Fortran code support Nvidia’s GPU hardware, there are two mainstream solutions. The first solution is to rewrite the original fortran code into the Compute Unified Device Architecture Fortran (CUDA-Fortran) version. This solution has high acceleration performance, but the code development workload is relatively large, and most numerical model experts cannot be proficient, which makes it difficult to upgrade and maintain the code. The second solution is to add some OpenACC guidance statements on the original fortran code. OpenACC is a directive-based programming model designed to provide a simple yet powerful approach to accelerators without significant programming effort. So it is easy for non-professionals to master, and it is convenient for maintenance and upgrade of the code. Therefore, we use OpenACC directives to accelerate the numerical model.

OpenACC directives are much like OpenMP directives. The main work of OpenACC acceleration lies in two aspects. On the one hand, it is to add some directives to the code with a large amount of calculation. For example, the subroutine for calculating gradient can achieve GPU acceleration only by adding some directives as follow:

subroutine grad_operator(f_cell, grad)

!$acc parallel present(variables needed to compute gradient)

!$acc loop independent

do i = 1, nEdges

Calculating the gradient…

end do

!$acc end parallel

end subroutine grad_operator

where “!$acc parallel present” indicates that the variable has been copied to the GPU device and can be used directly; “!$acc loop” indicates that the do loop uses GPU calculation; “!$acc end parallel” indicates the end of the calculation on the GPU. On the other hand, do a good job of copying data between the CPU host and the GPU device. The variables of the model initial field and atmospheric drive field need to be copied from the CPU host to the GPU device during model initialization. When the integrated data needs to be saved, the variables need to be copied from the GPU device to the CPU host.

On the basis of serial code, we added only 500 lines of OpenACC directives. Compared with the serial version, the GPU parallel acceleration ratio on an NVIDIA V100 reaches 80−100 times. In the simulation of China’s coastal waters (the area shown in Fig. 4), the 3-d forecast only takes 110 s. However, the operational forecast system based on the ADCIRC model uses a computing configuration of 28 CPU cores on a single node, and it takes 10 minutes to complete a similar task. The operational organizations purchase with 2−3 GPU graphics, the ensemble forecasting task of typhoon storm surge can be completed within half an hour.

5. Model evaluation

5.1 Description of typhoon case

To evaluate the simulation effect of the model, the typhoon-induced storm surge processes that have had the great impact on China’s offshore in recent years were selected for simulation, and the accuracy of the model was evaluated (Fig. 3). The typhoon data were from the nearest path dataset of the China Meteorological Administration (CMA) Tropical Cyclone Data Center (Lu et al., 2021). Six typical typhoons in China shown in Fig. 3 were selected for hindcast experiments. They all made landfall in China and caused large water surges, and their paths essentially covered the typical paths of typhoons that make landfall in China.

Figure  3.  Six typical typhoons that made landfall in China. The star (★) is the location of the tide gauges.

DownLoad: Full-Size Img PowerPoint

5.2 Model configuration

A spherical centroidal Voronoi tessellation (SCVT) is a Voronoi tessellation whose generators coincide with the mass centroids of their corresponding Voronoi cells. There is a duality relation between Delaunay Triangulations and Voronoi Diagrams, each vertex of the Delaunay triangulation is a site in the Voronoi diagram which gets associated with its Voronoi cell. So we can use JIGSAW-GEO tool to generate spherical triangle mesh, and then use the duality relation to form SCVT mesh (Engwirda, 2017).

Figure  4.  The resolution of the SCVT grid used in the simulation transitions from 20 km offshore to 2 km nearshore.

DownLoad: Full-Size Img PowerPoint

The grid configuration is shown in Fig. 4. The spatial resolution of the grid transitions from 20 km offshore to 2 km nearshore, covering the whole China offshore. The number of grid cells is 124 588. This grid ensures the simulation accuracy of the nearshore areas while saving computational resources.

In NMEFC, there is an operational typhoon storm surge forecasting system based on the ADCIRC model. The spatial resolution of the grid transitions from 25 km offshore to 1 km nearshore. In order to compare with our model, we use the same Holland empirical wind field to conduct hindcast experiments. The empirical parameters such as bottom friction and wind drag coefficient are consistent.

5.3 Result discussion

The distributions of the maxinum storm surges during typhoon-induced storm surge processes (Fig. 5) show that due to the use of the parametric model wind field as the driving field and due to the typhoon movement direction and land distribution, the right semicircular area of typhoon was dominated by the onshore wind, resulting in generally large storm surges. This distribution characteristic indicates that the model can reasonably simulate typhoon wind field-driven storm surge processes. The six selected typical typhoon-induced storm surge processes are strong, and the areas with the maximum storm surge exceeding 2 m are large. For example, Typhoon 1415 (Kalmaegi) caused large-scale storm surges exceeding 2 m in the coastal areas of Guangxi and on the east coast of Leizhou Peninsula in Guangdong, northern Hainan Island. Typhoon 1909 (Lekima) also brought severe storm surges to the coastal areas of Zhejiang, Shanghai, and Jiangsu as it moved northward.

DownLoad: Full-Size Img PowerPoint

The tide gauges with the largest storm surges during the six typhoon-induced storm surge processes were selected as the representatives, and the time series of storm surges were plotted (Fig. 6). The time series diagrams of storm surges show that the model can simulate typhoon-induced storm surge processes. The simulated maximum storm surges and their occurrence times are consistent with the observation data at the tide gauges, and the simulation results can reflect some precursor waves of storm surges before the typhoon affected these areas.

Figure  6.  Storm surge processes at the representative stations. The black scatter shows the observed values, and the green line and red line show the simulated values of ADCIRC model and our model respectively (unit: cm).

DownLoad: Full-Size Img PowerPoint

The tide gauges with the largest storm surges were selected for comparing the observed and simulated storm surges of the six typhoon-induced storm surge processes (Table 1). The comparison results show that the mean absolute error (MAE) of simulated maxinum storm surges are between −3 cm and 7 cm and that the MAE of simulated phases are mostly within 1 h. The storm surge curves and the results of the hindcast experiments indicate that the model can well characterize the storm surge processes at the tide gauges during the typhoon-affected period. Compared with the operational system built by ADCIRC model, the MAE of the maxinum storm surge and phase simulated by our model is smaller. Of course, it needs to be explained that the empirical parameters set in the hindcast experiment may not be the most suitable for the ADCIRC model, but this result shows that the performance of our model has reached the level of advanced and mature models (such as ADCIRC).

Table  1.  Comparison of simulated and measured maxinum storm surge and phases

Typhoon	Tide gauges	Maximum storm surge/cm				Maximum storm surge error/cm			Phase error/h
		OBS	ADCIRC	NMEFC		ADCIRC	NMEFC		ADCIRC	NMEFC
1323	Longgang	399	381	396		–18	–3		0	1
1415	Nandu	495	518	495		23	0		2	3
1522	Shuidong	232	235	236		3	4		0	0
1614	Shijing	288	269	291		–19	3		1	0
1822	Sanzao	339	347	343		8	4		–1	0
1909	Haimen	312	332	319		20	7		2	1

 | Show Table

DownLoad: CSV

6. Conclusions and discussion

This paper introduces the development of a new-generation numerical storm surge forecast model based on the unstructured SCVT grid. This model is based on the vector-invariant, depth-averaged two-dimensional shallow water equations. The finite volume method is used for spatial discretization, and the RK3 scheme is used for temporal discretization. The Arakawa C grid is used for variable configuration. And the sponge boundary condition is used for the regional open boundary. The use of the unstructured SCVT grid enables refinement in key areas and can well fit the coastline information. In addition, the simulation speed of the model is accelerated by using the openACC-based GPU acceleration technology to meet the timeliness requirements of operational applications. It only takes 37 s to simulate a day in the coastal waters of China, and users can forcast storm surge on their personal computers.

Using a grid configuration with a resolution of 20 km offshore to 2 km nearshore and the Holland typhoon model wind field, the hindcast experiments on the storm surge processes of six typhoons affecting China show that the model can reasonably simulate the distribution characteristics of storm surges under the influence of typhoons, namely, the higher storm surges in the right semicircular areas of the typhoons. The simulated maxinum storm surges and their occurrence times are consistent with the observed data at the representative tide gauge stations, and the mean absolute errors are 3.5 cm and 0.6 h respectively, showing high accuracy and application prospects.

In this study, we mainly considered the ability of the model to simulate strong typhoon-induced storm surges. In the future, we will investigate its application in the simulation of weak typhoon-induced storm surges and use the numerical wind field as the driving field to evaluate the performance of the model for extratropical storm surge simulation. To broaden the application scope of the model, we will continue to develope the model, including the coupling of astronomical tides and floodplains process simulation. Currently, only single GPU acceleration is supported, and multi GPU acceleration will continue to be developed in the future.