Assessment of Extreme Storm Surges over the Changjiang River Estuary from a Wave–Current Coupled Model

Disastrous storm surges and waves caused by typhoons are major marine dynamic disasters affecting the east China coast and the Changjiang River Estuary, especially when they occur coincidentally. In this study, a high-resolution wave–current coupled model consisting of ADCIRC (Advanced Circulation) and SWAN (Simulating Waves Nearshore) was established and validated. The model shows reasonable skills in reproducing the surge levels and waves. The storm surges and associated waves are then simulated for 98 typhoons affecting the Changjiang River Estuary over the past 32 years (1987–2018). Two different wind fields, the ERA reanalysis and the ERA-based synthetic wind with a theoretical typhoon model, were adopted to discern the potential uncertainties associated with winds. Model results forced by the ERA reanalysis show comparative skills with the synthetic winds, but differences may be relatively large in specific stations. The extreme surge levels with a 50-year return period are then presented based on the coupled model results and the Gumbel distribution model. Higher risk is presented in Hangzhou Bay and the nearshore region along the coast of Zhejiang. Comparative runs with and without wave effects were conducted to discern the impact of waves on the extreme surge levels. The wave setup contributes to 2–12.5% of the 50-year extreme surge level. Furthermore, the joint exceedance probabilities of high surge levels and high wave height were evaluated with the Gumbel–logistic statistic model. Given the same joint return period, the nearshore region along the coast of Zhejiang is more vulnerable with high surges and large waves than the Changjiang River Estuary with large waves and moderate surges.


Introduction
Storm surges are perturbations in sea surface level caused by strong atmospheric disturbance. They are the major geophysical risks in coastal regions and are often associated with significant losses of life and property. The risk of storm surges is primarily linked to wind, but is also related to various other factors such as atmospherics pressure, tides, river flow and the local topography and bathymetry. When an advancing surge of water coincides with astronomical high tide or when it interferes with river flooding, the resulting increase in coastal sea levels can be further exacerbated. The Changjiang River Estuary (CRE) is one of the most economically active area of China. However, it has suffered severe damages from storm surges over the past 30 years, e.g., typhoon Chan-hom in 2015, which caused more than RMB 9.8 billion in direct economic losses [1]. Due to the large Changjiang River flooding, the complex coastline and topography and the numerous offshore islands, storm surges near the Changjiang River Estuary are especially complex and difficult to predict.
Both observational data and numerical models have been used to assess the extreme storm surges worldwide. Long-term in situ observations are capable of giving accurate estimates of sea level variation, but are usually spatially limited [2]. The well-validated numerical models, on the other hand, are powerful in producing the spatial and temporal variation of storm surges and uncover the underlying dynamics. Multiple models have been developed and configured for both global and regional applications [3][4][5][6][7][8]. The spatial resolution of the global surge model has increased to as high as 5 km in coastal regions [3]. For regional applications, the meter range resolution has been frequently adopted [9]. The performance of the simulations of the extreme sea levels, however, depends on multiple factors, including the accuracy of the meteorological forcing [10][11][12], the representation of coastal features and topography [13], the wave-current-surge interaction [8,14] and so forth.
Storm surges often occur alongside extreme wind waves. Positive storm surges combined with high tides and wind waves can cause coastal floods and saltwater intrusion, which are probably the most destructive natural hazards of geophysical origin [9,15]. The wave-current interaction during an extreme system may also change both the surge level and the wave height substantially. The wave-induced radiation stress acts as an additional force in the momentum equations that modify coastal water circulations, resulting in wave setup/down near the coast. Funakoshi et al. [16] revealed the importance of the wave radiation stress in the prediction of storm surges and found that it contributes to 10-15% of the peak water level. Using a coupled model, Xie et al. [17] showed a wave setup of 0.2 m to the surge level in the Gulf of Maine during the Patriots Day storm, and confirmed the importance of wave radiation stress. In a case study in the Taiwan Strait, Yu et al. [18] found that wave setup may contribute up to 24% of the total storm surge with a maximum of 0.28 m along the east China coast during typhoon Morakot. The wave setup can potentially affect the surge level over the CRE. Wang et al. [14] found that the wave setup reached a maximum of~0.4 m in Jiaojiang Bay right south of the CRE during typhoon Saomai and typhoon Chan-hom.
Many efforts have been made to understand the dynamics of storm surge in regions adjacent to the CRE and revealed the importance of several different factors. Hu et al. [19] found that the effects of remote wind should be considered when stimulating storm surges adjacent to the CRE. Guo et al. [20] investigated the model sensitivity of storm surges to the conduction of theoretical typhoon wind and found that the radius to maximal wind plays a key role in predicting peak surge levels in Hangzhou Bay. Pan and Liu [8] also demonstrated the importance of the radius to maximal wind in setting up the surge level in Hangzhou Bay. Wang et al. [14] demonstrated the importance of waves on surge levels. Among these studies, the model with unstructured mesh, due to its efficiency in resolving the complex geometry and bathymetry, has been widely used [14,21]. The unstructured circulation model ADCIRC and wave model SWAN have been coupled and used in this study. The effects of wave-current interaction on surges and waves are thus being considered.
Estimations of extreme surge levels and wave heights are important for coastal management and design of coastal structures. In general, the extreme waves and surges with a given return period can be estimated using a particular statistical method. For example, the Weibull distribution [22] and the Gumbel distribution [23] are commonly used to estimate extreme wave heights or surges. Feng and Jiang [24] indicated that the Gumbel distribution performs well on estimation of the extreme surge level along the coast of the Northwestern Pacific Ocean. The coastal engineering design commonly needs to face many complex coastal phenomena, such as overtopping resulting from combination of two physical processes; thus, the univariate extreme models should be extended to multivariate joint probability method [25]. There have been several studies on the optimalization of the multivariate joint probability method (JPM). It has been shown that the JPM is more suitable for consistency of estimates for different designs and is more complicated than the univariate extreme models [26]. The JPM thus has been widely used in estimating extreme water levels and waves [27][28][29][30].
There are also several studies analyzing the long-term hazard near the CRE. Li et al. [31] gave the broad scale distributions of the extreme surges and waves over the Northwest Pacific using the univariate Gumbel method. Chen et al. [32] applied the JPM to nine selected stations along the coast of China and found that the traditional empirical method overestimates the joint extreme surge levels and waves as compared to the JPM, and suggested the use of the JPM to achieve realistic estimations. Based on coupled model results and the Quadrature Probability Method, Yin et al. [29] estimated the 100-year and 200-year annual maximum water levels at three selected stations in the CRE. However, the generation of synthetic typhoon tracks were relatively limited. Though waves are considered in places, their contribution to the estimation of extreme surge levels is still unclear. Additionally, surges and waves experience different variations when approaching nearshore, especially when complex topography and islands are placed [33]; thus, it is more interesting to evaluate the joint exceedance probabilities of high surge levels and large waves. Given this, high-resolution wave-current coupled models are needed to resolve the irregular coastline and bathymetry adjacent to CRE. Long-term wave-current hindcast simulations adopting numerous typhoon tracks are needed to discern the effects of waves on the estimation of extreme surge levels. The JPM will be adopted to reveal the joint exceedance probability of extreme surges and waves.
In this study, 98 typhoon storm surges and the associated waves affecting the CRE from 1987 to 2018 are simulated using the coupled circulation and wave model, ADCIRC (Advanced Circulation) and SWAN (Simulating Waves Nearshore). Different atmospheric forcings, the ERA-interim reanalysis wind and the ERA-based synthetic wind are compared and adopted to simulate the surges and waves adjacent to CRE. Both coupled and uncoupled runs are conducted to discern the effects of wave-current interactions on surge levels and on estimation of extreme surge levels. Both the traditional Gumbel model and the joint Gumbel-logistical model are used to estimate the extreme surge levels and waves, and the joint exceedance probabilities over the CRE and the adjacent region.
The rest of the paper is organized as follows. Section 2 gives a brief introduction of the coupled ADCIRC-SWAN model and statistical methods used to deduce the extreme surge levels and waves. The model performance is assessed in Section 3. Results and discussions are presented in Sections 4 and 5, respectively. Conclusions are summarized in Section 6.

Study Area
This study focuses on the CRE and the adjacent coastal regions ( Figure 1). The Shanghai city is situated on the eastern fringe of the Changjiang River Delta, with the Jiangsu province to the north and the Zhejiang province to the south. The Hangzhou Bay is located right south of the CRE, with multi-islands sitting at the bay mouth, i.e., the Zhousan islands. The CRE experience a strong spring neap tide, with an average range of about 3 m [34]. The tidal range may be distorted by the Changajing River discharge. The Changjiang River discharge varies significantly from the maximum of about 5 × 10 4 m 3 /s in summer to about 1.0 × 10 4 m 3 /s in winter [35,36]. A high river discharge may reduce the tidal range and increase the tidal deformation compared to low river discharge [37], thus affecting the sea level prediction. The wave field is generally dominated by local wind seas, but swells from the outer shelf might occasionally be important [38]. Climate over the CRE is mainly governed by the East Asian monsoon. During the summer season, the southeast wind dominates and waves have a dominant northwest direction but are relatively small, except during the occasional typhoon events. During winter, the strong northwest wind generates larger waves mainly in the southeast direction [39]. Precipitation over the CRE is also dominated by the East Asian Monsoon and shows strong seasonal variations, with an average amount of approximately 1000-1400 mm yr −1 [40,41]. However, the precipitationdriven river discharge during typhoon events may be quite large, contributing significantly to water levels inside the estuary [42]. strong seasonal variations, with an average amount of approximately 1000-1400 mm yr −1 [40,41]. However, the precipitation-driven river discharge during typhoon events may be quite large, contributing significantly to water levels inside the estuary [42]. The study area and locations of tide gauge stations (squares), wave buoys (stars) and stations with surge levels available (dots). A1−A8 (dots) represents Lvsi, Dajishan, Tanhu, Zhenhai, Shipu, Dachen, Kanmen and Wenzhou, respectively. B1−B2 (stars) denote the two wave buoys. C1−C4 (squares) represent Baozhen, Wusongkou, Gaoqiao station and Luchaogang, respectively.

Model Description
The oceanic circulation model used is the ADvanced CIRculation (ADCIRC) Model [43,44]. ADCIRC is a depth-integrated, shallow water, finite-element model capable of simulating tidal circulations and storm surge propagations. It solves the Navier-Stokes equation over unstructured meshes with triangular finite elements of varying sizes to represent complex coastal features, barrier islands and internal barriers. It has been widely used in simulations of astronomical tides, wind-induced circulations and storm surges [9,14,45].
The wave model is Simulating Waves Nearshore (SWAN) [46]. SWAN is a thirdgeneration shallow water spectral wave model including wave propagation, refraction due to currents and depth, generation by wind, dissipation (white-capping, bottom friction, depth induced breaking) and nonlinear wave-wave interactions [46]. It is designed to obtain realistic estimates of wave parameters in coastal areas, lakes and estuaries from given wind and current conditions, and is effective in resolving complex coastline and bathymetry by using the unstructured triangular grid.
The ADCIRC and the SWAN models are coupled on an identical, unstructured mesh, share the parallel computing infrastructure and run sequentially in time [47]. The coupled ADCIRC-SWAN has been widely used to simulate the surge levels during extreme weather events, to investigate the effects of wave-current interactions, to reveal the role of shelf geometry on surge-tide-wave interactions and so forth [13,14]. Its efficiency and superior performance to the ADCIRC-only model have been demonstrated elsewhere [14,48]. The wave setup has been identified to be an important parameter that must be incorporated in an operational storm surge forecasting system [11,49].
In the present application, the computational model domain covers the Changjiang River Estuary, the East China Sea and part of the North Pacific Ocean (Figure 2). The unstructured mesh includes 113,570 computational nodes and 216,950 triangle units. High

Model Description
The oceanic circulation model used is the ADvanced CIRculation (ADCIRC) Model [43,44]. ADCIRC is a depth-integrated, shallow water, finite-element model capable of simulating tidal circulations and storm surge propagations. It solves the Navier-Stokes equation over unstructured meshes with triangular finite elements of varying sizes to represent complex coastal features, barrier islands and internal barriers. It has been widely used in simulations of astronomical tides, wind-induced circulations and storm surges [9,14,45].
The wave model is Simulating Waves Nearshore (SWAN) [46]. SWAN is a thirdgeneration shallow water spectral wave model including wave propagation, refraction due to currents and depth, generation by wind, dissipation (white-capping, bottom friction, depth induced breaking) and nonlinear wave-wave interactions [46]. It is designed to obtain realistic estimates of wave parameters in coastal areas, lakes and estuaries from given wind and current conditions, and is effective in resolving complex coastline and bathymetry by using the unstructured triangular grid.
The ADCIRC and the SWAN models are coupled on an identical, unstructured mesh, share the parallel computing infrastructure and run sequentially in time [47]. The coupled ADCIRC-SWAN has been widely used to simulate the surge levels during extreme weather events, to investigate the effects of wave-current interactions, to reveal the role of shelf geometry on surge-tide-wave interactions and so forth [13,14]. Its efficiency and superior performance to the ADCIRC-only model have been demonstrated elsewhere [14,48]. The wave setup has been identified to be an important parameter that must be incorporated in an operational storm surge forecasting system [11,49].
In the present application, the computational model domain covers the Changjiang River Estuary, the East China Sea and part of the North Pacific Ocean (Figure 2). The unstructured mesh includes 113,570 computational nodes and 216,950 triangle units. High resolution is placed in the Changjiang River Estuary to resolve the complex coasts and island shorelines. The resolution is about 200 m in the estuary but increases gradually to 50 km at the open ocean. The bathymetry over the inner shelf is provided by Shanghai Municipal Oceanic Bureau, while for the open ocean, data from the GEBCO (General Bathymetric Chart of the Oceans) were used (Figure 1a). The circulation model is forced by tides at the offshore boundary and by river flow at the upstream of the Changjiang River Estuary. Tidal forcing at the open boundary is specified using the Oregon State University (OSU) global inverse tidal model of TPXO7.2 [50]. Eight major tidal constituents, M 2 , S 2 , N 2 , K 2 , K 1 , O 1 , P1 1 and Q 1 , are considered. At the upstream boundary of the Changjiang River Estuary, the momentum boundary condition is imposed on the river discharge divided by the upstream cross-sectional area. The Changjiang River discharge is set to a typical value in summer of about 5 × 10 4 m 3 /s. The SWAN wave model is run and coupled to the same grid as the ADCIRC model. Twenty-five frequencies (0.03-1.0 Hz) and 36 directional bands are used. Wave growth by wind is computed with the exponential term of Komen et al. [51] with a default coefficient. Different time steps are used for ADCIRC and SWAN. In ADCIRC, a 3 s barotropic timestep is used. SWAN is run in non-stationary mode with a time step of 300 s. The coupling between ADCIRC and SWAN is carried out synchronously with a 30 min time interval. The ADCIRC model offers the water level and the oceanic currents to the wave model, and takes the wave parameters to calculate the radiation stress, which further affect the water levels and oceanic currents [47]. resolution is placed in the Changjiang River Estuary to resolve the complex coasts and island shorelines. The resolution is about 200 m in the estuary but increases gradually to 50 km at the open ocean. The bathymetry over the inner shelf is provided by Shanghai Municipal Oceanic Bureau, while for the open ocean, data from the GEBCO (General Bathymetric Chart of the Oceans) were used (Figure 1a). The circulation model is forced by tides at the offshore boundary and by river flow at the upstream of the Changjiang River Estuary. Tidal forcing at the open boundary is specified using the Oregon State University (OSU) global inverse tidal model of TPXO7.2 [50]. Eight major tidal constituents, M2, S2, N2, K2, K1, O1, P11 and Q1, are considered. At the upstream boundary of the Changjiang River Estuary, the momentum boundary condition is imposed on the river discharge divided by the upstream cross-sectional area. The Changjiang River discharge is set to a typical value in summer of about 5 × 10 4 m 3 /s. The SWAN wave model is run and coupled to the same grid as the ADCIRC model. Twenty-five frequencies (0.03-1.0 Hz) and 36 directional bands are used. Wave growth by wind is computed with the exponential term of Komen et al. [51] with a default coefficient. Different time steps are used for ADCIRC and SWAN. In ADCIRC, a 3 s barotropic timestep is used. SWAN is run in non-stationary mode with a time step of 300 s. The coupling between ADCIRC and SWAN is carried out synchronously with a 30 min time interval. The ADCIRC model offers the water level and the oceanic currents to the wave model, and takes the wave parameters to calculate the radiation stress, which further affect the water levels and oceanic currents [47].

Diagnostic Typhoon Model
The wind field and the pressure field used in this study are from the global atmospheric reanalysis products ERA-interim provided by the European Center for Medium Range Weather Forecast, with a spatial resolution of 0.125° × 0.125° and temporal resolution of 6 h. To account for the possible underestimation of the typhoon wind in ERA reanalysis, a diagnostic typhoon model was also adopted and merged with ERA reanalysis

Diagnostic Typhoon Model
The wind field and the pressure field used in this study are from the global atmospheric reanalysis products ERA-interim provided by the European Center for Medium Range Weather Forecast, with a spatial resolution of 0.125 • × 0.125 • and temporal resolution of 6 h. To account for the possible underestimation of the typhoon wind in ERA reanalysis, a diagnostic typhoon model was also adopted and merged with ERA reanalysis wind (ERA-Based Synthetic Wind, ERA-BSW hereafter). The diagnostic typhoon model used in this study is a circular typhoon model introduced by Jelesnianski [52,53]. Previous studies indicated that it is proper to give a good estimation of typhoon wind in the East China Sea [31,54]. The Jelesnianski circular typhoon model calculates the air pressure and the wind fields as follows: where (x c , y c ) and (x, y) represent the coordinates of typhoon center and the target point; R max is the maximum wind radius; V max is the maximum wind speed; r is the distance between typhoon center and the target grid point; P 0 is the typhoon central pressure; P a is pressure at the target point, P ∞ is peripheral pressure setting to 1013.25 hPa; θ is the flow angle setting to 20 • ; and V ox and V oy are the typhoon's moving speed. The information of typhoons, i.e., the moving track (x c , y c ) and central pressure (P 0 ), are derived from Best Track Dataset provided by the Japan Meteorological Agency [55]. Following Atkinson and Holliday [56], the maximum wind speed V max is calculated as: The maximum wind radius is calculated as follows [57]: where ϕ denotes the latitude of the typhoon center. The aforementioned diagnostic wind and pressure are then merged to the ERA reanalysis through the following method [58]: where V Q represents the ERA reanalysis wind/pressure, V M is the diagnostic wind/pressure and e is the weight coefficient defined as where n is an adjustable coefficient, set to 10 in the East China Sea [59].

Extreme Value Analysis
The Gumbel distribution provides a mathematical framework for univariable longterm factor statistical analysis and has been widely used for predicting occurrences of natural extreme events such as water levels and wave heights in coastal engineering [60][61][62]. It assumes that the maximum data come from a large number of observations. A Gumbel distribution function is defined as in Castillo et al. [63]: where A and B are the location and the scale parameter, respectively. The corresponding cumulative distribution function is given as: The formula for a given return period T and corresponding probability is: where λ is the average number of samples used each year. The annual maximum taken from either observations or models is fitted to the Gumbel distribution to estimate the maximum surge level of different return periods.
In addition to the univariate extreme analysis, the joint probability analysis of surges and waves are also presented. The Gumbel-logistic model is adopted in this study. This approach is based on the standard Gumbel distribution as a variable marginal distribution, and the logistic correlation model is used to correlate two variables of extreme surges and waves [30,64,65]. This approach has been used with success in the estimation of joint exceedance probability of extreme surges and waves in the Bohai Sea and the Shandong Peninsula [30,66].
The joint probability function is defined as: where P is the joint probability value when the annual maximum surge ζ and the corresponding significant wave height H s are larger than x 1 and x 2 , respectively. α denotes the association between ζ and H s , and is estimated as follows [54]: where r 12 is the correlation coefficient between ζ and H s . The marginal distribution F(x i ) is defined as:

Tide and Wave Validation
The model predicted tidal heights were first validated against observations.  [67]). We have also calculated the co-tidal and co-amplitude charts for the M 2 tide K 1 tide (not shown). The predicted co-tidal charts are in agreement with observations obtained from 10 years of TOPEX-Poseidon altimetry measurements [68]. The significant wave heights from the wave-current coupled models are also compared with observations. The significant wave height is derived from the ADCIRC-SWAN coupled run forced by ERA reanalysis wind. Due to the limited buoy wave observations, waves are compared during the summer season of 2018, when observations are available. There are five typhoons affecting the Changjiang River Estuary during the time of comparison. Figure 4 shows the observed and predicted significant wave height The significant wave heights from the wave-current coupled models are also compared with observations. The significant wave height is derived from the ADCIRC-SWAN coupled run forced by ERA reanalysis wind. Due to the limited buoy wave observations, waves are compared during the summer season of 2018, when observations are available. There are five typhoons affecting the Changjiang River Estuary during the time of comparison. Figure 4 shows the observed and predicted significant wave height at two buoys, whose locations are shown in Figure 1b. The simulated significant wave heights agree with the observations closely, especially at buoy B2. The maximum significant wave heights during typhoons are reproduced reasonably well by the wave model. The averaged rms error, relative error, correlation coefficient and skill are 27.6 cm, 10.72%, 0.90 and 0.94 at buoy B1, respectively, and 21.5 cm, 5.5%, 0.95 and 0.97 at buoy B2, respectively. It should be noted that the high-frequency fluctuation of waves at buoy B1 is not presented in the wave-current simulation, indicating that it originated not only from the tidal oscillation but also the sea breeze, which is not resolved by the ERA reanalysis wind.

Storm Surges Validation
In addition to the tides and waves, we also examined the model prediction of su levels. There are about 98 typhoons affecting the CRE during1987 to 2018. Accordin their tracks, the typhoons affecting the CRE can be divided into four general catego [69]. Most (46 of 98, 47.4%, track I) of the typhoons passed through the east of the C without landfall, e.g., typhoon 0205 in Figure 5, and 41.8% (41 of 98, track II) made la fall south of the CRE, and then moved northwestward into the mainland, e.g., typh 0509 in Figure 5. A third category (7 of 98, 7.1%, track III) landed shortly but re-ent the sea and moved further north along the east China coast, e.g., typhoon 1509 in Fig  5. Very few (4 of 98, 4.1%, track IV) typhoons made landfall right at the CRE (typh 1810 in Figure 5). Four different typhoons were selected for comparison accordin their tracks. The surge levels are derived from the difference between the ADCI SWAN coupled run forced by ERA reanalysis wind and the tide-only run without w and waves. Figure 6 shows the comparison of observed and simulated surge levels at eight gauge stations along the east China coast near the Changjiang River Estuary. It ca seen from Figure 6 that the simulated peak surges are in good agreement with obse

Storm Surges Validation
In addition to the tides and waves, we also examined the model prediction of surge levels. There are about 98 typhoons affecting the CRE during1987 to 2018. According to their tracks, the typhoons affecting the CRE can be divided into four general categories [69]. Most (46 of 98, 47.4%, track I) of the typhoons passed through the east of the CRE without landfall, e.g., typhoon 0205 in Figure 5, and 41.8% (41 of 98, track II) made landfall south of the CRE, and then moved northwestward into the mainland, e.g., typhoon 0509 in Figure 5. A third category (7 of 98, 7.1%, track III) landed shortly but re-entered the sea and moved further north along the east China coast, e.g., typhoon 1509 in Figure 5. Very few (4 of 98, 4.1%, track IV) typhoons made landfall right at the CRE (typhoon 1810 in Figure 5). Four different typhoons were selected for comparison according to their tracks. The surge levels are derived from the difference between the ADCIRC-SWAN coupled run forced by ERA reanalysis wind and the tide-only run without winds and waves.
simulations show reasonable agreement with observations. The relative error is generally less than 18%, with a relatively larger error at Lvsi, and the mean is about 10%. The maximum surges produced by ERA reanalysis perform slightly better than the ERA-BSW. Model results forced by ERA reanalysis are presented in the following sections; results forced by the ERA-BSW are also presented where necessary. Potential effects of errors due to different winds are discussed in Section 5.    (Figure 6g,h). Table 1 presents the mean relative error of the peak surge levels for all 98 typhoons whenever the observations are available. The relative error is generally less than 15%, with a mean of about 7%.

Extreme Surge Levels
We first examined the simulated maximum surge levels over the past 30 years, i.e., the highest surge levels for all typhoon storms. Maximum surge levels for each category of typhoon tracks are given in Figure 7 for comparison. As shown in Figure 7a, track I, which passes through the east of CRE without landfall, generates higher surges along the whole east China coast, and is the only track generating higher surges north of the CRE. The maximum surge level in track I is about 2 m. Track II, which makes landfall south of the CRE, generates the largest surges of about 3.2 m over In addition to the run forced by ERA reanalysis wind, model results forced by ERA-BSW are also compared with available observations. As shown in Table 1, both model simulations show reasonable agreement with observations. The relative error is generally less than 18%, with a relatively larger error at Lvsi, and the mean is about 10%. The maximum surges produced by ERA reanalysis perform slightly better than the ERA-BSW. Model results forced by ERA reanalysis are presented in the following sections; results forced by the ERA-BSW are also presented where necessary. Potential effects of errors due to different winds are discussed in Section 5.

Extreme Surge Levels
We first examined the simulated maximum surge levels over the past 30 years, i.e., the highest surge levels for all typhoon storms. Maximum surge levels for each category of typhoon tracks are given in Figure 7 for comparison. As shown in Figure 7a, track I, which passes through the east of CRE without landfall, generates higher surges along the whole east China coast, and is the only track generating higher surges north of the CRE. The maximum surge level in track I is about 2 m. Track II, which makes landfall south of the CRE, generates the largest surges of about 3.2 m over the Hangzhou Bay (Figure 7b). Surges over the CRE are generally lower than in Hangzhou Bay and are about 1-2 m in track II (Figure 7b). Track III, which lands shortly but re-enters the sea, has a similar pattern of maximum surge levels as track II, but the magnitude is about 1 m lower than track II (Figure 7c). Surprisingly, typhoons making landfall right at the CRE generate the lowest surge levels. The historical maximum surges due to track IV are generally less than 1.0 m over the CRE (Figure 7d). the Hangzhou Bay (Figure 7b). Surges over the CRE are generally lower than in Hangzhou Bay and are about 1-2 m in track II (Figure 7b). Track III, which lands shortly but re-enters the sea, has a similar pattern of maximum surge levels as track II, but the magnitude is about 1 m lower than track II (Figure 7c). Surprisingly, typhoons making landfall right at the CRE generate the lowest surge levels. The historical maximum surges due to track IV are generally less than 1.0 m over the CRE (Figure 7d).    and higher south of the CRE, indicating the higher risk along the coast of Zhejiang. The maximum 50-year extreme surge level is located in the upstream of Hangzhou Bay, but several other peaks are also presented, with a range approximately the same as in Hangzhou Bay. The locations of the surge peaks are also plotted in Figure 8. It is noted that the surge peaks are mainly located in regions where the narrowest channel upstream is reached. This highlights the importance of accurately representing the channel geometry in surge forecasting models.   and higher south of the CRE, indicating the higher risk along the coast of Zhejiang. The maximum 50-year extreme surge level is located in the upstream of Hangzhou Bay, but several other peaks are also presented, with a range approximately the same as in Hangzhou Bay. The locations of the surge peaks are also plotted in Figure 8. It is noted that the surge peaks are mainly located in regions where the narrowest channel upstream is reached. This highlights the importance of accurately representing the channel geometry in surge forecasting models.   The 50-year extreme surge levels along the east China coast are further shown in Figure 9. The vertical lines in Figure 9 indicate the location of the Changjiang River Estuary. As shown in Figure 9, the 50-year extreme surge levels are lower north of the CRE and higher south of the CRE, indicating the higher risk along the coast of Zhejiang. The maximum 50-year extreme surge level is located in the upstream of Hangzhou Bay, but several other peaks are also presented, with a range approximately the same as in Hangzhou Bay. The locations of the surge peaks are also plotted in Figure 8. It is noted that the surge peaks are mainly located in regions where the narrowest channel upstream is reached. This highlights the importance of accurately representing the channel geometry in surge forecasting models.

Wave Effects on Evaluation of Extreme Surge Levels
As waves propagate from the continental shelf into the shallow waters, the excessive wave momentum results in wave setup near the coast. The importance of wave setup to surge levels has been recognized in coastal regions frequently affected by typhoons [13,18,70,71]. However, their roles in the estimation of the extreme surge levels have not been well documented, and are considered to account for the large model-data discrepancy [5].
To discuss the effects of wave-current interaction on surge levels, all 98 storm surge processes affecting in the CRE were further run without waves, i.e., with ADCIRC only. The wave setup is calculated as the difference of water level between the ADCIRC-SWAN coupled run and the ADICIRC only run. Figure 10 gives an example of wave setup during typhoon Rammasun (0205). As shown in Figure 10, the wave setup locates in the very nearshore area, while the wave set-down distributes over the offshore area, consistent with former studies [18]. The mean wave setup is about 0.15 m along the coast of Zhejiang, with a maximum of about 0.18 m around latitude 29 • N. The wave setup at CRE is lower than the coast of Zhejiang but has discernable magnitude of about 5 cm.

Wave Effects on Evaluation of Extreme Surge Levels
As waves propagate from the continental shelf into the shallow waters, the excessive wave momentum results in wave setup near the coast. The importance of wave setup to surge levels has been recognized in coastal regions frequently affected by typhoons [13,18,70,71]. However, their roles in the estimation of the extreme surge levels have not been well documented, and are considered to account for the large model-data discrepancy [5].
To discuss the effects of wave-current interaction on surge levels, all 98 storm surge processes affecting in the CRE were further run without waves, i.e., with ADCIRC only. The wave setup is calculated as the difference of water level between the ADCIRC-SWAN coupled run and the ADICIRC only run. Figure 10 gives an example of wave setup during typhoon Rammasun (0205). As shown in Figure 10, the wave setup locates in the very nearshore area, while the wave set-down distributes over the offshore area, consistent with former studies [18]. The mean wave setup is about 0.15 m along the coast of Zhejiang, with a maximum of about 0.18 m around latitude 29° N. The wave setup at CRE is lower than the coast of Zhejiang but has discernable magnitude of about 5 cm.  Figure 11a, and differences in the estimation with waves are shown in Figure 11b. The 50-year extreme sea levels without waves are underestimated in the nearshore area, and are well confined by the 20 m isobath. The contribution due to wave setup reaches as high as 0.25 m along the coast of Zhejiang, but is less than 0.1 m at Hangzhou Bay and the Changjiang The 50-year extreme surge levels without wave effects are shown in Figure 11a, and differences in the estimation with waves are shown in Figure 11b. The 50-year extreme sea levels without waves are underestimated in the nearshore area, and are well confined by the 20 m isobath. The contribution due to wave setup reaches as high as 0.25 m along the coast of Zhejiang, but is less than 0.1 m at Hangzhou Bay and the Changjiang River Estuary. Wave effects north of the Changjiang River Estuary are further weakened, with a wave contribution of less than 0.05 m. The effects of the wave setup along the east China coast are also shown in Figure 9. As shown in Figures 9 and 11b, the wave setup contributes to about 2-12.5% of the overall extreme surge level along the east China coast. In contrast to the raised extreme surge levels in the nearshore region due to wave setup, the extreme surge levels are lowered in the offshore region due to wave set-down. However, the offshore drop is less than the nearshore raise, with a maximum of about −0.1 m, approximately two times smaller than the nearshore uplift. However, the offshore drop is less than the nearshore raise, with a maximum of about −0.1 m, approximately two times smaller than the nearshore uplift.

Joint Probability Analysis of Extreme Wave Heights and Surges
Oceanic dynamic hazards having the largest impact are events that have both high surge levels and large wave heights. The two generally occur synchronously. Thus, it is more interesting to evaluate the joint exceedance probabilities. Joint probability analysis provides additional information that cannot be obtained by single variable storm surge analysis, such as the joint return periods, and the conditional return periods for these variables [30,64,72]. The Gumbel logistic model is used to obtain the joint probability distribution of surge levels and significant wave heights from the wave-current coupled models. Eight points along the East China coast are selected to show their joint return periods, with their location marked in Figure 1b. Figure 12 shows contours of joint return periods of storm surge levels and significant wave heights at the selected stations. Given the return period of a storm event, it can be seen from Figure 12 that both extreme surge levels and significant wave heights increase with the increasing return period. However, the changing rate varied at different points. For example, when the extreme surge level of 1.2 m with 10-year return period increases to 1.7 m with 100-year return period, the corresponding significant wave height increases from 4 m to 5.5 m at Kanmen. While in Wenzhou, given the similar change of extreme surge level from 1.1 m with 10-year return period to 1.7 m with 100-year return period, the corresponding significant wave height only increases from 2.2 m to 3.2 m. This indicates that Kamen is more likely to be affected by both tremendous surges and waves, while tremendous surges are accompanied with moderate waves at Wenzhou.
In addition, the 50-year extreme surge levels when the co-occurred significant wave height reached 1.0 m are shown in Figure 13a. Higher surges of about 1.5-2.0 m are found at the Hangzhou Bay and coasts along Zhejiang province. Moderate surge levels (1.0-1.5 m) occur over the CRE region. The conditional joint return periods with the

Joint Probability Analysis of Extreme Wave Heights and Surges
Oceanic dynamic hazards having the largest impact are events that have both high surge levels and large wave heights. The two generally occur synchronously. Thus, it is more interesting to evaluate the joint exceedance probabilities. Joint probability analysis provides additional information that cannot be obtained by single variable storm surge analysis, such as the joint return periods, and the conditional return periods for these variables [30,64,72]. The Gumbel logistic model is used to obtain the joint probability distribution of surge levels and significant wave heights from the wave-current coupled models. Eight points along the East China coast are selected to show their joint return periods, with their location marked in Figure 1b. Figure 12 shows contours of joint return periods of storm surge levels and significant wave heights at the selected stations. Given the return period of a storm event, it can be seen from Figure 12 that both extreme surge levels and significant wave heights increase with the increasing return period. However, the changing rate varied at different points. For example, when the extreme surge level of 1.2 m with 10-year return period increases to 1.7 m with 100-year return period, the corresponding significant wave height increases from 4 m to 5.5 m at Kanmen. While in Wenzhou, given the similar change of extreme surge level from 1.1 m with 10-year return period to 1.7 m with 100-year return period, the corresponding significant wave height only increases from 2.2 m to 3.2 m. This indicates that Kamen is more likely to be affected by both tremendous surges and waves, while tremendous surges are accompanied with moderate waves at Wenzhou.   Figure 13b. As shown in Figure 13b, the broad offshore area in red indicates that the conditional joint return periods are longer than 100-year period. The conditional joint return periods are longer than the 100-year period right over the CRE and north of CRE, but shortened to about 10-30 years along the coasts of Zhejiang, indicating the higher risk over this region.

Discussion
Due to limitations of long-term observations of both sea levels and waves, models are widely used to evaluate the single and joint extreme events. However, the model performance depends on model dynamics, the representation of topography, the atmospheric forcing, etc. Both satellite retrievals and models tend to underestimate the extreme wind during typhoon storms [30,32]. Synthetic winds are commonly adopted to overcome the deficiency. The synthetic method, though widely used, fails to resemble the typhoon asymmetric structure [73]. As shown in Section 3.2, comparisons between model simulations forced by ERA-BSW improve the model skill at Tanhu, Dachen, but become worse in the other three stations. Figure 14 shows the annual extreme surge levels at Tanhu, which have relatively long observations. As shown in Figure 14, the simulation forced by synthetic wind overestimates the annual extreme surge levels in some years, with the maximum discrepancy exceeding 0.7 m in 1988, but may also underestimate it for other years, e.g., 2007. The discrepancy would pass to the statistical analysis of long-term extreme surge levels. The five tidal gauge stations with relatively long-term observations are analyzed to investigate the potential effects due to the error transmission of winds. Figure 15 gives the estimation of the extreme surge levels with a 100-year return period derived from both observations and models. The extreme surge level forced by ERA reanalysis resembles the observations slightly better than the ERA-BSW. The ERA-BSW overestimates the extreme surge level in some stations, but underestimates it in others, and so does ERA. The uncertainty of the extreme surge levels resultant from different wind forcings can be relatively large at specific locations. For example, the extreme surge levels from simulation with ERA reanalysis are slightly higher than the estimation from ERA-BSW at Dajishan and Tanhu. In contrast, the ERA-BSW wind overestimates the 100-year surge levels at Kanmen by about 0.6 m, but underestimates it by about 0.7 m at Lvsi, with a relative error of about 30%.

Discussion
Due to limitations of long-term observations of both sea levels and waves, models are widely used to evaluate the single and joint extreme events. However, the model performance depends on model dynamics, the representation of topography, the atmospheric forcing, etc. Both satellite retrievals and models tend to underestimate the extreme wind during typhoon storms [30,32]. Synthetic winds are commonly adopted to overcome the deficiency. The synthetic method, though widely used, fails to resemble the typhoon asymmetric structure [73]. As shown in Section 3.2, comparisons between model simulations forced by ERA-BSW improve the model skill at Tanhu, Dachen, but become worse in the other three stations. Figure 14 shows the annual extreme surge levels at Tanhu, which have relatively long observations. As shown in Figure 14, the simulation forced by synthetic wind overestimates the annual extreme surge levels in some years, with the maximum discrepancy exceeding 0.7 m in 1988, but may also underestimate it for other years, e.g., 2007. The discrepancy would pass to the statistical analysis of long-term extreme surge levels. The five tidal gauge stations with relatively long-term observations are analyzed to investigate the potential effects due to the error transmission of winds. Figure 15 gives the estimation of the extreme surge levels with a 100-year return period derived from both observations and models. The extreme surge level forced by ERA reanalysis resembles the observations slightly better than the ERA-BSW. The ERA-BSW overestimates the extreme surge level in some stations, but underestimates it in others, and so does ERA. The uncertainty of the extreme surge levels resultant from different wind forcings can be relatively large at specific locations. For example, the extreme surge levels from simulation with ERA reanalysis are slightly higher than the estimation from ERA-BSW at Dajishan and Tanhu. In contrast, the ERA-BSW wind overestimates the 100-year surge levels at Kanmen by about 0.6 m, but underestimates it by about 0.7 m at Lvsi, with a relative error of about 30%.  To further investigate the effects of model errors on the estimation of extreme surge levels, we constructed a series of samples with ±1-20% offset of the original surge levels,  To further investigate the effects of model errors on the estimation of extreme s levels, we constructed a series of samples with ±1-20% offset of the original surge le To further investigate the effects of model errors on the estimation of extreme surge levels, we constructed a series of samples with ±1-20% offset of the original surge levels, and re-estimated the extreme surge level with a 100-year return period. Figure 16 shows the relative offset of the 100-year extreme surge level to the original assessment when positive or negative errors are added to the original samples at Tanhu. As shown in Figure 14, the observed annual maximum surge level ranges between 0.26 m and 1.79 m, and the original 100-year extreme surge level is estimated to be 1.9 m. Errors increase with the increased error ratio (y-coordinate), and model errors in higher surge events result in the largest uncertainty. The relative error of the 100-year extreme surge level will increase from 4% to 8% when the error with the highest surge level increases from 10% to 20%. Storm events with moderate surge levels (0.8-1.4 m) would also affect the estimation of the extreme surge levels, but the uncertainty is generally less than 4%. Thus, specific attention should be paid to typhoon events with higher surge levels, which dominate the estimation of extreme surge levels. and re-estimated the extreme surge level with a 100-year return period. Figure 16 shows the relative offset of the 100-year extreme surge level to the original assessment when positive or negative errors are added to the original samples at Tanhu. As shown in Figure 14, the observed annual maximum surge level ranges between 0.26 m and 1.79 m, and the original 100-year extreme surge level is estimated to be 1.9 m. Errors increase with the increased error ratio (y-coordinate), and model errors in higher surge events result in the largest uncertainty. The relative error of the 100-year extreme surge level will increase from 4% to 8% when the error with the highest surge level increases from 10% to 20%. Storm events with moderate surge levels (0.8-1.4 m) would also affect the estimation of the extreme surge levels, but the uncertainty is generally less than 4%. Thus, specific attention should be paid to typhoon events with higher surge levels, which dominate the estimation of extreme surge levels.

Conclusions
In this study, we conducted a high-resolution wave-circulation coupled model for the Changjiang River Estuary and the continental shelves. The coupled model is forced by tidal forcing at the open boundary and atmospheric forcing at the surface. The model performance in reproducing sea levels and waves was assessed with in situ observations. Two different atmospheric products, the ERA reanalysis wind and the ERA-based synthetic wind, were used to estimate the surge levels and to discern the potential uncertainties associated with winds. Comparison with observed surge levels show that both winds reproduce the storm surge levels reasonably well. The synthetic winds, which were considered to produce higher surge levels, underestimate the surge levels in

Conclusions
In this study, we conducted a high-resolution wave-circulation coupled model for the Changjiang River Estuary and the continental shelves. The coupled model is forced by tidal forcing at the open boundary and atmospheric forcing at the surface. The model performance in reproducing sea levels and waves was assessed with in situ observations. Two different atmospheric products, the ERA reanalysis wind and the ERA-based synthetic wind, were used to estimate the surge levels and to discern the potential uncertainties associated with winds. Comparison with observed surge levels show that both winds reproduce the storm surge levels reasonably well. The synthetic winds, which were considered to produce higher surge levels, underestimate the surge levels in some places. This is likely due to the less realistic representation of typhoon structure in the synthetic wind, which consists of a theoretical typhoon wind. In contrast, the simulation forced by the ERA reanalysis shows comparative skills with the synthetic winds and is adopted for estimations of extreme surge levels.
There were 98 typhoon storms affecting the Changjiang River Estuary from 1986 to 2018, and are generally divided into four categories, with 47% having passed the East China Sea without landfall, 41.8% making landfall south of the Changjiang River Estuary, 4.1% making landfall right around the Changjiang River Estuary and 7.1% making landfall shortly and then turning northeast. Among them, typhoons making landfall south of the Changjiang River Estuary have the largest surges over the Changjiang River Estuary, and typhoons passing through the East China Sea without landfall are in second place. Historical higher surge levels are mainly located in Hangzhou Bay and coasts along Zhejiang province. Typhoons passing through the East China Sea without landfall produce higher surge levels north of the Changjiang River Estuary than others. Based on Gumbel distribution, the extreme surge levels with a 50-year return period were estimated. The extreme surge levels are about 1.0 m to 3.5 m in Hangzhou Bay and coasts along Zhejiang, and are about 1.0-2.0 m over the Changjiang River Estuary and the northern coasts.
Comparative runs with and without wave effects were conducted to discern the effects of waves on the extreme surge levels. Wave setup of about 0.05-0.25 m was found in the nearshore inner shelf with depth shallower than 20 m. The predominant wave setup is located along the coast of Zhejiang province. Significant wave setup was also found right outside the Changjiang River Estuary. The wave effects in the Hangzhou Bay are relatively weak due to the block of multiple Zhoushan islands. The wave setup contributes to 2-12.5% of the 50-year extreme surge levels along the east China coast. It is clear from this study that the wave-current interaction will not only affect the temporary sea levels, but also the estimation of long-term extreme surge levels.
The joint exceedance probabilities of high surge levels and high wave heights were evaluated with the Gumbel-logistic statistic model. Given the joint return period, high surges and large waves occur synchronously along the coasts of Zhejiang, while over the Changjiang River Estuary, large waves are accompanied by moderate surges. Thus, given the same hazard levels, the conditional return periods are shorter along the coasts of Zhejiang. These results indicate that the joint probability analysis can contribute meaningfully to hydrological engineering design and management.