Assimilation Research of Wind Stress Drag Coefficient Based on the Linear Expression

Abstract

The wind stress drag coefficient plays an important role in storm surge models. This study reveals the influences of wind stress drag coefficients, which are given in form of formulas and inverted by the data assimilation method, on the storm surge levels in the Bohai Sea, Yellow Sea, and East China Sea during Typhoon 7008. In the process of data assimilation, the drag coefficient is based on the linear expression C_d = (a + b × U₁₀) × 10⁻³ (generally speaking, a and b are empirical parameters determined by observed data). The results showed that the performance of the data assimilation method was far superior to those of drag coefficient formulas. Additionally, the simulated storm surge levels obviously changed in the neighborhood of typhoon eye. Furthermore, the effect of initial values of a and b in the C_d expression on the storm surge levels was also investigated when employing the data assimilation method. The results indicated that the simulation of storm surge level was the closest to the observation when a and b were simultaneously equal to zero, whereas the simulations had slight differences when the initial values of a and b were separately equal to the drag coefficients from the work of Smith, Wu, and Geernaert et al.. Therefore, we should choose appropriate initial values for a and b by using the data assimilation method. As a whole, the data assimilation method is much better than drag coefficient parameterization formulas in the simulation of storm surges.

1. Introduction

Storm surges are anomalous changes of sea water level induced by typhoons or tropical cyclones. The southeast coastal regions in China often suffer from tropical cyclones of the Northwestern Pacific Ocean. Such substantial destructive typhoons are able to cause severe economic loss and threaten life safety in the low-lying areas along the coast, especially in areas where typhoons pass through [1,2,3,4,5,6]. Chen et al. used a fully coupled tide–surge–wave model and three kinds of wind products to evaluate the effect of wind forcing on the hindcasting of typhoon-driven extreme waves and storm surges on the northeastern coast of Taiwan [5]. The results showed that the Cross-Calibrated Multi-Platform V2 (CCMPV2) combined with the parametric typhoon model had the best overall performance. Based on a coupled wave-circulation model, Hisao et al. computed significant wave heights by introducing a combination of wind field datasets and super Typhoon Nepartak (2016), and they investigated the model’s performance under varying spatial and temporal resolutions [6]. In 2020, the economic loss caused by storm surge disasters accounted for 97% of the total economic loss resulted from all kinds of oceanic disasters (Bulletin of China Marine Disaster 2020, http://www.mnr.gov.cn/sj/sjfw/hy/gbgg/zghyzhgb/ accessed on 14 August 2021). Therefore, it is necessary to improve storm surge models and the prediction of sea level for the sake of preventing and responding to storm surge disasters in advance.

As an important parameter in storm surge models, the wind stress drag coefficient has been given much attention. In early studies, the drag coefficient appeared in the form of linear expression [7,8,9,10]. In recent years, many scholars have been devoted to optimizing and improving the drag coefficient through different methods [11,12,13,14,15,16,17,18,19,20]. By combining observed data, Guan and Xie extended the drag coefficient to a near-linear function that depended on wave age [11]. Jarosz et al. evaluated the momentum transfer across the air–sea interface in the form of the drag coefficient; their results showed that the drag coefficient increased at first and then peaked at about 32 m/s for wind speeds between 20 and 48 m/s [12]. Moon et al. investigated the effects of spatial resolution on storm surge modeling and found that the RMS errors between the storm surge levels calculated by three wind stress drag coefficients and observations were totally different for three kinds of spatial resolutions. There was no uniform conclusion regarding the influence of spatial resolution on storm surges when using different drag coefficients [13]. Zhao et al. proposed a drag coefficient formulation that depended on wind speed and water depth during typhoon landfall in the South China Sea. The formulation of the drag coefficient was also applied to a typhoon forecast model, which suggested that predictions of surface wind speed and typhoon track were promoted [14]. Based on numerical experiments and observations in the South China Sea, Cao et al. formulated the drag coefficient related to wind speed in the form of a piecewise function [15]. This parameterization of the drag coefficient raised the model accuracy compared to that from Wu [2] used in SWAN. Zou et al. estimated the drag coefficient at high wind speeds using the turbulence closure and bulk model according to the observed ocean current and temperature profiles during Typhoon Megi, indicating that the drag coefficient increased to the peak at about 30 m/s and then declined as wind speeds increased [18]. When the inlets were closed and the effect of tide propagation was nullified, Mel et al. found that the action of wind forcing on a cross-lagoon setup was enhanced and became more clearly recognizable, allowing for a robust calibration of the wind stress drag coefficient with particular linear relationships [20]. They also achieved a robust calibration of the wind drag coefficient for moderate wind speed in the closed lagoon. To aid us in obtaining an optimal estimate of drag coefficient, they provided a recent overview of the linear relationship between the drag coefficient and wind speed. These researchers optimized the earlier drag coefficient expressions and even provided the drag coefficient parameterizations at high wind speeds.

Data assimilation is a significant tool to realize the combination of numerical models and observed data. In order to improve the accuracy of prediction in storm surge models, the data assimilation method has been widely used in recent years [21,22,23,24,25,26,27]. Peng and Xie constructed a four-dimensional variational (4D-Var) assimilation algorithm according to linear and adjoint three-dimensional Princeton Ocean Models, and they evaluated the effect of initial conditions on storm surge forecasting [22]. Using this adjoint optimal technique of 4D-Var, Peng et al. continued to adjust the initial and surface boundary conditions in the process of data assimilation, suggesting that it was better to adjust both conditions at the same time because of uncertain errors [23]. To study the linear and nonlinear bottom friction parameterizations in regional ocean tidal model, Zhang et al. simulated the M2 tide in the Bohai Sea and the Yellow Sea using the data assimilation method [25]. Zheng et al. adopted the adjoint-free data assimilation method to discuss the wind drag coefficient in the German Bight [27]; their results showed that the storm surge forecasting was more precise when the drag coefficient was adjusted by employing the data assimilation method. In addition, Flowerdew et al. developed and evaluated the Met Office Global and Regional Ensemble Prediction System, providing a range of atmospheric evolutions that were consistent with the latest observations [28]. Statistical verification and probabilistic verification ensured that ensemble forecasting system is applicable to simulating storm surges in coastal regions.

In the present study, we compared the influence of six wind stress drag coefficients on a storm surge model in the Bohai Sea, Yellow Sea, and East China Sea during Typhoon 7008. Here, the first five coefficients were calculated with the drag coefficient formulas provided by previous academics, and the last was the drag coefficient of spatial distribution inverted by the data assimilation method in this paper. We also discuss the influence of initial values of a and b in the C_d expression on storm surge levels in the process of data assimilation.

The paper is composed of four sections. The numerical model and adjoint assimilation model for storm surge, as well as the drag coefficient expressions and experimental designs, are introduced in the second section. In the third section, the influences of the drag coefficients and initial values of a and b in the C_d expression on the storm surge levels are discussed in detail. The conclusion and discussion are provided in the final section.

2. Materials and Methods

2.1. Typhoon and Stations

Typhoon 7008 was generated from the surface of the West Pacific on 21 August 1970. After reaching strong typhoon status, Typhoon 7008 began to travel northwest and entered the East China Sea and Yellow Sea. At about 14:00 on August 31, it weakened into a tropical storm and made landfall on the southwest region of North Korea with a maximum wind speed of 25 m/s and a minimum central pressure of 985 hPa. Finally, Typhoon 7008 vanished in the form of tropical depression on September 1 in North Korea. In the present study, the selected time periods and trajectory of Typhoon 7008 are shown in Figure 1, where the track and intensity data came from Typhoon Online (“http://www.typhoon.org.cn/” accessed on 14 August 2021). The distribution of the 9 tidal stations used in the model is also shown in Figure 1, and the longitude and latitude of these stations are given in

Table 1. Tidal stations used in the model.

Station	Lon (°E)	Lat (°N)
YingKou	122.15	40.63
HuLuDao	120.67	40.72
QinHuangDao	119.62	39.92
LongKou	120.32	37.65
YanTai	121.38	37.55
RuShan	121.48	36.80
QingDao	120.30	36.08
ShiJiuSuo	119.55	35.38
LianYunGang	119.45	34.75

.

2.2. Numerical Model and Adjoint Assimilation Model

In this study, the storm surge was simulated by the depth-averaged numerical model in two-dimensional form. With a right-handed Cartesian coordinate system, the governing equations of this numerical model consist of a vertically integrated continuity equation and momentum equations, as shown below [24,29].

$$\frac{\partial \zeta }{\partial t}+\frac{\partial \left[\left(h+\zeta \right)u\right]}{\partial x}+\frac{\partial \left[\left(h+\zeta \right)v\right]}{\partial y}=0$$

(1)

$$\begin{array}{c}\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}-fv+\frac{ku\sqrt{u^2+v^2}}{h+\zeta }-A\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)\\ +g\frac{\partial \zeta }{\partial x}+\frac{1}{{\rho }_w}\frac{\partial P_a}{\partial x}-\frac{{\rho }_a}{{\rho }_w}\frac{C_d{W}_x\sqrt{W_x^2+W_y^2}}{h+\zeta }=0\end{array}$$

(2)

$$\begin{array}{c}\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+fu+\frac{kv\sqrt{u^2+v^2}}{h+\zeta }-A\left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}\right)\\ +g\frac{\partial \zeta }{\partial y}+\frac{1}{{\rho }_w}\frac{\partial P_a}{\partial y}-\frac{{\rho }_a}{{\rho }_w}\frac{C_d{W}_y\sqrt{W_x^2+W_y^2}}{h+\zeta }=0\end{array}$$

(3)

where:

t = time;

x and y = longitude and latitude, respectively;

h = unperturbed water depth;

ζ = sea surface level relative to the unperturbed depth;

u and v = the current speeds pointing east and north, respectively;

f = the Coriolis parameter;

k = the bottom friction factor;

A = the eddy viscosity coefficient in the horizontal direction;

g = the gravitational acceleration;

ρ_w = the sea water density;

ρ_a = the air density;

C_d = the wind-stress drag coefficient;

W_x and W_y = the surface wind field;

P_a = the air pressure on the sea surface.

Generally, the storm surge level can be calculated by the governing equations. However, the simulated water level from governing equations may not agree very well with the observed data because of the errors thar can result from the equations, parameters, approximate hypotheses, and so on. To reduce the errors between the simulation and observation, the Lagrange multiplier method was introduced.

In the present paper, the objective function was constructed in the following form:

$$J\left(\zeta \right)=\frac{1}{2}{K}_{\zeta }{\displaystyle {\int }_{\Sigma }{\left(\zeta -{\zeta }_{obs}\right)}^{2}dxdydt}$$

(4)

where:

K_ζ = constant;

ζ = simulation;

ζ_obs = observation.

Correspondingly, the Lagrangian function was constructed as follows:

$$\begin{array}{c}L=J\left(\zeta \right)+{\displaystyle {\int }_{\Sigma }{\zeta }_a\left\{\frac{\partial \zeta }{\partial t}+\frac{\partial \left[\left(h+\zeta \right)u\right]}{\partial x}+\frac{\partial \left[\left(h+\zeta \right)v\right]}{\partial y}\right\}dxdydt}\\ +{\displaystyle {\int }_{\Sigma }{u}_a\left[\frac{ku\sqrt{u^2+v^2}}{h+\zeta }-A\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)+\frac{1}{\rho }\frac{\partial p_a}{\partial x}-\frac{{\rho }_a}{\rho }\frac{C_d{W}_x\sqrt{W_x{}^2+W_y{}^2}}{h+\zeta }\right]dxdydt}\\ +{\displaystyle {\int }_{\Sigma }{v}_a\left[\frac{kv\sqrt{u^2+v^2}}{h+\zeta }-A\left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}\right)+\frac{1}{\rho }\frac{\partial p_a}{\partial x}-\frac{{\rho }_a}{\rho }\frac{C_d{W}_x\sqrt{W_x{}^2+W_y{}^2}}{h+\zeta }\right]dxdydt}\end{array}$$

(5)

Based on the methods of He et al. [30], the corresponding adjoint equations we obtained:

$$\begin{array}{c}\frac{\partial {\zeta }_a}{\partial t}+u\frac{\partial {\zeta }_a}{\partial x}+v\frac{\partial {\zeta }_a}{\partial y}+\frac{ku\sqrt{u^2+v^2}{u}_a}{{\left(h+\zeta \right)}^2}+\frac{kv\sqrt{u^2+v^2}{v}_a}{{\left(h+\zeta \right)}^2}\\ +g\frac{\partial u_a}{\partial x}+g\frac{\partial v_a}{\partial y}=K_{\zeta }\left(\zeta -\widehat{\zeta }\right)\end{array}$$

(6)

$$\begin{array}{c}\frac{\partial u_a}{\partial t}-\left[f+\frac{kuv}{\left(h+\zeta \right)\sqrt{u^2+v^2}}\right]v_a-u_a\frac{\partial u}{\partial x}-v_a\frac{\partial v}{\partial x}+\frac{\partial }{\partial x}\left(u{u}_a\right)\\ +\frac{\partial }{\partial y}\left(v{u}_a\right)+\left(h+\zeta \right)\frac{\partial {\zeta }_a}{\partial x}+A\left(\frac{{\partial }^2{u}_a}{\partial x^2}+\frac{{\partial }^2{u}_a}{\partial y^2}\right)-\frac{k\left(2u^2+v^2\right)u_a}{\left(h+\zeta \right)\sqrt{u^2+v^2}}=0\end{array}$$

(7)

$$\begin{array}{c}\frac{\partial v_a}{\partial t}-\left[f+\frac{kuv}{\left(h+\zeta \right)\sqrt{u^2+v^2}}\right]u_a-u_a\frac{\partial u}{\partial y}-v_a\frac{\partial v}{\partial y}+\frac{\partial }{\partial x}\left(u{v}_a\right)\\ +\frac{\partial }{\partial y}\left(v{v}_a\right)+\left(h+\zeta \right)\frac{\partial {\zeta }_a}{\partial y}+A\left(\frac{{\partial }^2{v}_a}{\partial x^2}+\frac{{\partial }^2{v}_a}{\partial y^2}\right)-\frac{k\left(u^2+2v^2\right)v_a}{\left(h+\zeta \right)\sqrt{u^2+v^2}}=0\end{array}$$

(8)

where:

${\zeta }_a$, $u_a$ and $v_a$ = the adjoint variables of $\zeta$, $u$, and $v$, respectively.

The wind field and pressure field of typhoon adopted in this model were from the work of Jelesnianski [31]. The expression of the former is as follows:

$$\overrightarrow{W}=\{\begin{array}{c}\frac{r}{R+r}\left(V_{ox}\overrightarrow{i}+V_{oy}\overrightarrow{j}\right)+W_R\frac{1}{r}{\left(\frac{r}{R}\right)}^{\frac{3}{2}}\left(A\overrightarrow{i}+B\overrightarrow{j}\right),\begin{array}{cc}& 0<r\le R\end{array}\\ \frac{R}{R+r}\left(V_{ox}\overrightarrow{i}+V_{oy}\overrightarrow{j}\right)+W_R\frac{1}{r}{\left(\frac{R}{r}\right)}^{\frac{1}{2}}\left(A\overrightarrow{i}+B\overrightarrow{j}\right),\begin{array}{cc}& r>R\end{array}\end{array}$$

(9)

where:

$\overrightarrow{i}$ and $\overrightarrow{j}$ = the unit vector points of the x and y axes, respectively;

$V_{ox}$ and $V_{oy}$ = migration velocities of typhoon center;

$r$ = the distance between the grid center $\left(x,y\right)$ and typhoon center $\left(x_c,y_c\right)$;

R = the radius of maximum wind speed $W_R$;

$A=-\left[\left(x-x_c\right)\sin \theta +\left(y-y_c\right)\cos \theta \right]$;

$B=\left[\left(x-x_c\right)\cos \theta -\left(y-y_c\right)\sin \theta \right]$;

$\theta =\{\begin{array}{c}{20}^{\circ },\begin{array}{cc}& r\le R\end{array}\\ {15}^{\circ },\begin{array}{cc}& r>R\end{array}\end{array}$.

The latter (i.e., pressure field) is:

$$P_a=\{\begin{array}{c}{P}_0+\frac{1}{4}\left(P_{\infty }-P_0\right){\left(\frac{r}{R}\right)}^3,\begin{array}{cc}& r\le R\end{array}\\ P_{\infty }-\frac{3}{4}\left(P_{\infty }-P_0\right)\left(\frac{R}{r}\right),\begin{array}{cc}& r>R\end{array}\end{array}$$

(10)

where:

P_a = the pressure at r on the sea;

P₀ = pressure at the typhoon center;

P_∞ = ambient pressure.

2.3. Wind-Stress Drag Coefficient

In a storm surge model, the drag coefficient C_d is generally used to calculate the wind stress. The parameterization of C_d has been widely studied through laboratory experiments, model simulations, and oceanic measurements [7,8,9,11,12,13,15,18,32,33]. Most previous studies exhibited a linear increase of the drag coefficient with wind speed U₁₀ at 10 m above the sea surface. However, some quadratic expressions for the drag coefficient were developed in recent years, and in this expression, C_d achieves its maximum when the wind speed researches a certain value.

The five drag coefficient formulas [7,8,9,15,32] used in this paper are listed below.

$$C_d=\left(0.61+0.063\times U_{10}\right)\times {10}^{-3},\begin{array}{cc}& 6\le U_{10}\le 22\end{array}\mathrm{m}/\mathrm{s},$$

(11)

$$C_d=\left(0.8+0.065\times U_{10}\right)\times {10}^{-3},\begin{array}{cc}& 0\le U_{10}\le 50\end{array}\mathrm{m}/\mathrm{s},$$

(12)

$$C_d=\left(0.577+0.085\times U_{10}\right)\times {10}^{-3},\begin{array}{cc}& 4\le U_{10}\le 24\end{array}\mathrm{m}/\mathrm{s},$$

(13)

$$C_d=\{\begin{array}{c}\left(0.29+3.1/U_{10}+7.7/U_{10}^2\right)\times {10}^{-3},\begin{array}{cc}& 3\le U_{10}\le 6\end{array}\mathrm{m}/\mathrm{s}\\ \left(0.6+0.07\times U_{10}\right)\times {10}^{-3},\begin{array}{cc}& 6\le U_{10}\le 26\end{array}\mathrm{m}/\mathrm{s}\end{array},$$

(14)

$$C_d=\{\begin{array}{c}1.06\times {10}^{-3},\begin{array}{cc}& U_{10}\le 8\end{array}\mathrm{m}/\mathrm{s}\\ \left(0.396+0.083\times U_{10}\right)\times {10}^{-3},\begin{array}{cc}& 8\le U_{10}\le 24\end{array}\mathrm{m}/\mathrm{s}\\ \left(2.1+0.012\times U_{10}\right)\times {10}^{-3},\begin{array}{cc}& 24\le U_{10}\le 34\end{array}\mathrm{m}/\mathrm{s}\end{array}.$$

(15)

Here, the expressions of C_d from Smith [7], Wu [8], and Geernaert [9] have a simple linear relationship with the wind speed U₁₀. When the wind speed is larger than 6 m/s, the piecewise function C_d of Yelland and Taylor [32] also supports a linear increase with wind speed U₁₀. Analogously, the latter two parts from piecewise Function (15) are in linear forms, and the only distinction is the unequal coefficients. These five drag coefficients calculated from Equations (11)–(15) are shown in Figure 2.

Considering that most of drag coefficients present a linear relationship with wind speed U₁₀, we also adopted a linear expression in the numerical model, that is, C_d = (a + b × U₁₀) × 10⁻³, where a and b are empirical parameters determined by observations [11,20]. Then, coefficients a and b were inverted by the data assimilation method with Lagrangian multiplier. Finally, we provided the spatial distribution of the drag coefficient in certain seas. We also evaluated the effect of different initial values of a and b in the C_d expression on the storm surge model when employing the data assimilation method.

2.4. Experimental Designs and Model Setting

2.4.1. Experimental Designs

In order to achieve the goals of the present study, two sets of numerical experiments were designed. In the first experiment, five drag coefficient formulas from Equations (11)–(15) were used for the storm surge model (denoted by E1, E2, E3, E4, and E5, respectively) in order to evaluate the influence of the drag coefficient on the storm surge level (

Table 2. Experiments evaluating the influence of the drag coefficient on storm surge level.

Experiments	Formula of Cd
E1	Smith (1980)
E2	Wu (1982)
E3	Geernaert (1987)
E4	Yelland and Taylor (1996)
E5	Cao et al. (2017)

). In the second experiment, we estimated the influence of the initial values of a and b in the C_d expression on the storm surge level in the process of data assimilation. For this experiment, four initial values of a and b denoted by E6, E7, E8, and E9 were selected. In the E6, the initial values of a and b were equal to 0, and the remaining three initial values were the same as those from Equations (11)–(13) (

Table 3. Experimental results of the estimation of the influence of the initial values of a and b in the C_d expression on the storm surge level.

Experiments	a	b
E6	0	0
E7	0.61	0.063
E8	0.8	0.065
E9	0.577	0.085

).

2.4.2. Model Setting

In the present paper, the study area of the storm surge simulation was 24–44 °N and 116–130.5 °E with a 5′ × 5′ grid resolution, which covered the Bohai Sea, Yellow Sea, and East China Sea. The model was run with a cold start, the current velocity and sea surface level taken as 0, and a hypothesis of a constant bottom drag coefficient (0.0016) all seas. The open boundary conditions in the model were the Taiwan Strait and the first island chains. Along the closed boundaries, we assumed that no water flowed into or away from the coast. In this model, a staggered grid (Arakawa C grid) was used, and the wind stress and pressure were the driving forces. The time step was 60 s. To obtain the temporal variation of the drag coefficient, the migration process of Typhoon 7008 was separated into 12 periods, and every period lasted for 6 h. The observed data used for assimilation were from 9 tidal stations during Typhoon 7008.

3. Results

3.1. Drag Coefficients during Typhoon 7008

The drag coefficients calculated by the five C_d formulas from Equations (11)–(15) and the data assimilation method are compared in this section. Figure 3 shows the spatial distributions of the computed drag coefficients at 08:00 am on 31 August 1970 calculated in the six experiments E1–E6. Figure 3 shows both the magnitudes of C_d and their spatial distributions with the positions of the extreme values C_d during the typhoon. The calculated minimum C_d in E1–E5 (Figure 3a–e) was near the typhoon center, and the maximum was around the periphery of the typhoon eye. Finally, C_d decreased gradually from the middle to the margin. This spatial distribution of C_d was principally related to the adopted circular wind field from Jelesnianski [31], whereas the C_d calculated using the data assimilation method was obviously different from those in E1–E5. Figure 3f shows that the extrema of C_d mainly appeared in the Bohai Sea and the Yellow Sea, particularly in the coastal areas.

Different spatial distributions of drag coefficients can influence the estimation of wind stress and storm surge level. In next section, the comparison of time series between the simulated storm surge levels in E1–E6 and observations is discussed.

3.2. Influence of the Drag Coefficient on the Storm Surge Level

The influence of the drag coefficient on storm surge level was studied on the basis of experiments E1–E6. We used five drag coefficient formulas from Equations (11)–(15) (i.e., E1–E5) to calculate the storm surge levels during Typhoon 7008. Nevertheless, in experiment E6, the data assimilation method was applied and the initial values of a and b in the C_d expression were set to 0.

Table 4. Root mean square errors between the simulated storm surge levels and observations in every period of Typhoon 7008 (unit: cm).

Exp	1	2	3	4	5	6	7	8	9	10	11	12	Mean
E1	23	24	37	37	42	38	35	20	54	90	105	103	51
E2	25	27	39	42	49	43	38	22	58	97	114	112	55
E3	25	27	40	43	51	43	39	26	61	101	118	114	57
E4	24	25	38	39	45	40	37	24	56	94	110	107	53
E5	24	25	38	39	45	40	37	24	57	94	110	106	53
E6	11	10	12	14	13	15	14	10	8	7	5	6	10

lists the root mean square (RMS) errors between the simulated storm surge levels and observations in the 12 periods during Typhoon 7008. The mean RMS errors in E1–E5 were 51, 55, 57, 53, and 53, respectively. The E1 result from Formula (11) of Smith [7] was the best, and the E3 result from the formula of Geernaert [9] was the worst. However, after comparing the results in E1–E5 to that in E6, the superiority of the data assimilation method becomes apparent. The mean RMS error in E6 was only 10 cm, far less than the errors in E1–E5.

Table 5. Root mean square errors between the simulated storm surge levels and observations at 9 tidal stations during Typhoon 7008 (unit: cm).

Tidal Stations	E1	E2	E3	E4	E5	E6
YingKou	85	93	95	89	89	7
HuLuDao	80	87	89	83	83	11
QinHuangDao	70	76	77	73	73	15
LongKou	58	63	64	60	60	5
YanTai	49	53	54	51	51	6
RuShan	40	44	46	42	42	9
QingDao	37	41	43	39	39	7
ShiJiuSuo	40	44	46	42	42	9
LianYunGang	46	53	55	50	50	21
Mean	56	62	63	59	59	10

also shows the similar results at the nine tidal stations, that is, the performance of the data assimilation method was much better than the performance of the drag coefficient formulas.

Figure 4 demonstrates the spatial distributions of the simulated storm surge levels in E1–E6. It can be seen that the models with similar drag coefficients (Figure 3a–e) produced similar storm surge levels (Figure 4a–e), while the surge levels (Figure 4f) calculated by the model with the data assimilation method were entirely different from those shown in Figure 4a–e. Additionally, the storm surge levels obviously changed in the neighborhood of the typhoon eye but varied slightly in the open sea.

The comparison of time series and differences between the simulated storm surge levels and observations at the YingKou, LongKou, and YanTai stations (see Figure 5, Figure 6 and Figure 7) revealed that the influences of the drag coefficients on the storm surge levels were very distinct. Here, the results in E1–E5 obviously departed from the observations. Only the simulated levels in E6 were close to the observed levels. Additionally, the surge peak results from E1, E3, and E6 at the YingKou, LongKou, and YanTai stations were compared to the observed surge peak (see Figure 8). The formula of Geernaert [9] overrated the peaks of surges by about 23.37 cm at YingKou station, 69.59 cm at LongKou station, and 27.42 cm at YanTai station, but the overestimations in E6 were only 3.06, 1.58, and 11.35 cm at the same three tidal stations, respectively.

3.3. Influence of the Initial Values of a and b on the Storm Surge Level

In the present section, we discuss the influence of the different initial values of a and b on C_d expression in the storm surge model based on the data assimilation method. The initial values of a and b in E6–E9 are listed in Table 3. The RMS errors between the simulated storm surge levels and observations in the 12 periods (

Table 6. Root mean square errors between the simulation and observation of storm surge levels in every period of Typhoon 7008 (unit: cm).

Exp	1	2	3	4	5	6	7	8	9	10	11	12	Mean
E6	11	10	12	14	13	15	14	10	8	7	5	6	10
E7	11	11	13	13	14	10	9	8	13	26	31	31	16
E8	11	11	14	13	14	12	10	8	16	29	34	35	17
E9	12	12	15	13	16	14	12	8	16	30	37	37	18

) at the nine tidal stations (

Table 7. Root mean square errors between the simulation and observation of storm surge levels at 9 tidal stations during Typhoon 7008 (unit: cm).

Tidal Stations	E6	E7	E8	E9
YingKou	7	11	11	13
HuLuDao	11	20	22	23
QinHuangDao	15	31	35	37
LongKou	5	18	20	22
YanTai	6	16	18	19
RuShan	9	9	9	9
QingDao	7	10	11	11
ShiJiuSuo	9	12	13	13
LianYunGang	21	22	24	25
Mean	10	17	18	19

) revealed that the performance of the model in E6 was the best. Figure 9, Figure 10 and Figure 11 show the time series of simulated and observed storm surge levels, as well as their differences at YingKou, LongKou, and YanTai stations. This shows that the results in E6 were the closest to the observations. Comparisons of surge peaks among E6–E9 (see Figure 12) also exhibited analogous results. Therefore, the influence of initial values of a and b (zero or non-zero initial values) on the storm surge model was evident in the process of data assimilation. However, regardless of the initial values of a and b, the data assimilation method was superior to the drag coefficient formulas.

4. Discussion and Conclusions

The wind stress drag coefficient is a very important parameter in storm surge models. The purpose of this paper was to evaluate the influence of the different drag coefficients on storm surge levels and the effect of selecting initial values for a and b in the C_d expression on simulated storm surge levels when using the data assimilation method in a storm surge model.

In order to study the first issue, a series of experiments were designed for six drag coefficient parameterizations: (1) the linearly increasing drag coefficient formulas used by Smith [7], Wu [8], and Geernaert [9]; (2) the piecewise functions used by Yelland and Taylor [32] and Cao et al. [15]; and (3) the drag coefficient inverted by the data assimilation method. The experimental results for Typhoon 7008 illustrated that the performance of the data assimilation method was much better than that of the drag coefficient formulas. Hence, the superiority of the data assimilation method is obvious. In addition, the storm surge levels significantly varied around the typhoon eye but had no evident changes in the open sea.

The influence of the data assimilation method on a storm surge model depends on the selection of the initial values of a and b in the C_d expression. Here, when a and b were equal to zero, the performance of the model was the best. When a and b were taken from the coefficients of Smith [7], Wu [8], and Geernaert [9], the root mean square (RMS) errors between the simulations and observations of the storm surge levels presented small differences. Therefore, when using the data assimilation method, one should pay attention to the option of initial values of a and b. Regardless of whether the initial values of a and b are zero, the data assimilation method is superior to drag coefficient parameterization formulas.

In the present study, we designed these experiments for Typhoon 7008 and obtained a series of results. For other typhoons, the influences of drag coefficients and initial a and b values’ selection on storm surge models may be different. In future study, we will continue to explore related issues.