1 INTRODUCTION

The Bohai Sea, Yellow Sea, and East China Sea (BYECS) are shallow, except for the narrow Okinawa Trough near the Ryukyu Islands (Fig. 1). As the waters around and east of those islands are deep with rough topography, tidal waves experience strong internal tide dissipation when propagating from the Pacific toward the BYECS. Since satellite altimeter data have become available, knowledge of the distribution of tides in the BYECS has been greatly improved (Yanagi et al., 1997). Using 10 years of TOPEX/ Poseidon (TP) data, Fang et al. (2004) showed that the accuracy of their TP solutions in the BYECS was within 2-4 cm in amplitude and 5° in phase-lag for principal constituents. Numerical studies by a number of researchers have greatly improved knowledge of tidal dynamics in this area (Choi, 1980; Zhao et al., 1993; Wang et al., 1999; Lü and Fang, 2002; Li et al., 2005; Lü and Zhang, 2006; Han et al., 2006; Zhu et al., 2012, 2014). However, in previous numerical models, internal tide dissipation and effects of the self-attraction and loading (SAL) tide have not been considered.

Egbert and Ray (2000) showed that ~1 TW of tidal energy (~1/4 of global total tidal energy dissipation) is dissipated in the deep ocean through conversion from surface tides to internal tides. The East China Sea contains a deep water area in its southeast and thus tidal waves there should experience significant dissipation through barotropic to baroclinic tidal conversion. Niwa and Hibiya (2004) showed that 35 GW of M₂ surface tide energy is converted to internal tide energy near the Ryukyu Islands ridges and Luzon Strait. Zu et al. (2008) also showed considerable energy dissipation in the Luzon Strait, associated with the generation of internal tides. Direct simulation of internal tides requires a threedimensional model with very fine resolution. When a two-dimensional barotropic model is used, a parameterized term should be included to indicate energy loss from the conversion of surface to internal tides. In addition, it is well known that SAL tides are important in tidal modeling (e.g., Ray, 1998). For the BYECS, Fang et al. (2013) noticed that although amplitudes of semidiurnal SAL tides are smaller than those of equilibrium tides, because of shorter wavelengths, gradients of the former tides are greater and thus may have a stronger influence on ocean tides. The SAL tide has not been considered in numerical simulations of this area, although the equilibrium tide has often been treated.

In the present paper, we extend the study to include the internal tidal dissipation and SAL tide terms using the two-dimensional version of the Finite-Volume Community Ocean Model (FVCOM) to simulate constituents M₂ and S₂ simultaneously. Diurnal tides are not simulated together with semidiurnal tides because diurnal internal tidal waves behave very differently from semidiurnal ones; their dissipation parameters should be different, at least in principle. A devoted effort will be made in the future to simulate additional tidal constituents.

2 NUMERICAL MODEL AND ITS IMPLEMENTATION 2.1 Governing equations

We use the vertically integrated momentum and continuity equations, inclusive of internal tide dissipation and SAL tide as governing equations. These are written in spherical coordinates as follows.

where u, v are zonal and meridional components of the velocity, t is time, R is the Earth's radius, a=Rcosϕ, λ is longitude, ϕ is latitude, ζ is height of the free surface, ζ and are adjusted heights of the equilibrium tide (representing the tide-generating force) and SAL tide, respectively, and g is the acceleration due to gravity. H=h+ζ is total water column depth with h representing the undisturbed water depth, r is the bottom friction coeffi cient. A is the horizontal eddy viscosity calculated via the Smagorinsky formulation (Chen et al., 2006).

In Eqs. 2 and 3, is the parameterized term for internal tide dissipation in the form proposed by Jayne and St Laurent (2001), in which (κ, h) are the wave number and amplitude scales for topographic roughness. This form of parameterization has been successfully applied to the simulation of global ocean tides and South China Sea tides by Jayne and St Laurent (2001) and Gao et al. (2015), respectively. However, as pointed out by Jayne and St. Laurent in that work, in principle, the parameterization form, or at least the coefficient, should be dependent on tidal frequency. Thus, the present study only treats the semidiurnal tides. As in Jayne and St Laurent (2001), we regard κ as an adjustable parameter to be determined from numerical experiments. The value of h² at depth grid point (i, j) is calculated from the following equation, as done in Gao et al. (2015).

is the buoyancy frequency, with ρ and ρ₀ the water and reference densities. N_b can be calculated using the WOA13 temperature and salinity dataset (Locarnini et al., 2013). Because the internal tides are generated near the seabed, ∂ρ/∂z is taken as the vertical average of the lowest fi ve layers.

Equilibrium tidal heights adjusted for Earth's elasticity are expressed in the form

where subscript C stands for the tidal constituent with C=1 for M₂ and C=2 for S₂, ω is angular speed of the tidal constituents, f is the nodal factor, u is the nodal angle, V₀ is initial phase angle, and p=2 for semidiurnal tides. According to Wahr (1981), amplitudes H_C for M₂ and S₂ are equal to the following (in m):

The SAL tide is given as

where are amplitudes and Greenwich phase-lags of M₂ and S₂ SAL tides, respectively. Readers may refer to Pugh (1987) for further details on the adjusted equilibrium and the SAL tides.

2.2 Model setting and implementation

We used the FVCOM, which uses an unstructured grid, for numerical computation. The model domain was 117°-140°E and 18°-41°N. The finest horizontal resolution was 0.04°, as shown in Fig. 2. Along the closed boundary, the normal velocity was set to zero. Along the open boundary, tidal heights were calculated from the harmonic constants of M₂ and S₂, using

where H_C and G_C are harmonic constants obtained from global tidal model DTU10. Harmonic constants of the SAL tides were acquired from Fang et al. (2013) (see Fig.B1 of Appendix for their distribution). The initial condition was ζ=u=v=0 when t=0. The model was run for 30 days, and results of the last 15 days were used for the harmonic analysis.

3 NUMERICAL SIMULATION

In the numerical simulation, the terms in Eq.6 were all included, and the parameters in the equation were obtained through an optimization procedure.

3.1 Parameter optimization

In the study area, we selected 59 stations where observed tidal harmonics were available for parameter optimization and model verification. The locations of these stations are shown in Fig. 3. Along the coast of the China mainland, west coast of the Korean Peninsula, and Ryukyu Islands, the harmonic constant was taken from the database of the International Hydrographic Office. For the Taiwan Strait, the harmonic constants were taken from Jan et al. (2002). For the off shore area, the harmonic constants were obtained by analyzing TOPEX/Poseidon altimetry at the crossover points (Wang et al., 2012).

We used the RMS difference between model outputs and observations at 59 stations (Fig. 3) as the cost function to optimize model parameters and evaluate accuracy of the model results. The cost function was calculated by

where F_C is the RMS vector difference between computed and observed harmonics of the constituent C, taken as

where k is the station index, H and G are harmonic constants for amplitude and phase-lag, and Com and Mea are computed and observed results, respectively. The model produces tidal harmonic constants by T_tide (Pawlowicz et al., 2002).

We assumed that the bottom friction and internal tide dissipation coefficients are constants independent of time and space. Based on our previous experience, r was estimated to be ~0.001 5 for the study area, and κ was estimated to be of order of 10^-5/m. Thus, the optimal value of r should lie between 0.000 5 and 0.004, and that of κ between 0 and 0.000 1/m. To find the optimal (r, κ) combination, we first performed 48 numerical computations using coarser resolution. Comparing these 48 F values, we found that the minimum F appeared at (r, κ)=(0.001 5, 0.000 02). We took this point as the center, reduced the (r, κ) domain, increased the resolution, and performed a second set of computations. Those showed a minimum F at (r, κ)=(0.001 25, 0.000 03). Finally, we took this new point as the center, further reduced the domain of (r, κ), and executed a third set of computations. From that set, we found that when (r, κ)=(0.001 25, 0.000 035), the value of F was a minimum, 7.71 cm. The third set of computations also revealed no great change among F values at points neighboring that point, indicating that further refined computations were unnecessary. Therefore, we terminated the computation and took (r, κ)=(0.001 25, 0.000 035) as optimal estimates. Results of the three sets of computations are shown in Fig. 4, in which the units of κ are /m. The output obtained with the optimal (r, κ) combination is regarded as the optimal solution. For convenience, the simulation with all terms in Eq.6 included and parameters (r, κ)=(0.001 25, 0.000 035) is called the standard simulation in the following.

As stated above, the M₂ and S₂ mean RMS difference in the standard simulation was 7.71 cm. The mean absolute difference of amplitudes and phase-lags for M₂ were 6.36 cm and 5.27° respectively, and those for S₂ were 2.66 cm and 5.62°. A detailed station-by-station comparison of model-produced harmonic constants with observations is given in the Appendix.

3.2 Results of standard simulation

Cotidal charts for M₂ and S₂ based on the optimal solution are shown in Figs. 5 and 6, respectively. These charts agree well with observed ones from Fang (1986) and Fang et al. (2004), as well as with previous numerical products. This demonstrates that the model with optimized parameters is reliable. Characteristics of the semidiurnal tidal regimes have been ably described in previous works, and so are not repeated in the description here.

The energy flux across a section of unit width is called energy flux density (Fang et al., 1999). Vectors of energy flux density for M₂ and S₂ tides obtained are plotted in Figs. 7 and 8, respectively. These figures clearly show the strongest energy fluxes in passages of the Ryukyu Islands chain. Energy inflow into the East China Sea appears to have three main branches. The north branch is divided into two parts; most energy enters the Yellow Sea circumambulating Cheju Island and extends along the eastern shore, with a very small portion flowing toward the Japan/ East Sea, mainly along the Japanese coast. The middle branch flows northwestward and forms a counterclockwise rotating pattern with a southwest turn along the Shandong Peninsula, then extends southeastward and terminate at the middle North Jiangsu coast. The south branch first flows northward and then turns toward the Taiwan Strait, supporting strong tides there.

4 ENERGY BALANCE

To reveal the relative importance of each term in the governing equations, we analyzed the energy balance. The method used was similar to that of Fang et al. (1999). First, the study area was divided into five subareas as shown in Fig. 9. In this figure, G₁ through G₅ roughly represent the Bohai Sea, Yellow Sea, East China Sea shelf, East China Sea deep basin (including shelf margin and Okinawa Trough), and the east zone of the Ryukyu Islands. Boundaries separating these subareas are marked as B₁, B₂, B₅ and B₆. B₃ and B₄ are open boundaries of G₃ connecting the Taiwan Strait and Japan/East Sea, respectively. B₇ is the open boundary of G₅. G₁ to G₄ together represent the Bohai, Yellow and East China Seas, the area of focus of the present study. The analysis was conducted constituent by constituent. Because the governing equations contain some nonlinear terms, we first expanded these terms in Fourier forms:

Because model-produced hourly values of ζ, u and v are available, hourly values of each term on the left sides of the above equations are readily obtained. Then, the amplitudes and phase-lags of these terms (Z and γ in Eq.10) can be derived through harmonic analyses. The number of harmonic terms contained in Eqs.10-16 is larger than that in Eq.7, and we only used terms with frequencies equal to those of M₂ and S₂. For a specific constituent C, Eqs.1, 2 and 3 become correspondingly

Multiplying Eqs.17, 18 and 19 by ρgζ_C, ρgu_C and ρgv_C respectively and adding them, then taking the average for a period, and finally integrating over area G, we obtain

where subscript C has been dropped. The four terms on the left side represent energy input, and the three terms on the right side represent energy consumption. The energy flux entering G through its boundary is

where T is the period of the constituent; u_n is the velocity component normal to boundary B and is zero along the closed (solid) boundary, with B representing any of B_i (i=1, 2, …, 7) or their combinations. If Φ is calculated with reference to an area, the outward velocity (or inward flux) component is defi ned as positive. The rate of work done by equilibrium tides on an area G (G representing any of Gⁱ, i=1, 2, …, 5) can be calculated from

The rate of work by SAL tides can be calculated by

The values of N, D_b, D_i and D_l are calculated from

In comparison with the energy balance equation given by Fang et al. (1999), Eq.20 has two extra terms, i.e., the SAL term S and internal tide dissipation term D_i. Because the grid points of FVCOM are irregularly distributed, in practical calculation we first interpolate the model outputs to rectangular grid points, then calculate the integrands in the above equations at these points, finally integrating the integrands over subarea G. The interpolation can induce certain discrepancies, especially near the Ryukyu Islands chain, where islands are numerous and coastlines are very complicated. Thus, calculated values for energy terms in the following cannot completely satisfy the balance required by Eq.20.

The calculated energy fluxes F through boundaries B₁-B₇ are listed in Table 1. We see that tides in the Bohai, Yellow, and East China Seas are essentially maintained by energy fluxes from the Pacific Ocean through the Ryukyu Islands passages (B₆ in Fig. 9). The magnitudes of the fluxes (133.25 GW for M₂; 21.24 GW for S₂) agree well with the results of previous studies (see Zhu et al., 2014 and references therein). Energy fluxes through the Taiwan Strait (B₃) listed in Table 1 are southwestward and very large, implying strong southwestward semidiurnal waves propagating through the strait. In contrast, northeastward semidiurnal fluxes into the Japan/East Sea (B₄) are very small, consistent with weak tides in the sea.

Energy budgets for the five subareas and entire BYECS (G₁-G₄ total) are shown in Tables 2 and 3 for constituents M₂ and S₂, respectively. These tables show that maintenance of tidal motion in the study area almost entirely (exceeding 94% for both M₂ and S₂) relies on the energy flux from outer oceans. The tide-generating force contributes 1.2% M₂ power for the entire BYECS and up to 2.8% for the East China Sea deep basin. Strikingly, the SAL tide contributes much more power than the tide-generating force, 4.4% M₂ power for the BYECS and up to 9.3% for the East China Sea deep basin. The contribution of the advection term is generally negligible. The total dissipation D_b+D_i+D_l should be in balance with the total energy input Φ+E+S+N. However, Tables 2 and 3 show a small imbalance, which is mainly caused by model output interpolation, as mentioned above. Among the three dissipation terms, the contribution of lateral friction is very small. Bottom friction plays a major role in dissipating tidal energy in shelf regions G₁, G₂ and G₃. In contrast, the internal tide effect is important in deepwater regions G₄ and G₅. For the entire BYECS (G₁-G₄ total in Tables 2 and 3), the bottom friction still dissipates most tidal energy (~80% for both M₂ and S₂). Total internal tide dissipation for M₂ in the entire study area (G₁ to G₅ together) is estimated at 31 GW, similar to the 35 GW estimate of Niwa and Hibiya (2004) for a larger area. The internal tide dissipation is mostly near the Ryukyu Islands chain, which also agrees with Niwa and Hibiya.

5 NUMERICAL STUDY OF THE EFFECTS OF TIDE-GENERATING FORCE AND SAL TIDES

To further demonstrate the influences of the tidegenerating force and SAL tide on the ocean tide in the BYECS, we carried out two additional numerical experiments. In the first (Exp. 1), we artificially removed the tide-generating force in the BYECS (G₁ to G₄ in Fig. 9) by setting ζ=0 there, and retained the values in G₅ as in the standard simulation. The bottom friction, internal parameters and SAL tide were also unchanged. The difference between Exp. 1 and the standard simulation (vector difference of that simulation minus Exp. 1) for M₂ is shown in Fig. 10 (S₂ is not shown). In the second experiment (Exp. 2), we artificially removed the SAL tide in the BYECS (G₁ to G₄ in Fig. 9) by setting ₀ there, and retained the values in G₅ as in the standard simulation. The bottom friction, internal parameters and the tidegenerating force were unchanged. The difference between Exp. 2 and standard simulation (vector difference of that simulation minus Exp. 2) for M₂ is shown in Fig. 11.

Figure 10 reveals that the phase-lag distribution of the difference owing to removal of the tide-generating force in BYECS has a certain similarity to that of the cotidal chart (Fig. 5). However, the former lags the latter by ~160° in the southeast, decreasing to ~70° in the northern Yellow Sea, and ~30° in Liaodong Bay. The amplitude reaches its maximum, greater than 30 cm, in the Liaodong Bay. The figure gradually decreases further in Bohai Bay and northeastern North Yellow Sea (>25 cm), to Haizhou Bay (>20 cm) and finally Inchon Bay (~10 cm). Figure 11 shows that the phase-lag distribution of the difference caused by removal of the SAL tide in BYECS is very similar to that of the cotidal chart (Fig. 5). Also, their phase lags have no signifi cant difference. Maximum amplitude (>40 cm) appears in Liaodong Bay. The next largest (>35 cm) was in Bohai Bay, the northeastern North Yellow Sea, and Inchon Bay. Although the amplitude distribution (Fig. 11) of the difference between with and without SAL tide in the forcing term is very similar to that of ocean tide (Fig. 5), the relative magnitude (i.e., role of SAL tide) is spatially variable. For example, in Liaodong Bay, the form for M₂ is over 40 cm, amounting to ~1/3 of the latter; whereas in Inchon Bay that form is over 35 cm, only ~1/8 of the latter.

6 SUMMARY AND CONCLUSION

A parameterized internal tide dissipation term and SAL tide term were introduced for the fi rst time in a barotropic numerical model to investigate the dynamics of semidiurnal tidal constituents M₂ and S₂ in the BYECS. Optimal parameters for bottom friction and internal dissipation were obtained through a series of numerical computations. It was found that:

1. Maintenance of tidal motion in the study area almost entirely (> 94% for both M₂ and S₂) relies on the energy flux from outer oceans. The tide-generating force contributes 1.2% M₂ power for the entire BYECS and up to 2.8% for the East China Sea deep basin. The SAL tide contributes much more power than the tide-generating force, 4.4% M₂ power for the BYECS and up to 9.3% for the East China Sea deep basin (including shelf margin and Okinawa Trough);

2. Bottom friction plays a major role in dissipating tidal energy in the shelf regions, and the internal tide eff ect is important in the deepwater regions. For the entire BYECS, the bottom friction dissipates about 80% of tidal energy, and total internal tide dissipation accounts for the remaining 20%, for both M₂ and S₂;

3. Numerical experiments show that artificial removal of the tide-generating force in the BYECS can cause a signifi cant difference in model output. Maximum magnitude of the difference was in Liaodong Bay, exceeding 30 cm. The artificial removal of SAL tide in the BYECS can cause an even greater difference in model output. The maximum magnitude of the difference was also in Liaodong Bay, exceeding 40 cm. This clearly shows the importance of SAL tide inclusion in the numerical simulation of tides in the waters, even those with a greater proportion of shelf seas.

The main conclusions are as follows: a) internal tide dissipation is important in the shelf margin area and deep waters near the Ryukyu Islands, and can be parameterized by Jayne and St. Laurent's formula; b) the effect of the SAL tide is more important than the tide-generating force for ocean tides in the BYECS, so that tide should be accounted for in numerical simulations, especially if the tide-generating force is considered.

We used the simplest parameterization form for internal tide dissipation, proposed by Jayne and St Laurent (2001). Other forms, such as used by Green and Nycander (2013), should also be used in future tidal modeling of the BYECS. This is especially so for diurnal tides, which were not simulated in the present work.

7 ACKNOWLEDGEMENT

We are grateful to two anonymous reviewers for constructive comments and suggestions on an earlier version of this manuscript.