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

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

PDF

Yuanyong Gao¹, Fujiang Yu¹^,²^,^*, Cifu Fu¹^,^*, Jianxi Dong¹^,², Qiuxing Liu¹

Acta Oceanologica Sinica | 2024, 43(3) : 40 - 47

Less

Acta Oceanologica Sinica | 2024, 43(3): 40-47

• Physical Oceanography, Marine Meteorology and Marine Physics •

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

Full

Yuanyong Gao¹, Fujiang Yu¹^,²^,^*, Cifu Fu¹^,^*, Jianxi Dong¹^,², Qiuxing Liu¹

Affiliations

¹ National Marine Environmental Forecasting Center, Beijing 100081, China

² Key Laboratory of Marine Hazards Forecasting, Ministry of Natural Resources, Beijing 100081, China

Published: 2024-03-25 doi: 10.1007/s13131-023-2175-9

Outline

Abstract

Less

Storm surge is often the marine disaster that poses the greatest threat to life and property in coastal areas. Accurate and timely issuance of storm surge warnings to take appropriate countermeasures is an important means to reduce storm surge-related losses. Storm surge numerical models are important for storm surge forecasting. To further improve the performance of the storm surge forecast models, we developed a numerical storm surge forecast model based on an unstructured spherical centroidal Voronoi tessellation (SCVT) grid. The model is based on shallow water equations in vector-invariant form, and is discretized by Arakawa C grid. The SCVT grid can not only better describe the coastline information but also avoid rigid transitions, and it has a better global consistency by generating high-resolution grids in the key areas through transition refinement. 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 ensemble forecast. It only takes 37 s to simulate a day in the coastal waters of China. The newly developed storm surge model was applied to simulate typhoon-induced storm surges in the coastal waters of China. The hindcast experiments on the selected representative typhoon-induced storm surge processes indicate that the model can reasonably simulate the distribution characteristics of storm surges. The simulated maximum 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.

Key words

typhoon-induced storm surge / numerical model / GPU acceleration / unstructured grid / spherical centroidal Voronoi tessellation (SCVT)

Cite this Article

Yuanyong Gao, Fujiang Yu, Cifu Fu, Jianxi Dong, Qiuxing Liu. A typhoon-induced storm surge numerical model with GPU acceleration based on an unstructured spherical centroidal Voronoi tessellation grid[J]. Acta Oceanologica Sinica, 2024 , 43 (3) : 40 -47 . DOI: 10.1007/s13131-023-2175-9

Full Text

Less

1 Introduction

Less

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

Less

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.

(1)

$ \frac{{\partial h}}{{\partial t}} + \nabla \cdot (H {\boldsymbol{u}}) = 0 ,$

(2)

$ \frac{{\partial {\boldsymbol{u}}}}{{\partial t}} - qH{{\boldsymbol{u}}^ \bot } = - \nabla (gh + K) - \frac{1}{{{\rho _{\rm{w}}}}}\nabla {p_{\rm{a}}} + \frac{1}{{{\rho _{\rm{w}}}H}}({{\boldsymbol{\tau}} _{\mathbf{s}}} - {{\boldsymbol{\tau}} _{\mathbf{b}}}), $

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

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

(3)

$ \left\{ \begin{aligned} & {h_{{n}}} = {w_{{n}}}{h_{{\rm{inner}}}} + (1 - {w_{{n}}}){h_{{\rm{outer}}}}, \\ & {w_{{n}}} = 1 - 0.1{{n}}, \;\;{{n}} = 1,2,\;\cdots,7 , \end{aligned} \right. $

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.

(4)

$ h = \frac{{{P_{\rm{e}}} - {P_{\rm{a}}}}}{{{\rho _{\rm{w}}}g}} , $

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.

(5)

$ P(r) = {P_{\rm{o}}} + ({P_{\rm{e}}} - {P_{\rm{o}}})\exp \left[ - {({R_{\rm{max}}}/r)^B}\right] , $

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

Less

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.

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:

(6)

$ \frac{{\partial {h_i}}}{{\partial t}} + {[\nabla \cdot {F_e}]_i} = 0 , $

(7)

$ \frac{{\partial {{ {\boldsymbol{u}}}_e}}}{{\partial t}} - {\hat q_e}{F_e}^ \bot = - {[\nabla {E_i}]_e} + \frac{1}{{{\rho _{\rm{w}}}{{\hat H}_e}}}{({{\boldsymbol{\tau}} _{\rm{s}}} - {{\boldsymbol{\tau}} _{\rm{b}}})_i} . $

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}}}$

(8)

$ {[\nabla \cdot {F_{{e}}}]_{{i}}} = \frac{1}{{{A_{{i}}}}}\sum\limits_{{{e}} \;\in\; EC({{i}})} {{n_{{{e}}, i}}{l_{{e}}}{F_{{e}}}} , $

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:

(9)

$ {[\nabla {E_i}]_e} = \frac{1}{{{d_e}}}\sum\limits_{i \;\in\; {\rm{CE}}(e)} { - {n_{e,i}}{E_e}} , $

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:

(10)

$ {\boldsymbol{k}} \cdot {[\nabla \times {\boldsymbol{u}}]_{{v}}} = \frac{1}{{{A_{{v}}}}}\sum\limits_{{{e}} \;\in\; {\rm{EV}}({{v}})} { - {t_{{{e, v}}}}{{\boldsymbol{u}}_{{e}}}{d_{{e}}}} , $

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:

(11)

$ \frac{{\partial {h_{{i}}}}}{{\partial t}} + \frac{1}{{{d_{{e}}}}}\sum\limits_{{{i}} \;\in\; CE({{e}})} { - {n_{{{e, i}}}}{E_{{e}}}} = 0 , $

(12)

$ \begin{split} & \frac{{\partial {{{\boldsymbol{u}}}_e}}}{{\partial t}} - \frac{{{{\hat H}_e}{{ {\boldsymbol{u}}}^ \bot }}}{{{A_v}}}\sum\limits_{e\; \in\; EV(v)} { - {t_{e,v}}{{u}_e}{d_e}} =\\& - \frac{1}{{{d_e}}}\sum\limits_{i\; \in\; CE({{e}})} { - {n_{e,i}}{E_e}} + \frac{1}{{{\rho _{\rm{w}}}{{\hat H}_e}}}{({{\boldsymbol{\tau}} _{\rm{s}}} - {{\boldsymbol{\tau}} _{\rm{b}}})_i} .\end{split}$

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:

(13)

$ \frac{{\partial {h_i}}}{{\partial t}} + L(h,{\boldsymbol{u}}) = 0 , $

(14)

$ \frac{{\partial {{ {\boldsymbol{u}}}_e}}}{{\partial t}} + F(h,{\boldsymbol{u}}) = 0 . $

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

(15)

$ \left\{ \begin{gathered} h_i^{t + 1} = h_i^t + \frac{{\Delta t}}{6}({K_1} + 4{K_2} + {K_3}) , \\ {K_1} = L({h^t},{{\boldsymbol{u}}^{{t}}}), \\ {K_2} = L\left({h^t} + \frac{{\Delta t}}{2}L({h^t},{{\boldsymbol{u}}^{{t}}}),{{{u}}^t} + \frac{{\Delta t}}{2}L({h^t},{{\boldsymbol{u}}^{{t}}})\right) , \\ {K_3} = L\left({h^t} - \Delta t{K_1} + 2\Delta tL({h^t},{{\boldsymbol{u}}^{{t}}}),{u^{{t}}} - \Delta t{K_1} + 2\Delta tL({h^t},{{\boldsymbol{u}}^{{t}}})\right) , \\ \end{gathered} \right. $

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

Less

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

Less

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.

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

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.

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.

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

6 Conclusions and discussion

Less

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.

Acknowledgements

Less

We thank Lilong Zhou and Li Dong for developing the TMCORE model code and providing the benchmark code for the development of this model, and we also thank Ye Yuan for his guidance and help in GPU programming.

Funding

Less

The National Natural Science Foundation of China under contract No. 42076214.

References

Less

Chen Changsheng, Liu Hedong, Beardsley R C. 2003. An unstructured grid, finite-volume, three-dimensional, primitive equations ocean model: application to coastal ocean and estuaries. Journal of Atmospheric & Oceanic Technology, 20(1): 159–186

Dietrich J C, Zijlema M, Westerink J J, et al. 2011. Modeling hurricane waves and storm surge using integrally-coupled, scalable computations. Coastal Engineering, 58(1): 45–65, doi: 10.1016/j.coastaleng.2010.08.001

Dube S K, Rao A D, Sinha P C, et al. 1997. Storm surge in the Bay of Bengal and Arabian Sea: the problem and its prediction. Mausam, 48(2): 283–304, doi: 10.54302/mausam.v48i2.4012

Engwirda D. 2017. JIGSAW-GEO (1.0): JIGSAW-GEO (1.0): Locally orthogonal staggered unstructured grid generation for general circulation modelling on the sphere. Geoscientific Model Development, 10(6): 2117–2140, doi: 10.5194/gmd-10-2117-2017

Frank N L, Husain S A. 1971. The deadliest tropical cyclone in history. Bulletin of the American Meteorological Society, 52(6): 438–445, doi: 10.1175/1520-0477(1971)052<0438:TDTCIH>2.0.CO;2

Higaki M, Hayashibara H, Nozaki F. 2009. Outline of the storm surge prediction model at the Japan Meteorological Agency. RSMC Tokyo-Typhoon Center Technical Review, 2009(11): 25–38

Hubbert G D, Holland G J, Leslie L M, et al. 1991. A real-time system for forecasting tropical cyclone storm surges. Weather and Forecasting, 6(1): 86–97, doi: 10.1175/1520-0434(1991)006<0086:ARTSFF>2.0.CO;2

Kohno N, Dube S K, Entel M, et al. 2018. Recent progress in storm surge forecasting. Tropical Cyclone Research and Review, 7(2): 128–139

Lu Xiaoqin, Yu Hui, Ying Ming, et al. 2021. Western North Pacific tropical cyclone database created by the China meteorological administration. Advances in Atmospheric Sciences, 38(4): 690–699, doi: 10.1007/s00376-020-0211-7

Muis S, Verlaan M, Winsemius H C, et al. 2016. A global reanalysis of storm surges and extreme sea levels. Nature Communications, 7: 11969, doi: 10.1038/ncomms11969

Rappaport E N. 2014. Fatalities in the United States from Atlantic tropical cyclones: new data and interpretation. Bulletin of the American Meteorological Society, 95(3): 341–346, doi: 10.1175/BAMS-D-12-00074.1

Ringler T D, Thuburn J, Klemp J B, et al. 2010. A unified approach to energy conservation and potential vorticity dynamics for arbitrarily-structured C-grids. Journal of Computational Physics, 229(9): 3065–3090, doi: 10.1016/j.jcp.2009.12.007

Skamarock W C, Klemp J B, Duda M G, et al. 2012. A multiscale nonhydrostatic atmospheric model using centroidal voronoi tesselations and C-grid staggering. Monthly Weather Review, 140(9): 3090–3105, doi: 10.1175/MWR-D-11-00215.1

Thuburn J, Ringler T D, Skamarock W C, et al. 2009. Numerical representation of geostrophic modes on arbitrarily structured C-grids. Journal of Computational Physics, 228(22): 8321–8335, doi: 10.1016/j.jcp.2009.08.006

Yu Fujiang, Fu Cifu, Guo Honglin, et al. 2020. Modern Technologies and Application in Storm Surge Forecasting (in Chinese). Beijing: China Science Publishing & Media, 43–47

Zhou Lilong, Feng Jinming, Hua Lijuan. 2020. Extending square conservation to arbitrarily structured C-grids with shallow water equations. Geoscientific Model Development, 13(2): 581–595, doi: 10.5194/gmd-13-581-2020

Appendix

Less

Year 2024 volume 43 Issue 3

PDF

Cite this Article

BibTeX

Article Info

doi: 10.1007/s13131-023-2175-9

Receive Date：2022-11-07
Online Date：2025-11-18
Published：2024-03-25

Article Data

Affiliations

History

Received：2022-11-07
Accepted：2023-02-14

Funding

The National Natural Science Foundation of China under contract No. 42076214.

Affiliations

¹ National Marine Environmental Forecasting Center, Beijing 100081, China

² Key Laboratory of Marine Hazards Forecasting, Ministry of Natural Resources, Beijing 100081, China

Corresponding:

* E-mail: yvfujiang_2022@163.com

fucf@nmefc.cn

References

Share

https://castjournals.cast.org.cn/joweb/aos/EN/10.1007/s13131-023-2175-9

Share to

Scan QR to access full text

Cite this article

BibTeX

Citations

表12种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT

Typhoon	Tide gauges	Maximum storm surge/cm		Maximum storm surge error/cm		Phase error/h
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

Fig. 1. Schematic diagram of the open boundary.

Fig. 2. Variable locations on the SCVT grid.

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

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

Fig. 5. Maximum storm surges during typhoon-induced storm surge processes.

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

Articles: Latest Articles; Most Read; Collections

Updates: Events; News; Multimedia

About: About Us

Contact

No. 86 Xueyuan South Road, Haidian District, Beijing

100081

010-62199257

qkjq@cast.org.cn

Copyright © 2025 China Association for Science and Technology. All rights reserved. For all open access content, the relevant licensing terms apply.
Sponsored by the Office of the Leading Group for Cybersecurity and Informatization of CAST, and supported by Science and Technology Review Publishing House