Super typhoon Rammasun was the strongest recorded typhoon to land in China since 1949. It was generated in the western Pacific Ocean on July 7, 2014, rapidly moved westward through the Philippines into the SCS, gradually intensified to super typhoon status, and successively struck three provinces of southern China ( Figure 1 ). The typhoon track data were downloaded from the Chinese Meteorological Administration (CMA) in the Typhoon track real-time release system of the Zhejiang Water Conservancy Department (http://typhoon.zjwater.gov.cn). These track data showed that Rammasun moved at a velocity of 20 km/h to the northeast coast of Hainan Island with wind speeds exceeding 60 m/s and a central pressure of 910 hPa at 15:00 on July 18. Then, it continued to move westward and passed over the Leizhou Peninsula of Guangdong Province and landed in Guangxi Province. Finally, it rapidly weakened to a severe tropical storm and faded away in mainland China. The track of typhoon Rammasun represents one of the most frequent TC tracks in the coastal areas of SCS. It was estimated that Rammasun caused 7.2 billion dollars of property damages and resulted in the death of 88 people (Yang et al., 2019). Thus, typhoon Rammasun was selected as an example of severe TCs. In this study, because of the hourly time resolution of typhoon track data starting on the 16th at 14:00, the model is run from 16th at 14:00 to 20th at 02:00, spanning a total of 84 h (3.5 days), covering almost the entire passage of the typhoon in the SCS.

Storm surge simulations were validated against observations from six tide gauges ( Figure 1 ). The South China Sea Forecast Center, State Oceanic Administration granted hourly storm surge data recorded from 00:00 local standard time on the 17th to 10:00 local standard time on the 19th, July 2014, the 59 h from storm surge emergence to end. The entire observation period was covered in the model runtime. The capability of each model was evaluated only using those times when observations were available.

To simulate the storm surge forerunner caused by the remote wind effect at the open ocean, a large domain model grid was created as shown in Figure 1 . This model domain encompasses the entire SCS and a portion of the western Pacific Ocean. Resolutions of the model grid were set to range from 200 m near the shore to 20 km at oceanic open boundaries with 212,063 nonoverlapping triangular cells and 121,113 triangular nodes to enable superior flexibility in representing the bathymetry. For the 3-D model simulations, we employed 11 uniform terrain-following S-coordinates to properly represent the bottom boundary layer in the vertical. The nearshore bathymetry data were obtained from unpublished sea charts provided by the State Oceanic Administration of China in this study. For the open sea, the water depth data of the model grid were obtained from the ETOPO1 global relief dataset with a resolution of 1 min (1.5–2 km) from NOAA (http://www.ngdc.noaa.gov/mgg/global/global.html). These datasets were interpolated into model grid nodes. The model datum was referenced to mean sea level. Open boundaries of the model run without tidal and wave forcings in this study.

The study area (17°N–22°N, 108°E–113°E) is in the northwestern region of the SCS, including Hainan Island, and a part of the coast of Guangdong and Guangxi Provinces where storm surges frequently occur. There are six tide gauge stations located along the coastline. Haikou and Boao stations are located in Hainan Island on the left side of typhoon Rammasun track and the rest of the stations are on the western coast of Guangdong Province ( Figure 1 ).

The SLR data used in simulations is derived from the regional and local projections data provided by the IPCC 6th Assessment Report (AR6), which updates and expands the related assessments from previous reports, and provides improved projections and uncertainty estimates of future change (Fox-Kemper et al., 2021). In addition, the projections incorporate all potential processes that would affect sea level change with medium and low confidences. For the low confidence condition, the time series of the projection are up to the year 2300. The regional projections on a regular global grid are available from the NASA Sea Level Projection Tool (https://sealevel.nasa.gov/data_tools/17). In this study, we focus on four cases of future SLR by years 2050, 2100, 2200, and 2300 for the RCP 4.5 scenario ( Figure 2 ), which relates to medium greenhouse gas emissions. To reduce the uncertainties caused by multiple contributing factors of SLR, total integrated values of all components of SLR will be selected, which contain the assessed medium-confidence and low-confidence processes (total_ssp245_low_confidence_values). Furthermore, we used the 50% quantile (median) in the range of projections from the dataset in this study. The values of regional SLR are relative to the 1995–2014 baseline, close to the occurrence time of super typhoon Rammasun (2014). The values are suitable to be directly applied to the study without further modification.

The SCHISM model solves the standard Navier–Stokes equations with hydrostatic and Boussinesq approximations, and uses semi-implicit, finite-volume, and Eulerian–Lagrangian algorithms (Zhang et al., 2016). We executed the simulation by using the 2-D and 3-D barotropic modes of SCHISM. In 3-D mode, full 3-D cells were utilized in the entire domain with 11 vertical layers. The governing equations of the SCHISM model are composed of continuity equation and momentum equation. Continuity equation in 3-D (Eq. (1)) and 2-D depth-integrated forms (Eq. (2)) (note that bold texts in all equations represent vectors):

Momentum equation:

where vectors are expressed by using boldface type, ∇ = ( ∂ ∂ x , ∂ ∂ y ) is the horizontal gradient operator, u (x,y,z,t) is the horizontal velocity vector with Cartesian components (u,v),w is the vertical velocity, z is the vertical coordinate, n(x,y,t) is the free-surface elevation and represents storm surge height, h(x,y) is bathymetric depth, u is the depth-averaged velocity in the 2-D model, D/Dt is the material derivative, m _z is the vertical eddy viscosity term, and g is the acceleration of gravity. F represents the additional forcing terms, including horizontal viscosity, Coriolis, and atmospheric pressure. The baroclinic gradient, earth tidal potential, and radiation stress are neglected in this study. The vertical eddy viscosity term m_z denotes different forms in 3-D and 2-D as:

where ι is the vertical eddy viscosity and τ _w is the surface wind shear stress calculated based on the wind velocity vector U₁₀ on the 10-m height above the sea surface. τ _b is the bottom frictional stress. ρ is the water density and P_a is the air density. H=h+n is the total water depth. As boundary conditions, τ _w and τ _b are essential for the 2-D model and the 3-D model. At the sea surface, the internal Reynolds stress is balanced by the surface wind shear stress:

C_w is the drag coefficient, which is derived from the following formula (Pond and Pickard, 2013):

At the bottom boundary layer, the internal Reynolds stress is balanced by the bottom frictional stress:

The bottom stress is also computed by a similar quadratic drag law expressed as:

The u _b is the near-bottom velocity vector for the 3-D model. In this study, we calculated velocities by the 3-D model at 11 layers from the bottom of the water column to the surface, and the velocity at the first layer above the bottom represents the near-bottom velocity. However, the u _b of the 2-D model is different from the 3-D model and is the depth-averaged velocity vector. In the 3-D model, the C_d is generally the bottom drag coefficient computed according to the thickness of the logarithmic bottom boundary layer:

where K is the von Karman constant (κ = 0.4) and z ₀ is the bottom roughness parameter, which we set to 0.02, the typical value obtained at the oceanic boundary layer (Lueck and Lu, 1997). δ_b is the thickness of the bottom boundary layer that is closest to the bottom. In the 2-D model, C_d may be calculated by using a function of a Manning coefficient:

where n is Manning’s friction coefficient, which could be constant or variable in 2-D flow. The 2-D model is highly dependent on the Manning coefficient, and the previous studies concluded that the temporally and spatially varying Manning coefficient may significantly improve the accuracy of the 2-D model but could cause the storm surge model to be much more cumbersome and less robust (Lapetina and Sheng, 2015). Thus, the Manning coefficient is generally set as a constant of 0.02 in the 2-D model (Xu et al., 2010; Garzon and Ferreira, 2016). In this study, we used typical values for n and z₀ to compute the bottom drag coefficient under 2-D and 3-D formulations. There are significant differences of bottom drag coefficient between both formulations that may lead to different bottom stresses and have an impact on storm surge heights.

For assessment of model skill, accurately catching the peak values of storm surge is considered the criterion. Here, we calculated percentage errors (PE) of peak surge between model results and observations at six tide gauge stations to evaluate the capabilities of different wind forcings in hindcasting coastal storm surge. The formula was as follows:

where n _mod and n_obs are the modeled and observed peak surges, respectively. The mean absolute percentage error (MAPE) is commonly used as a measurement of simulation accuracy in model evaluation due to its very intuitive interpretation in terms of relative error. It is usually expressed as a ratio defined by the formula:

where n is the number of tide gauge stations (n = 6).

To compare time series of storm surges, model performance can be quantified with various skill metrics (Guillou, 2014; Akbar et al., 2017; Drost et al., 2017; Ding et al., 2020). The statistical errors of water surface elevations such as root-mean-square error (RMSE) and bias were frequently calculated for assessment of model capability. In this study, the metrics RMSE and bias were used to assess the absolute error between the modeled and observed values, with an ideal value of zero. In addition, Willmott skill (WS), which represented the degree of agreement, was also used for the quantitative evaluation of model performance (Willmott, 1984). A WS of 1 represents the best possible agreement between model simulation and observations, and complete disagreement was a value of 0. These metrics are calculated as follows:

where X _mod and X_obs are the modeled and observed values, respectively, and N is the number of observed samples. X mod ¯ and X o b s ¯ are the average values. These skill metrics are adopted to test the performance of the storm surge model.

In this study, we employed a parametric vortex model to generate sea surface pressure fields and wind velocity components. The real-time track data provide the hourly location of the TC center, the central pressure, and maximum wind speed at the height of 10 m. The wind velocity used in the model is composed of rotational and translational velocities. The rotational velocity is calculated from an analytical model of the radial profiles of sea surface pressures and winds, which is proposed by Holland (1980). Holland’s model function expression is as follows:

where r is the distance from the TC center,p is the local pressure at the location of r, P_c is the central pressure, P_n is the ambient pressure (1,010 hPa), R_m is the radius of maximum wind velocity, V(r) is the gradient wind velocity at r,P_a is the air density (1.15 kg/m³), e is Euler’s number, f is the Coriolis parameter (f=2ΩsinΦ), Ω is the angular frequency of the Earth rotation, Φ is the latitude of the TC center, and B is pressure profile parameter. B can be calculated from Eq. (18) by neglecting the Coriolis parameter and setting r=R_m as follows:

where V_m is the maximum gradient of wind velocity. This method was extensively used because the required maximum wind velocity and radial pressure gradient can be easily obtained from meteorological agencies (Hu et al., 2012). However, R_m is not provided with real-time track data. To supplement this, we used the best track data that were provided by the Joint Typhoon Warning Center (JTWC) (https://www.metoc.navy.mil/jtwc/jtwc.html). In their data, the radius of maximum wind speed had a 6-h temporal resolution, which was interpolated linearly to hourly resolution in this study. The direction of wind velocity components was adjusted by the inflow angle (α_F ), which is about 25° outside and decreases to 0 at the TC center (Mattocks and Forbes, 2008):

To remove the impact of the TC motion, a translational velocity vector (V_mov ) caused by the movement of the TC was added to the rotational velocity vector (V) for the stationary symmetric TC to create an asymmetric wind velocity vector (Ueno, 1981), and the translational velocity vector is calculated by the exponential function form:

in which V_mc is the moving velocity vector of the typhoon center. The wind field is calculated by the following formula:

where i and j are the unit vectors in the x and y directions, respectively; V_inital is the wind vector obtained from Holland’s parametric vortex model [Eq.(18)], and C ₁, C ₂ were, respectively, 1.0 and 0.8 in the parametric vortex model, which was confirmed by multiple comparative calculations (Yang and Yan, 2021).

According to previous research, the magnitude of wind simulated by the parametric vortex model decreased much faster away from the radius of maximum wind than the observed one and thus may underestimate the actual background wind (Lin and Chavas, 2012). In contrast, the reanalysis datasets are better at reconstructing the wind field in the areas farther away from the TC center but worse near the TC center (Pan et al., 2016). To improve the accuracy of the wind field, we used a combination of Holland’s parametric vortex model and the reanalysis data in this study. The method is considered to be an optimized combination of the two wind fields. We compared the combined data with single data to verify the conjecture that the combined data generally perform well.

The global ERA5 reanalysis dataset provided hourly data with a spatial resolution of 0.25°(latitude) × 0.25°(longitude) on single levels by the European Centre for Medium-Range Weather Forecasts (ECMWF) (https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels). A smooth transition is very important between the wind field of the parametric vortex model and the wind field obtained from the ERA5 reanalysis data (hereinafter abbreviated to ERA5 data). To avoid discontinuity between the two wind fields, we used a weighted average formula to smooth the transition between the two wind fields:

where V_c is the combined wind vector, V_b is the background wind vector from the ERA5 dataset, C is an intermediate variable, α is the weight coefficient, and n=4 (Li et al., 2020).

Figure 3 presents the horizontal wind speed profiles and spatial distributions of typhoon Rammasun when it landed on Hainan Island. The main differences between the parametric wind profile and the ERA5 wind profile were reflected in the differences in the peak wind. The Holland parametric TC model can reproduce the maximum wind speed (V_m ) with >50 m/s, which was close to the observed value. However, wind speeds of ERA5 data are lower than 30 m/s and the radius of maximum wind is much larger than the radius provided by the JTWC track data. Thus, this reconstruction method was used to improve the simulated results by combining the better part of the two wind profiles.

In order to compare the performance of the SCHISM model under three different forcings, three tests (Tests 1–3) were carried out. In Test 1, the ERA5 pressure and Holland’s forcing wind fields were applied. In Test 2, we used the ERA5 data to execute the storm surge simulation. In Test 3, the reconstructed wind fields (see in Eqs. (23)–(25)) and the ERA5 pressure were assembled to drive the model. In this section, we ran 2-D simulations and compared storm surge hindcasts in three tests in the study area.

Figure 4 shows comparisons of observed and simulated storm surges at six tide gauge stations, and their PEs are given in Table 1 . It can be seen that Test 3 gave a more accurate storm surge hindcast than the others. The simulation in Test 1 underestimated the peaks of storm surge and poorly modeled the surge forerunner caused by the remote wind at most stations. The result of Test 2 fit the observations better for the surge forerunner and its peak surges were also closer to the observed peak values than Test 1, except the Haikou station. The low-resolution reanalysis data (0.25°×0.25°) were unable to accurately reproduce wind speed at the cyclone center that led to the largest PE (−58.04%) at the Haikou station during Test 2. Test 3 had the lowest PEs in magnitude at all stations with the peak surges best matching observations; the combined wind field retained the advantages of both wind fields in Tests 1 and 2. Overall, the MAPEs ( Table 1 ) show that the best performance is that of Test 3 (7.57%), the second is Test 2 (24.85%), and the worst is Test 1 (31.93%).

To further evaluate model skill, statistical metrics (RMSE, Bias, and WS) about the time series of simulated and observed storm surge at all tide gauge stations were given ( Figure 5 ). The three statistics show that the result of Test 3 performed best among the three tests. The averaged RMSE (15.9 cm) and Bias (−9.3 cm) values of Test 3 are smaller than results in Test 1 (41.9 cm and −33.8 cm) and Test 2 (17.9 cm and −10.3 cm). The averaged WS value (0.95) in Test 3 also indicated that it is more consistent with observation than Test 1 (0.71) and Test 2 (0.93). Overall, hindcast of storm surge under Test 3 shows the best agreement with observations among three tests. Thus, in the subsequent comparison of 3-D and 2-D models, we selected the ERA5 pressure and the combined wind field as the atmospheric forcings for the model.

In this section, the above statistical parameters were used to compare the performance of the 2-D (Test 4) and 3-D (Test 5) models. Figure 6 shows that although the simulation results of the 2-D and 3-D models are both in good agreement with the observed surges during typhoon Rammasun, there are some differences expected from the setup of both models. The time series of surges of the 2-D model at all stations were smaller than those of the 3-D model’s during the period. To clarify that, the PE and MAPE of the two tests were calculated to quantify the model reliability ( Table 2 ). The 3-D model (4.22% MAPE) performed slightly better than the 2-D model (7.57% MAPE) overall. The PE values in Table 2 indicate that the 2-D model tended to underestimate the storm surges at most stations due to overestimation of bottom stress in the 2-D model (Weisberg and Zheng, 2008). The performances of the 3-D model at Zhapo, Zhanjiang, Naozhou, and Haikou stations (0.23%, 0.66%, −5.22%, and 0.13%) were better than the 2-D model (−5.14%, −3.65%, −10.49%, and −4.80%), while they were close at the remaining two stations. Overall, compared to the 2-D model, the peak predictions of the 3-D model had better quality, and their discrepancies between observation and simulation were considered acceptable.

To further evaluate the performance of the two models, RMSE, Bias, and the WS values of both models were calculated and plotted in Figure 7 . The average RMSE (13.4 cm) and average Bias (−3.8 cm) values of the 3-D model are smaller than results from the 2-D model (15.9 cm and −9.3 cm). Individually, the 3-D model has slightly lower accuracy than the 2-D model at Boao station. The WS values also demonstrated that the 3-D model (0.96) performed slightly better than the 2-D model (0.95) on average. From the above comparison, we can find that the 3-D SCHISM model reproduced the storm surges better than the 2-D model in the study area.

The spatial discrepancies between 2-D and 3-D simulated peak surges are described in this section. We firstly plotted the spatial distribution of the peak surges under the 2-D condition ( Figure 8A ). The high storm surge mainly occurred along the eastern Leizhou Peninsula coast and offshore of the northeastern Hainan Island, with heights of >2 m. Offshore of the northeastern Hainan Island experienced especially severe water levels of about 3-m height due to being close to the track of typhoon Rammasun. In this region, it can be found that the 3-D model generates higher surges than the 2-D model, and the largest discrepancies between 2-D and 3-D simulated results were over 0.2 m ( Figure 8B ). Since the 2-D model generally underestimates the storm surges, the 3-D model is suggested as a better choice for surge simulation in the study area.

Four tests with different TC’s translation velocities are conducted in this study, named Tests 6–9 with 4, 2, 1/2, and 1/4 times initial moving velocity (V_mc ), respectively. Figure 9 plots the discrepancies of peak surge heights between four tests and Test 5 (base case hereafter). When translation velocity was doubled, the changes of peak surge in most areas were insignificant except for the coastline of the western Guangdong Province in which peak surges increased from 0.1 to 0.5 m ( Figure 9B ). However, peak surges in most areas were lower than the initial scenario (extreme scenarios yield discrepancies of −1.6 m) when TC’s translation velocity was quadruple ( Figure 9A ). Halving the translation velocity produced lower peak surges (below −0.2 m) along the eastern Leizhou Peninsula coast ( Figure 9C ). Reducing translation velocity to a fourth, peak surges were slightly lower than the scenario with the 1/2 translation velocity ( Figure 9D ). Increased peak surges along the Guangxi Province may be caused by the terrain where sea water easily piled up in semienclosed bays. As the duration of the typhoon was longer, the phenomenon became more significant. For the region of complex terrain, due to the propagation speed of long waves below TC’s translation velocity, the typhoon wind cannot effectively affect the sea water redistribution to weaken storm surge heights. Conversely, when the typhoon has the lowest translation speed, higher surges occur in the semienclosed bay.

In most regions, the peak surge at Test 7 (doubled translation velocity with 30 mph) is the highest, and if we continue to increase or decrease translation velocity, the peak surge will decrease. When considering the spatial distribution of peak surges, the TC’s translation velocity moving across the shoreline is more likely to be the main factor for affecting surge heights. The relatively faster TCs within limits (V_mc ≤13.4 m/s) are more dangerous and cause higher surges.