Investigation of Storm Tides Induced by Super Typhoon in Macro-Tidal Hangzhou Bay

Typhoon-induced storm tides can cause serious coastal disasters and considerable economic losses. Understanding the mechanisms controlling storm surges helps the prevention of coastal disasters. Hangzhou Bay (HZB), a typical macro-tidal estuary, is located on the east coast of China, where typhoons frequently occur. The funnel-shaped topography makes this macro-tidal bay even more sensitive to storm tides. Super Typhoon Chan-hom was used as an example to study the characteristics and dynamic mechanisms of storm surges using a well-validated numerical model. The model considers the two-way coupling of waves and tides. The wind strength for the model was reconstructed using multi-source wind data and was refined by considering different rotating and moving wind fields. The Holland–Miyazaki model was used to reconstruct the local wind-field data with a good performance. The model results show that the total water level of HZB during typhoon Chan-hom was mainly dominated by tides, and the storm surge was closely related to the wind field. Surface flow was mostly influenced by winds, followed by tides. The spatial and temporal distributions of the significant wave height were controlled by the wind and local terrain. Wind stress was the largest contributor to storm surges (91%), followed by the pressure effect (15%) and the wave effect (5%). Both wind and wave-induced surges occurred during low slack waters. The tide-surge interaction changes (enhance or suppress) the surge by approximately 0.5 m during the typhoon, comprising approximately 50% of the total surge. Tides interacted with surges through various mechanisms, from the bay mouth (local acceleration and friction) to the bay head (friction and advection). The Coriolis force had a relatively minor effect. The findings of this study provide useful information for studies on sediment dynamics and coastal structures under extreme weather conditions.


INTRODUCTION
Storm surges are an abnormal rise in sea level caused by strong atmospheric turbulence, such as strong winds or sudden changes in atmospheric pressure (Feng, 1998). Storm tides combined with heavy rain induce coastal inundation, especially in spring in high slack waters. Storm waves are the main force involved in coastal and seabed erosion and thus cause damage to coastal structures. Storm events can persistently alter the morphological evolution of intertidal flats. The magnitudes of some bed-level changes are comparable to years of continuous evolution (de Vet et al., 2020). Therefore, understanding the mechanisms of storm surges during extreme weather conditions is of great importance.
The mechanism of storm surges has been investigated extensively. The factor affecting them can be divided into two categories: internal and external. Internal factors include pressure, wind stress, rainfall, waves, and the interaction between astronomical tides and typhoon-induced storm surges (Zhang and Sheng, 2015;Chu et al., 2019;Wang and Sheng, 2020;Shankar and Behera, 2021). External factors include the Coriolis force, pressure gradient, bottom stress, and sea-level rise Yang et al., 2019;Marsooli et al., 2019;Thuy et al., 2020). In numerical simulations, attention has been paid to wind field reconstruction, as wind and pressure are key factors affecting storm surge simulation. Pressure and wind fields are the direct driving forces of storm surges. During typhoons, highaccuracy pressure and wind fields are key factors for reproducing the evolution of storm surges (Lin and Chavas, 2012;Torres et al., 2019). Reanalysis of wind field data based on numerical models and data assimilation techniques by the European Centre for Medium-Range Weather Forecasts (ECMWF) and the National Centers for Environmental Prediction (NCEP) are open sources and are widely used by researchers to study ocean and offshore storm surges (Brenner et al., 2007;Lv et al., 2014). However, the wind speed near the cyclone center of this reanalyzed data is usually lower than the measured value, which cannot accurately depict typhoon structures (Signell et al., 2005;Cavaleri and Sclavo, 2006). Therefore, empirical models of the pressure and wind fields require optimization.
Considerable attention has been paid to the characteristics and nonlinear effects of storm surges. As early as 1974, Darbyshire discovered that the interaction of astronomical tides and storm surges is nonlinear, and that this nonlinear interaction can form additional storm surges (Darbyshire, 1974). The time of typhoon landfall, tidal cycle, water depth, air pressure, typhoon track, and geographic location affect the magnitude of the nonlinear effect (Rego and Li, 2010;Zhang et al., 2019). The nonlinear effects of tides and storm surges mainly originate from the friction term, shallow water effect, and advection term. Many scholars (Wolf, 1978;Jones and Davies, 2007) have conducted relevant research to determine which of the three is dominant. The results show that the combined effect of wind stress, bottom friction term, and advection term are the main factors, whereas the term related to local acceleration and Coriolis force affects the nonlinear level (Jones & Davies, 2007;Yang et al., 2019). The nonlinear effect of tide-surge interaction enhanced/reduced the surge during low/high slack water (Rego and Li, 2010b). The closer the bay area is to the interior of the bay, the stronger the nonlinear effect (Zhang et al., 2010).
An important nonlinear process during typhoons is the tidesurge interaction, especially in macrotidal estuaries. There are two main physical mechanisms involved in the interaction between tides and storm surges. The first is the generation of storm surges that shift the tidal phases (Olbert et al., 2013;Zhang et al., 2017). Second, tides can affect the storm surge magnitude. Changes in water depth lead to changes in storm surges, which are larger in shallow waters than in deep waters (Horsburgh and Wilson, 2007). He et al. (2020) discussed the contributions of wind, pressure, and waves to storm surges, in which wind stress played a major role. Zhang et al. (2019) applied the ADCIRC model to simulate South China Sea storm surges caused by typhoons along different paths. The results showed that tidal phase modulation was mainly affected by tidal phase interaction, and that the surge peak time was not sensitive to tidal phase changes.
This study uses the super typhoon Chan-hom as an example to study the characteristics and mechanisms of storm tides in macro-tidal Hangzhou Bay. First, a numerical model considering tides, waves, and surges was developed and calibrated using field data. The wind field in the model was reconstructed using an optimized wind model. The characteristics and mechanisms of storm tides were then examined using model results and numerical tests. The numerical model is described in Section 2 and validated in Section 3. The model results and discussion are provided in Sections 4 and 5, respectively.

Tide-Wave Coupled Model
The hydrodynamic model, finite volume coastal ocean model (FVCOM) is an unstructured ocean model using a finite volume method to discretize three-dimensional original equations (Chen et al., 2003). To fit the complex shoreline better, the model adopts an unstructured triangular mesh in the horizontal direction. In the vertical direction, the model adopts s coordinates. The continuity and momentum equations considering the wave radiation stress in the s coordinate system are as follows: where x, y, and s are the north and vertical coordinates of the middle eastern area in the s coordinate system, respectively. t is time, and u, v, and w are the velocity components in the x, y, and z directions, respectively. D is the total water depth, D = h + z; f is the Coriolis force parameter, r is the density of seawater, P atm is the atmospheric pressure, g is the acceleration due to gravity, K m is the vertical eddy viscosity coefficient, and K h denotes the vertical thermal eddy friction coefficient. t x , t y and R x , R y denote the turbulent and wave radiation stress terms, respectively. The items on the left-hand side of the formula are the local and horizontal convection terms in x, y, and z directions and the Coriolis force term. The right-hand side of the formula is the hydrostatic atmospheric pressure strength turbulent stress horizontal diffusion radiation stress, which is a static pressure model that does not consider hydrodynamic pressure. Thermohaline is assumed as a constant as there is no in-depth research on thermohaline distribution, so relevant formulas are not provided.
At the sea surface level s=0: At the bottom of the sea s=-1 Where, the subscripts x and y indicate the x and y directions, respectively. t sx and t sy are the components of the surface wind stress; t bx and t bx are the components of the base stress; u s and v s are the components of the surface wind speed; and u b and v b are the components of the bottom velocity. k is the Karman constant, taking 0.4; near the sea surface, Z ab and Z 0 are the sea surface roughness values. Near the seabed, Z ab is the distance between the nearest grid and the seabed, and Z 0 is the seabed roughness. C s is wind drag coefficient. C d is the bed friction coefficient. FVCOM-SWAVE is based on SWAN (the third generation wave model), and retains the physical process of SWAN in the grid division. It adopts unstructured triangular grid, and selects the finite volume method in the discrete method. Waves are described using the wave action density spectrum, which is defined as the energy density divided by the natural frequency: The conservation equation of wave action is described in Cartesian coordinates as: Where s, q, x, y, t are the coordinates of frequency wave direction in space and time. N is the wave action density. The first term on the left side of the equation represents the local variation of wave action density in time. The second and third terms represent the propagation of wave action density in space. The fourth term represents the change of the wave action density in the frequency shift s space due to changes in water depth and flow. The fifth term shows the variation of wave action density in q space by refraction and shallower action caused by water depth and flow. S tot is the energy source term on the right side of the equation, which represents the generation, dissipation and redistribution of wave energy. The ocean dynamic model and the wave model use the same horizontal unstructured grid variable transfer diagram as shown in Figure 1F. Among them, the wave model provides parameters such as significant wave height, wave direction, wavelength, spectral peak period, bottom orbital speed, bottom orbital period and other parameters to the ocean model. The ocean model provides the depth-weighted average Doppler flow rate and water level to the wave model. The bottom boundary layer model obtains the flow velocity from the hydrodynamic model, the wave elements from the wave model, and the bottom bed roughness and other variables from the bottom boundary layer model. After calculation, the bottom stress is fed back to the ocean model. Current-wave coupling is approached through radiation stress, bottom boundary layer, surface stress, and morphology. The radiation stresses are added into the momentum equations of FVCOM to include the wave-driven motions (FVCOM user manual).

Wind Field Re-Construction Model
Pressure and wind field models are essential for accurate storm surge predictions. Parameterized wind field models are widely accepted by oceanographers owing to their high efficiency in calculation (Zhong, 2019). In classic atmospheric models, the wind field comprises a rotating and moving wind field. The former is due to the balance of the pressure gradient, Coriolis, and centrifugal forces, and friction, whereas the latter is mainly related to the movement of the cyclone center (Pan et al., 2016).
The maximum wind speed radius (Table 1), R max , is the distance between the typhoon center and the maximum wind speed. This is the most important parameter in typhoon models (Chang et al., 2015). The method of Graham and Nunn (1959) (GN59) is a commonly used empirical formula for calculating the maximum wind speed radius in the early stage of typhoon research. Based on a large amount of field data, Knaff et al. (2018) applied the empirical formula of R max to the western Pacific Ocean, Eastern Pacific Ocean, and Western Atlantic Ocean.   Li (2007) studied the correlation between typhoon parameters in the southeast coastal areas of China using the typical correlation analysis method and believed that there was a high correlation between R max and the central pressure difference. Jiang et al. (2008) proposed an empirical formula for R max in the form of a power function. This formula is applicable to the northwest Pacific Ocean based on a previous statistical empirical formula and the measured pressure data of the typhoon center. Typhoon Chan-hom was used as an example to verify the accuracy of the above methods for R max . The temporal variation in R maxv , calculated using the above equations, is shown in Figure 2A. The data of geophysical locations, pressure, and moving speed of the typhoon center were obtained from tcdata.typhoon.org.cn (Ying et al., 2014;Lu et al., 2021). According to Zhang (2015), the stronger the typhoon, the smaller the R max value. R max calculated according to the method of Knaff (2007), Knaff et al. (2018) has only small temporal variations. The R max calculated by the GN59 method (Graham and Nunn, 1959) was small, with a smaller R max occurring when the typhoon was first generated (< 30 h) than that at the peak of the typhoon (75-90 h). These characteristics are not consistent with the physical law of typhoons, indicating that the methods of Graham and Nunn (1959) and Knaff (2007), Knaff et al. (2018) were not suitable for this study. F is the latitude of the cyclone center, P c is the central pressure of the cyclone (hPa), V c is the moving speed of the cyclone center. DP is the pressure difference between the peripheral pressure and cyclone center (hPa).P min is the lowest central pressure measured off the coast of China, 895 hPa; R min is the mean of the maximum wind speed radius corresponding to Pmn, which is 27.9 km. Knaff et al. (2007) R max = 70.376-0.1124V c -0.0074 (f-25) Li (2007) R max =exp(-0.0511*DP 0.7515 +4.9213) Jiang et al. (2008) R max = 1119*(1010-P c ) -0805 Zhou (2018) R max = 29.178exp[0.0158(P 0 -900)] Zhong The pressure field and rotating wind field Author (year)

Formula Annotations
Takahashi, (1939) P(r) = P ∞ − P ∞ − P c 1 + r=R max P ∞ is the pressure outside the typhoon (1013 hPa), r is the distance between the calculated point and typhoon center, and V max is the maximum gradient wind speed. f is the Coriolis parameter, r a is the air density (r a = 1.205 kg/m 3 at 20°C). Fujita (1952) P Moving wind field Author (year) Formula Annotations Miyazaki et al. (1962) W mov is the speed of the moving wind field, V c is the moving speed of the cyclone center. R G is the empirical attenuation coefficient, which is 500 km or 10Rmax. Ueno (1981) W To compare the accuracy of the R max methods, we used typhoon data from the west Pacific in 2013-2015, which were published by the Joint Typhoon Warning Center (JTWC), to perform validations ( Figure 2B). The R max values calculated by Zhou (2018) were relatively large. The R max calculated by Jiang et al. (2008) model exhibited the best performance when the center pressures were low. However, when the typhoon intensity decayed, the pressure differences decreased, and the deviation of R max calculated using the method of Jiang et al. (2008) increased quickly. R max calculated by the methods of Li (2007) (SS=0.737) and Zhong (2019) (SS=0.682) showed better performance. As R max calculated according to the method of Li (2007) has a higher SS value, we used this method to calculate R max in this study.
Compared with the moving wind field, the rotating wind field contributed more to the total wind field. It is important to choose a suitable pressure field model for the calculation of the rotating wind field (Table 1). At present, the commonly used air pressure wind field models in China include the Takahashi (Takahashi, 1939), Fujita (Fujita, 1952), Jelesnianski (1966Jelesnianski, 1965), Fujida-Takahashi (Wang et al., 1991), and Holland models (Holland, 1980). The Jelesnianski pressure model was proposed by Jelesnianski (1965), and is the pressure field model used by the sea, lake, and overland surges from hurricanes (SLASH). Introducing vorticity theory, Jelesnianski (1965) proposed and improved an empirical formula for the wind field. Through many comparative analyses,  Wang et al. (1991) used the Fujita formula (Fujita, 1952) to calculate the distribution of atmospheric pressure in the areas within the distance of 2R max , and used the Takahashi formula (Takahashi, 1939) to reflect the change of atmospheric pressure in the areas exceed the distance of 2R max . When air particles move in the atmosphere, the gradient wind formula can be obtained by solving the momentum equation and considering the balance of the pressure gradient force, Coriolis force, and centrifugal force and ignoring the sea surface roughness. Holland (1980) introduced the shape coefficient B and proposed a pressure field model with more general applicability.
The above pressure field model was dimensionless, as shown in Figure 2C, where the shape coefficient B=1.2 was taken from the Holland model. Within the maximum wind speed radius R max , the growth rate of the Takahashi model was steep and the calculated value was high. The Jelesnianski model had a relatively stable growth rate. The simulation results of the Fujita and Holland models are relatively reasonable. Outside R max , with an increase in the distance from the typhoon center, the pressure values of all models gradually increased, with a gradually decreasing changing rate, and finally stabilized at P ∞ . The Holland model had a relatively large growth rate outside R max , and the calculated value of the external pressure outside 4R max was the largest among the five models.
Taking the typhoon Chan-hom as an example, the characteristics of the different rotating wind field models were examined. The Jelesnianski model simulates low wind speeds within R max , and relatively fast attenuation outside R max . The maximum wind speed calculated by the Takahashi model was smaller than those of the other models. The representative wind field distribution of the Takahashi model increased with increasing distance from the typhoon center. The maximum wind speed of the Fujita model was closest to reality. The Fujida-Takahashi model absorbs the advantages of both the Takahashi and Fujida models, but the calculated wind speed changes at 2R max . The wind speed of the Holland model near the typhoon center was smaller than that of the Fujita-Takahashi model, and the Holland wind field outside 2R max was the largest among the five models.
If the reconstructed wind field contains only the rotating wind field ( Figure 2D), then the wind field usually has a spatially symmetric distribution. In the wind field in the northern hemisphere, owing to the movement of typhoons, the rotating wind on the right-hand side and the moving wind superimpose on each other, while canceling each other on the left-hand side. Therefore, the typhoon wind field had a spatially asymmetric distribution. Commonly used models for the calculation of moving wind fields (Table 1) are the Jelesnianski model (Jelesnianski, 1966), Miyazaki model (Miyazaki et al., 1962), and Takeo Ueno model (Ueno, 1981). The moving wind field calculated by the Miyazaki Masawei model (1962) was affected by the empirical attenuation coefficient, and the wind speed decreased with the distance from the typhoon center. The Ueno model (1981) used the maximum wind speed radius to calculate the moving wind field of typhoons.
The moving wind field model above is dimensionless, as shown in Figure 2E. Both the Jelesnianski and Ueno models obtained the maximum wind speed at the maximum wind radius. In the Jelesnianski moving wind field model, the wind speed at the center of the typhoon was zero, R max reached 0.5V c , and then the wind speed decreased slowly. In Ueno's model, the central moving wind speed is approximately 0.456V c and reaches the maximum value V c at R max . Then the wind speed decreases rapidly, and the shifting wind speed was less than 0.1V c at about 4R max . The maximum moving wind speed V c calculated using the Miyazaki model was taken at the center of the typhoon, and the moving wind speed gradually decreased with increasing distance from the typhoon center. The calculated value of V c using the Miyazaki model at about 5.6R max was weaker than that calculated using the Jelesnianski model.
In the Jelesnianski moving wind field model, the wind speed at the center of the typhoon was zero, R max reached 0.5V c , and then the wind speed decreased slowly. In Ueno's model, the central moving speed was approximately 0.456 V c and reached the maximum value V c at R max . Then the wind speed decreased rapidly, and the shifting wind speed was less than 0.1 V c at about 4R max . The maximum transit wind speed V c of the Miyazaki Showei model was measured at the center of the typhoon, and the transit wind speed gradually decreased with increasing distance. The calculated value using the Miyazaki model at about 5.6R max was weaker than that calculated by the Jelesnianski model.
To obtain a more accurate wind field, this study combined the rotation of the wind field, the enhancement, and the NCEP wind field (rda.ucar.edu/datasets/ds084.1). The wind field near the typhoon center mainly depended on the calculated values of the model, whereas the NCEP wind field was applied to a large ocean area. The calculation formula is as follows (Wei et al., 2019): Where, W mov_x and W mov_y are the x and y components of the moving wind field, and W NCEP_x and W NCEP_y are the x and y components of the NCEP wind field, respectively. q represents the angle between the line connecting the calculation point with the typhoon center and the east (counter-clockwise). b in is the inflow angle, representing the deflection angle of the gradient wind vector across the isobars, which is approximately 20° ( Wang et al., 1991); and C 1 is the correction coefficient.
According to the theory of the atmospheric boundary layer, it is 0.71 (Wei et al., 2017). n = 9 is a constant.
We selected Shipu (121.95°, 29.20°N, 127 m), Dinghai (122.12°, 30.03°, 37 m), Baoshan (121.47°, 31.40°, 4 m), and Qidong (121.60, 32.07°, 10 m) as the field stations to illustrate the reconstructed wind data. The measured wind field data at 10 m was compared with that calculated by different structural wind field models during the influence of typhoon Chan-hom. The data was obtained from the Real-time Database of Atmospheric Environment (envf.ust.hk/dataview) of the Institute of Environmental Science, Hong Kong University of Science and Technology. Near the sea surface, the vertical profile of the wind speed was close to the logarithmic law. The vertical profile of the wind speed was close to the power law from 100 m above the sea surface to the top of the friction layer: Where, W z1 and W z2 are the wind speeds at altitudes z 1 and z 2 , respectively; z 0 is the surface roughness, 0.03 m on the ground and 0.003m on the sea; and P is the empirical coefficient, taking 0.14.
The SS (skill scores) for the different wind models are listed in Table 2. The definition and calculation method of SS are shown in Eq. 17. The rotating wind data dominated the characteristics of the wind field, whereas the moving wind data only slightly adjusted the spatial and temporal distributions of the wind field. The rotating wind data calculated using the Jelesnianski model matched the moving wind data. The rotating wind data calculated using the Fujita-Takahashi model matched the moving wind data of Ueno better. The center wind speed calculated by the Holland method was smaller than that of the field data. The moving wind data calculated using the Miyazaki model enhances the center wind speed calculated using the Holland model. Hence, during typhoon Chan-hom, the Holland-Miyazaki model performed better and was used to reconstruct the wind field data in this study. Wind data reconstructed using the Holland-Miyazaki model are shown in Figures 2F, G.

Super Typhoon Chan-hom
The super typhoon Chan-hom (International Code: 1509, Joint Typhoon Warning Center: 09 W) was the ninth storm of the 2015 Pacific typhoon season ( Figure 1B). At 12:00 on June 30, 2015 (CST, which is China Standard Time), Chan-hom developed over the northwest Pacific Ocean and then moved to the northwest, becoming stronger and reaching supertyphoon level. Chan-hom peaked on July 10 and then entered the East China Sea. At 09:00 on the 11th day, Chan-hom downgraded to a strong typhoon level and landed on the coast of the East China Sea. After landing, Chan-hom moved northward and eastward, and gradually lost its strength.

Model Domain and Configurations
The model domain ranges from 120°E in the west, 125.5°E in the east, 27.5°N in the south, and 34°N in the north, covering the entire Hangzhou Bay, Changjiang River Estuary, and the Zhoushan Islands. To improve computational efficiency and accuracy, the resolution of the model grid was set to 30 km at the open boundary of the ocean. Inside the bay, the resolution of the grid in the coastal area was refined to 200 m and 100 m because of the complex shoreline and island structure ( Figures 1A, C).
Shoreline data was obtained from high-precision data provided by the National Geophysical Data Center (NGDC) (https://www. ngdc.noaa.gov/mgg/shorelines/gshhs.html). Bathymetric data was obtained from ETOP1, a global model developed by the National Geophysical Data Center (NGDC). Open-boundary tide data was derived from the global ocean tide prediction model TPXO developed by Oregon State University. The Yangtze River discharge data were obtained from the Datong station of the water level management system (http://yu-zhu.vicp.net/), and sediment discharge data were obtained from the Yangtze River Sediment Bulletin. The flow data of the Qiantang River were obtained from the flow data of the Fuchun River hydrology station (Runoff temperature constant 20°, salinity constant 0‰). Wind and pressure field data under calm weather conditions are derived from reanalysis data provided by the ECMWF (http://apps. ecmwf.int/datasets/data/interim-full-daily/levtype=sfc/). Wind speed data at 10 m above the sea surface was selected as wind forcing data, and sea surface pressure was selected for the air pressure data. The temporal resolution of the two datasets was 6 h and the spatial resolution was 0.125 degree. In addition to the ECMWF data, the wind field and pressure field data during typhoons also included the latitude and longitude, the minimum pressure, and the movement speed of the typhoon center from tcdata.typhoon.org.cn.

Numerical Tests
Eight numerical experiments were designed based on the tidewave coupling numerical model ( Table 3). Case 0 is the control run. Case 1 removes the wind field based on the control condition (Case 0). Case 2 removes the influence of air pressure based on Case 0. Case 3 removes the wave action There is a nonlinear relationship between the increase in wind stress and astronomical tidal level. To explore the relationship between winds and astronomical tidal level, a numerical simulation (Case 5) was designed. In Case 5, only wind forcing was considered and tidal forcing was removed. To explore the relationship between wave-induced surges and astronomical tide levels, numerical simulation Case 6 was designed. In Case 6, only wave was considered and tide was removed. Numerical simulation Case 7 was set up to further investigate the nonlinear effect of the wind field and astronomical tide. Case 7 considered wind, wave and air pressure, but ignored tidal forcing.

Field Data
The distribution of field stations is shown in Figure 1B. The measured data were obtained from the Shanghai Meteorological Administration Numerical Prediction Innovation Center and Zhejiang Ocean Monitoring and Prediction Center. The measured tidal level data are at five stations, Z1, Z2, Z3, Z4, and Z5, as shown in Figure 1B. The tidal flow verification only includes H3/W4 in Figure 1B. The wave measured data were obtained from five stations, H1, H2, H3, H4, and H5, from July 1 to July 15, 2015 (CST, Figure 1B). Root mean squared error (RMSE), correlation coefficient (CC) and skill scores (SS) were used to evaluate the reliability and accuracy of the model (Murphy, 1992): where m i and O i are the model simulation and measured data, respectively. m and O are the averages of the simulated and measured values, respectively. S m and S o are the standard deviations of the simulated and measured values. The closer the correlation coefficient (CC) is to one, the greater the correlation between the simulated and measured values. When the skill scores (SS) is greater than 0.2, the model has some reliability. When the skill scores (SS) is greater than 0.5, the model is highly reliable.

Model Validation
Tidal levels at the five sites Z1, Z2, Z3, Z4, and Z5 (red dot in Figure 1B). The data period was from July 1 to July 15, 2015 (CST). The CC and SS of the statistical error between the simulated and observed tidal levels at each station is shown in Table 4. The CC and SS of tidal levels at each station is between 0.96 and 0.98 ( Table 4). The model simulates the changes better in the tidal level during typhoons ( Figure 3A).
The tidal flow verification data of the H4 station are from July 1 to July 15, 2015 (CST, Figure 1B). The flow speed validation SS is 0.88 and the CC is 0.81. The flow direction validation SS is 0.82 and the CC is 0.67 ( Table 4). The tidal current verification is within an acceptable accuracy range ( Figure 3B). More information on model validation was provided by Ye (2019) and Ren (2022).
Significant wave height data were obtained from stations H1 to H5 (green points in Figure 1B) at the verification site. The data period was from July 1 to July 15, 2015 (CST). The SS for the significant wave height at each site was higher than 0.95, and the minimum CC was 0.93 ( Table 4). The simulated wave height was in good agreement with the measured value ( Figure 3C). The model can accurately simulate the processes of wave generation and extinction during typhoons. Datong, Fuchunjiang M 2 , S 2 , K 1 , O 1 , N 2 , P 1 , Q 1 , K 2 Initial conditions water level, current and wave is set to zero; temperature is set to 20°C; salinity is set to 30‰ Wind data sources U and V components provided by ECMWF. The time and spatial resolution is 6 h and 0.125°respectively. Air pressure data sources NCEP-FNL pressure data. The time and spatial resolution is 6 h and 1°respectively.

Tides and Storm Surges
Three sections (C1, C2, and C3) and six points (S1, S2, S3, S4, S5, and S6) were used to illustrate the characteristics of the tides and surges in the bay. The water depth of the six feature points is 9.0 m, 14.7 m, 6.5 m, 12.9 m, 8.0 m, and 11.0 m, respectively. The surges during the storm tides is represented by the differences between the sea surface level in Case 0 (the control run, Table 5) and the sea surface level in calm weather (Case 4, Table 5).
Typhoon Chan-hom passed through Hangzhou Bay and its nearby areas from July 10 -11, 2015 (CST, Figures 4A, B). As the typhoon gradually approached, the sea surface level at the mouth of the bay was initially affected. It began to increase, and the maximum sea surface level at section C3 at the bay mouth reached 2.5 m. At 09:00 on July 11, the sea surface level at section C1 at the bay head reached 4.0 m, and then the maximum sea surface level areas continued to advance upstream, and finally reached a peak of 5.0 m at Yanguan. Similarly, the area of peak increase in sea surface level also began to advance from the bay mouth to the bay head. The maximum sea surface level at section C3 at the mouth of the bay was about 1.2 m, and the maximum sea surface level at section C2 at the middle of the bay was 1.5 m. Finally, at 09:00 on July 11, the maximum increase value of sea surface level (3.0 m) was reached at the head of the bay, and the variation trend of the two was basically the same.
In the bay mouth, the surge always remained small, the bay mouth was relatively wide, and the backwater effect caused by winds and waves was small. The surge increases earlier than the sea surface level, and the maximum surge moment occurs earlier than the moment of the high slack water level.
The maximum sea surface level and surge occurred when Hangzhou Bay was at the radius of the maximum wind speed ( Figures 4C, D). The wind speed in the bay was relatively uniform. Due to the influence of the funnel-shaped topography of the bay, the maximum water level and maximum water surge both occur at the narrow bay head and the maximum water level is up to 5.0 m in the upstream Yanguan station (YG in Figure 1B). The surge at the bay head was generally more than 2.0 m, and the surge at Yanguan was the largest, up to 3.0 m.
In counter-clockwise cyclones, sea surface water was driven by wind from the northern areas to the southern area and blocked by the southern coastlines. Larger sea surface level and surge occurred near the southern bank. Spatially asymmetric distributions of the sea surface level and surge occurred owing to the impacts of typhoons and geomorphology. The sea surface level in the bay was approximately 1 m and 0.5 m higher than in the open ocean and bay mouth, respectively.
The surge is mainly related to the wind field. Figure 4E shows the time series of the total water level and surge at each feature point. The maximum water level at the six feature points is 4.03 m, 4.06 m, 2.93 m, 2.96 m, 2.21 m, and 2.44 m, respectively ( Figure 4E). In all three sections, the water level on the southern bank (points S2, S4, S6) is lower than that on the northern bank (points S1, S3, S5). The maximum water level of the six feature points occurred in high slack water, and the difference between them was small.
The maximum surge of the six feature points is 1.95 m, 2.04 m, 1.67 m, 1.74 m 1.12 m and 1.50 m, respectively. The peak surge at the bay mouth was less than that at the head of the bay, and there was no difference between the peak surge in the southern and northern areas. The maximum surge time of the six feature points was gradually delayed from the mouth to the bay head. With the movement of typhoons, the peak surge rapidly advances into the bay. With the advancement of typhoons, the extreme value of the surge rapidly advances into the bay. At extreme water levels, the tide level dominates, and the tide levels at all points in the bay are relatively synchronized. Therefore, the maximum water level at each point is similar, and the extreme water level advances slowly into the bay.
Under the influence of typhoon Chan-hom, the surge in Hangzhou Bay is mainly of standard and mixed types. The surge curve at S1 and S2 has three obvious stages: forerunner, storm surge, and residual vibration. The standard type occurred at four other stations. Owing to the small peak value of the surge, the surge curve was relatively mild, generally in the form of tidal waves. There are small fluctuations at some moments, so they are a mixture of standard and fluctuation types.

Currents
According to the time series in Figure 3B, the main period for the storm to affect the hydrodynamic environment in the bay was from 02:00 to 13:00 on July 11, 2015 (CST). Therefore, four moments were selected during the tidal period to illustrate the  currents during the typhoon: ebb slack (02:00), flood peak (05:00), flood slack (08:00), and ebb peak (11:00). Surface currents are mostly controlled by the typhoon winds. The flow field in the typhoon center is consistent with the typhoon wind field and circulated counterclockwise around the typhoon center.
The flow field had one one-hour time lag compared with the wind field. The velocities of the flow field on the left and right sides of the typhoon center were large, whereas the velocities on the upper and lower sides were relatively small. This is because the wind field to the right of the typhoon center moves in the same direction as the typhoon, and the two are superimposed onto each other to create a higher velocity. The left side of the typhoon was near the shoreline, and water movement was impacted by the shorelines. The water channel was relatively narrow; therefore, the velocity in this area was high. The typhoon was far from the bay at the ebb slack ( Figure 5A). The coastal area of the southern East China Sea was affected by the typhoon, and the coastal current extended from northeast to southwest, with a velocity of up to 1.5 m/s. There are still ebbing tides in the inner part of the bay, but the estuary and open sea started to flood. The upstream ebbing currents meet open-sea flooding currents, forming a southward tidal current with a low velocity of approximately 0.2 m/s. At the flood peak ( Figure 5B), the offshore tidal current moved from the southeast to northwest directions, with a small velocity of approximately 0.7 m/s. When the tidal current enters the bay, the velocity is large in the narrow tidal channel and is relatively uniform inside the bay, ranging from 1.0 m/s to 1.5 m/s. At the flood slack ( Figure 5C), except for the area near the typhoon center where the velocity was high, the velocity in other sea areas was below 0.5 m/s. The offshore flow direction is basically northward, whereas the nearshore current turns northwest. Ebb tides are still dominant in Hangzhou Bay. At the peak ebb tide ( Figure 5D), coastal waters were strongly affected by the typhoon, with a maximum velocity of up to 2 m/s. The water moved southward along the coast with larger offshore currents. The water in the bay ebbed into the East China Sea with a velocity of 1.5 m/s in the narrow tidal channels at the bay mouth owing to the constraints of island topography. The flow velocity inside the bay was approximately 1 m/s. The flow directions were restricted by topography.

Waves
The distribution of the significant wave heights of typhoon waves is closely related to the distribution of typhoon wind fields ( Figure 6). Inside the bay, the significant wave height (Hs) was lower than that in the open sea. The significant wave height was higher near the southern bank and near the bay mouth during the typhoon (Figures 6A, B). The peak value of the significant wave height occurred near the mouth of HZB at 8:00 on July 11, 2015 (CST). With the approach of the typhoon wind field, the wave heights in the bay gradually increased, and the maximum wave height was 3.5 m in the bay mouth. Inside the bay, the wave height is mostly below 1.0 m. With the approach of the typhoon, the contour line of the Hs=2 m gradually expanded and extended into the bay, and the contour line of Hs=3 m was always concentrated in the central area of the bay mouth.
In the open sea, the significant wave heights were distributed in an elliptical shape, and the wave heights were the largest at the center of the ellipse and decreased towards the periphery. The wave height in the sea area to the right of the typhoon center (along the direction of typhoon movement) was generally higher than that of the sea area to the left of the typhoon center ( Figures 6C-F). The direction of the wind speed on the right side of the typhoon center was the same as the typhoon moving direction, and the wind speed after superposition of the two was greater than that on the left side of the typhoon center. Therefore, in the open sea, the development of wind waves on the right side of the typhoon center was better than that on the left side. The wave height was greater than that on the left side of the typhoon center. The distribution of the significant wave height was also restricted by the terrain. Because the center of the typhoon is close to the shoreline, the typhoon waves are easily blocked by the shoreline during the growth of the typhoon waves along the direction of the wind speed ( Figure 6D), so that the water body accumulates, and the significant wave height increases. Therefore, at this moment, large significant wave heights occurred on the right and upper sides of the typhoon center. Because the center of the typhoon is far from the coastline and the nearby sea area is open, a larger wave height occurs only on the right side of the typhoon center ( Figure 6F).

Wind-Induced Surge
During the typhoon (09:00-13:00 on July 11), the surges caused by wind stress in the inner area of Hangzhou Bay (S5 and S6) reached their maximum values of 1.27 m and 0.99 m, respectively ( Table 6). At 14:00 and 15:00 on July 11, the surge reached the central area of the bay (S3). Finally, at 17:00 on July 11, the storm surge induced by wind stress in the bay head reached a maximum value. The peak surge caused by wind stress at each station were approximately 96%, 96%, 92%, 93%, 86%, and 90% of the total storm surge, respectively. In general, there are small differences in the magnitudes of storm surges between the  northern and southern banks. The peak value of storm surge increased from the mouth to the head of the bay. This is because the reduction in the width and water depth from the mouth to the head of the bay reduces the tidal prism and amplifies the shallow water effect.
Wind stress is critical for storm surges ( Figure 7I). The windinduced surge time series coincided with the total surge curve, except for divergence during typhoon landing. The total surge was slightly larger than the wind-induced surge during typhoons. The average ratio of the peak values of wind-induced surges at the six feature points was 91%. The peak values of wind stress and surge at the six feature points were consistent with the peak values of the total surge, indicating that wind stress was the main factor leading to storm surges.
The spatial distribution of wind stress and storm surge at 17:00 on July 11 is illustrated ( Figure 8A), when the peak values of wind stress and surge occurred. At this moment, the typhoon center is located on the southeast side of the HZB. A northeasterly wind prevails in the bay, and the wind speed is approximately 20 m/s. The southern part of the typhoon center was dominated by a decrease in water, whereas the northern part was dominated by an increase in water. This pattern was mostly controlled by the wind and local morphology.
The peak of wind-induced storm surges always occurs during the low tide period (Figure 7). Case 5 was designed with the tidal level ignored at the open boundary to explore the nonlinear effect between the wind-induced storm surges and astronomical tides. In regions with a large water depth, tidal fluctuation has little

Li et al. Storm Tides in Hangzhou Bay
influence on the wind-induced surge mechanism. Therefore, at feature points S3-S6, the wind-induced surge in Case 5 was nearly the same as that in control, Case 0. For shallow regions, the water level is closely related to astronomical tides. The vertical distribution of the wind-induced shear on the water body was also affected by the water level. Therefore, when the wind speed is high and changes significantly, the surges in Cases 0 and 5 at stations S1 and S2, which are in shallow water areas, are different. The surge without tidal forcing (Case 5) was relatively gentle compared with the surge with tidal forcing (Case 0). The surge with tidal forcing showed multiple peaks, and peak values occurred at the low slack waters. This is because the depth of low slack waters is shallow, and the wind stress is stronger, leading to a larger storm surge. The wind-induced surge was controlled by both wind speed and tidal level. During typhoons, the maximum times of stress and storm surges were slightly different. In the absence of tidal forcing, the maximum wind-induced surge occurred at the time of maximum wind speed. Under the influence of tides, the maximum value of the wind-induced surge occurs in low slack waters near the peak wind speed moment.

Air Pressure-Induced Surge
During the typhoon, the air pressure exhibited clear spatial variations. The air pressure was low in the center of the typhoon and increased in the surrounding areas. By comparing the water levels of Cases 0 and 2 (Table 5), the intensity of the air pressure-induced surge was evaluated ( Table 6; Figure 7III).
The pressure-induced storm surge curve was relatively stable with no obvious peak values. The typhoon center was south of the HZB, and the southern bank was first impacted by the lowpressure area of the typhoon center. Therefore, the peak surge moment of each feature point is obviously different: the southern bank is approximately 05:00 on July 11, and the northern bank is 08:00 on July 11. During the typhoon, the pressure change in the HZB was small, so the pressure-induced storm surge was mostly uniform. The increase in water pressure caused by the decrease in air pressure was small, and the peak value of the air pressure-induced surge was relatively small compared to the peak value of the total surge. The proportions of each feature point (S1-S6) were 12%, 13%, 15%, 14%, 16%, 18%, respectively, and the average value 15%. The typhoon passes through the HZB from the outer sea, and the bay mouth is closer to the typhoon center. Therefore, the pressure-induced surges at S5 and S6 near the bay mouth were higher than those at the other four points.
Taking the time of the maximum storm surge at the bay head caused by atmospheric pressure as the reference time (08:00 on July 11), atmospheric-induced surges generally occurred in Hangzhou Bay ( Figure 8B). As the pressure was lowest at the center of the typhoon, it increased from the center to the periphery. At this time, the typhoon cyclone center was in the southeast of Hangzhou Bay, so the sea area south of Hangzhou Bay near the typhoon cyclone center had a large increase in atmospheric pressure and a water level of more than 0.3 m. The pressure-induced surge in Hangzhou Bay was largely uniform, approximately 0.2 m. The maximum pressureinduced surge was approximately 0.27 m at the head of the bay.

Wave-Induced Surge
The water levels of Cases 0 and 3 ( Table 5) were compared to evaluate the surge caused by waves in the storm surge ( Table 6; Figure 7III). The wave-induced surges were weak ( Table 6) feature points were 4, 5%, 4%, 5%, 4%, and 7%, respectively. The changes in significant wave height during typhoons are caused by wind, therefore, the peak time of the wave-induced storm surge is consistent with the peak time of wind-induced and total surges. The peak value of the wave-induced surge occurred first at the mouth of the bay and then in the middle and head of the bay. This process is consistent with the passage of the maximum wind radius. The peak moment of the wave-induced surge at the head of the bay was selected as the reference time (21:00 on July 10), and  its spatial distribution was illustrated ( Figure 8C). The maximum wind speed radius of the typhoon was at the top of the bay, and the wave surge in HZB decreased from the top to the mouth of the bay. The wave-induced surge at the top of the bay was 0.1 -0.15 m. In the middle and mouth of the bay, the coastline is relatively smooth, and the wave surge was almost the same between the northern and southern banks. The wave-induced surge in the bay was about 0.07 m and 0.04 m at the head and the mouth.
Wave-induced surges peaked in low-slack waters. Case 6 was designed while ignoring the astronomical tides to explore the nonlinear effect between wave-induced surges and astronomical tides ( Figure 7IV). The wave-induced surge without tidal forcing (case 6, blue line in Figure 7IV) was consistent with the variation in the significant wave height (dotted purple line in Figure 7IV). The greater the significant wave height, the greater the waveinduced surge. The peak of the wave-induced surge occurred slightly later than that of the significant wave height. In this test, it is assumed that the tide level is always zero and that there is always wave forcing at the boundary, which leads to water imbalance in the model domain, and the wave-induced surge is always positive. In the case of tidal forcing (Case 0), when the significant wave height was small, the wave-induced surge was close to zero. With an increase in the significant wave height during typhoons, the wave-induced surge begins to increase. There are two peaks of wave-induced surges, which all occur in low-slack waters.

Tide-Surge Interaction
Hangzhou Bay is funnel-shaped, with strong tides. To study the nonlinear interaction between tides and surges, equations (2) and (3) can be rewritten as follows: The terms on the left side are the local acceleration, advection, and Coriolis force along the estuary. The terms on the right side are the surge, barotropic, baroclinic, air pressure, friction dissipation, horizontal momentum diffusion, and wave radiation stress.
To simplify the discussion, the momentum equation was derived and transformed by ignoring the diffusion term. When only tidal forcing is considered, the equations can be rewritten as Considering only winds, air pressure, and waves, the equations can be rewritten as follows: Considering tides, wind fields, air pressure, and waves, the equations are as follows: The effects of the barotropic, baroclinic, horizontal momentum diffusion, and wave radiation stress terms are relatively small; therefore, they are ignored, and Eqs.
where ∂ U NS ∂ t and ∂ V NS ∂ t are the nonlinear local acceleration terms in the x and y directions, respectively; y x (U NS , U NS , W NS ) and y x (U NS , V NS , W NS ) are the nonlinear advection terms; -fV NS and fU NS are the nonlinear Coriolis force terms; and ( − t NS x , − t NS y ) are defined as a combination of wind stress and bottom friction terms resulting from the difference between the wind stress and bottom friction. Thus, the theoretical derivations are as follows: The above theoretical derivation establishes a relationship between the nonlinear residual levels and various influencing factors in the momentum equation. To numerically study the effects of the tide-surge interaction, we designed case 7, which was driven only by wind, waves, and astronomical tides. Cases 0, 4, and 7 were compared to obtain the time series relationship of the nonlinear surge at points P1 to P3 ( Figure 1D). The nonlinear factors based on Eqs. (31) and (32) for point P were calculated from the model results to quantitatively analyze the tide-surge interaction (Figure 9).
In the bay mouth (station P3), the nonlinear effect changes (enhance or suppress) the surge by approximately 0.5 m during the typhoon ( Figure 9A1), taking approximately 50% of the total surge. The peaks of the nonlinear effects (orange line in Figure 9A1) occur in the low slack waters of the tides (green dashed line in Figure 9A1). In the x-direction ( Figure 9C1), the nonlinear local acceleration and bottom friction terms dominate the nonlinear effect, followed by the convection term. In the ydirection ( Figure 9D1), the nonlinear local acceleration, convection, and friction terms were dominant and had comparable magnitudes. The Coriolis force term accounts for a relatively small proportion in both the x-and y-directions (yellow line). Summarizing the components in the x-and y-directions, we obtain the sums of the nonlinear increments in the two directions ( Figure 9B1). Comparing the results in the x-and ydirections with the tidal level ( Figure 9A1), the nonlinear local acceleration term enhances the nonlinear surge during low slack waters and inhibits it during high slack waters.
In the middle of the bay (station P2, Figures 9A2-D2), a similar pattern occurred at station P2, but with larger magnitudes (approximately doubled), compared with those at station P1 in the bay mouth. The Coriolis term (orange line) overtook the advection term (yellow line) in the y direction ( Figure 9D2).
In the upstream water channel of the bay (station P1, Figures 9A3-D3), the nonlinear local acceleration term decreased (blue line, Figures 9C3-D3). The bottom friction term (purple line) dominates the nonlinear effect, followed by the convection term (yellow line) in both the x-and y-directions ( Figures 9C3-D3). The Coriolis force term accounts for a very small proportion in both the x-and y-directions.

Sensitivity Test of the Model Domain
The model domain used in this study covers less than 400 km offshore. Given a typhoon radius of~50 km, the impact area is usually more than 500 km. Thus we tested the sensitivity of the model domain using two numerical tests. We run the model again using the same configurations in He et al. (2020) and Tang (2018) in Case 8. Base on Case 8, we changed the large domain to the small domain ( Figure 1A) in Case 9. The comparison of the model results is illustrated in Figure 10.
In the open sea ( Figure 10I), the significant wave heights are higher in the large domain model ( Figure 10IA, Case 8) than those in the small domain model ( Figure 10IB, Case 9). The significant wave heights in the small domain model with re-constructed wind field ( Figure 10IC, Case 0) are slightly higher than those in Case 8. Inside the HZB ( Figure 10II), the magnitudes and spatial distribution of the significant wave heights are similar in both Case 8 with large domain ( Figure 10IIA) and Case 9 with small model domain ( Figure 10IIB). The reconstructed wind field in Case 0 ( Figure 10IIC) slightly increases the significant wave height inside the HZB.
The reconstructed wind field fit the coastal wind data better, as shown in the section 2.1.2, Yu (2020) and Ren (2022). The small domain model with reconstructed wind field reproduced the storm tide and surge well, as indicted by the model validation. Hence, we used the small model domain with reconstructed wind field (Case 0) to do the simulation in this study.  Figure 1A (wind field not reconstructed), (c) in Case 0 with the small domain in Figure 1A (re-constructed wind field). Vectors indicate wind field and contours indicate significant wave heights. (I) and (II) are for the large area and the HZB, respectively.

CONCLUSIONS
This study investigates the impacts of the super typhoon Chanhom on hydrodynamics in the macro-tidal Hangzhou Bay (HZB) using a tide-surge-wave coupling numerical model. The wind field data of typhoon Chan-hom was reconstructed considering rotating and moving wind data. During typhoon Chan-hom, the Holland-Miyazaki method performed better in reconstructing the wind field in this study. The numerical model was fully validated using field data on sea surface levels, currents, and wave heights.
In HZB, wind stress was the main factor that dominated the storm surge, followed by the pressure surge effect, with the wave effect being the smallest. The surge caused by the wind stress of the three factors is the largest, and can reach 91.09% of the total surge, followed by air pressure, which is approximately 14.69% of the total surge, with the surge caused by waves being the smallest, only about 4.83%. The wind-and wave-included surges have an obvious tendency to gradually advance from the mouth to the top of the bay. At the mouth of the bay, they first peaked, followed by the middle of the bay and finally at the top of the bay. The pressureinduced surge was more evenly distributed in the HZB, and peaks appeared simultaneously. This is because the wind-and waveinduced surges are related to the change in the maximum wind speed radius of the typhoon, whereas the pressure surge is determined by the minimum pressure in the center of the typhoon.
Both the wind-and wave-induced surges are affected by the tidal level, and the peaks of both surges appear in the low-slack waters. The maximum wind-induced surge occurred near the peak moment of wind speed. The wave-induced surge peaks at the moment of maximum significant wave height. Both windand wave-induced surges occur during low slack waters. The tide-surge interaction changes (enhance or suppress) the surge by approximately 0.5 m during the typhoon, comprising approximately 50% of the total surge. Tides interact with surges via various mechanisms. The local acceleration and friction terms dominated the tide-surge interaction in the bay head, followed by the advection term. Friction and advection terms dominated at the bay head. The local acceleration term enhanced the surge in low-slack waters and suppressed it in high-slack waters.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
LL and ZH: Manuscript writing, data analysis; methodology; modelling. ZL, ZY, and YR: Material collection; modelling; data analysis. All authors contributed to the article and approved the submitted version.

FUNDING
This study was partially supported by the National Key Research a n d De v e l o p m e n t P r o g r a m o f C h i n a ( G r a n t N o . 2021YFE0206200), the National Natural Science Foundation of China (41976157, 42076177), the Science Technology D e p a r t m e n t o f Z h e j i a n g P r o v i n c e ( 2 0 2 1 C 0 3 1 8 0 , 2020C03012, U1709204).