On the Simulation of Shallow Water Tides in the Vicinity of the Taiwan Banks

The Taiwan Banks (Formosa Shoals), a large NE-SW oriented bathymetric feature near the southern end (23°N, 118 119°E) of the Taiwan Strait, is a region of extremely shallow water that exerts a profound effect on the propagation of tidal waves. As such waves propagate over the Taiwan Banks, they become distorted and asymmetric due to bottom friction and contribute to the generation of shallow water tides. The POM model was used in present study to simulate the tides in the Taiwan Strait region. Shallow water tidal dynamics in the area of Taiwan Banks are focused. The numerical model was validated against sea level observations from 34 tidal stations located on the coast of Mainland China and Taiwan. Trajectory records from two SVP drifters are used to be compared with the simulations using wavelet-based rotary spectral analysis. Numerical simulations show that the Taiwan Banks divide, separate and guide the southward propagating branch of an M2 tide along the western boundary of the Taiwan Strait as Kelvin Waves. The eastern branch which dissipates its energy by the bottom friction around the Taiwan Banks and the Penghu Islands gradually diminished along the south-eastern boundary of the Taiwan Strait. The northward propagating M2 dominates the tides in the waters east of the point 119.5°E, 22.8°N along the Penghu Channel instead. Moreover, the patterns of the standing waves of the 1/4 diurnal constituents are identified, it could be therefore hypothesized that the amplification of semi-diurnal tides in the Taiwan Strait is primarily due to the condition of two open ended channel co-oscillations. Co-oscillations are demonstrated by numerical simulations as well as an SVP drifter record which indicates that the Taiwan Banks act as the generator of the shallow water tides. The shallow water constituents are amplified and propagate southward to the South China Sea and diminish rapidly. From the SVP drifter data and short-time rotary spectral analysis, it is also interesting to note that the 1/5 diurnal constituent becomes significant compared to other shallow water constituents in the local vicinity of Taiwan Banks.


InTrOducTIOn
The Taiwan Strait (Fig. 1), about 180 km wide, 350 km long and located between Taiwan and mainland China, is part of the East China Sea continental shelf and is generally shallower than 80 meters (Huang and Yu 2003). To the north, the strait is connected to the East China Sea, which is essentially a broad continental shelf extending hundreds of kilometers offshore and conversely the south is open to the South China Sea, where the shelf gradually narrows. A large submarine valley, the Penghu Channel, is located southwest of Taiwan and extends northward until it reaches the Changyun Rise, it is delimited by 60 m isobaths northeast of the Penghu Islands, past which it turns west and passes north of the Penghu Islands and continues northeastward (Fig. 1). Northwest of the Penghu Channel is the Wuchui Trough, which is connected to the Penghu Channel via the Southwest Depression, which runs approximately NW-SE. North of the Penghu Channel is the Changyun Rise, with depths between 30 and 50 m. Southwest of the Penghu Terr. Atmos. Ocean. Sci., Vol. 21, No. 1, 45-69, February 2010 Islands is the Taiwan Banks, formerly known also as the Formosa Banks (Fig. 1).
The Taiwan Banks is a large NE-SW oriented bathymetric feature near the southern end (23°N, 118 -119°E) of the Taiwan Strait. It is generally defined by 40 m isobaths with a large portion shallower than 20 m (Lin et al. 1998) and it covers an area in excess of 10000 km 2 . There are more than 30 shoals shallower than 15 m in this area. Using satellite SAR images, Huang et al. (2006) showed that over 150 wave-like sand dunes, typically several kilometers long, are distributed over the Taiwan Banks. These regular bed forms originate from the interaction of strong tidal currents with the sea-bed (Hulscher 1996), which is made of coarse and fine sand and conglomerates (Lin et al. 1998). High resolution digital sea floor topography with resolution as high as 500 m, i.e., the TaiDBM v6 from NCOR (Taiwan) and the ETOPO2 from NGDC of NOAA, show areas shallower than 10 m across the Taiwan Banks. Fishermen reported sightings of emerged sand dunes on the Taiwan Banks during ebb tides but these reports are not yet confirmed by insitu or remote sensing observations. The Taiwan Strait is situated between major populated cities of Taiwan and mainland China and therefore the prediction of the tides in this area is of considerable importance for the navigation and the dispersion and transport of pollutants.
Since the early 1940's, British Navy charts and Japanese waterway reports have shown that the tides in the Taiwan Strait are primarily semidiurnal, with a maximum amplitude of 2.5 m in the middle part of the strait (Defant 1961). It is known today that the tidal range in excess of 6 m is relatively small, about 0.5 and 1.2 m at the north and the south ends respectively (Liu 1996). Simulations with a primitive three-dimensional numerical model suggest that the tides in the Taiwan Strait result from the interaction of two incident tidal waves coming from both ends of the strait (Zheng et al. 1982;Chen 1983). Based on speculations of shoaling and diffraction effects, Hwung et al. (1986) proposed that the southward propagating tide is dominant but they could not explain the observed tidal amplitudes in the central part of the Taiwan Strait. Lin et al. (2000Lin et al. ( , 2001 compared the observed tides along the western and eastern coasts of the Taiwan Strait with the results of a mild-slope, finite-element model without background rotation. In their study, the width of the Taiwan Strait together with the East China continental shelf was demonstrated to be about half the wavelength of the semidiurnal tides and the standing tides were found to be oscillating and causing the anomalous amplification of the M 2 tide in the central part of the Taiwan Strait. More recently, two and half years of current observations in the Taiwan Strait from moorings and bottom mounted ADCP's (TS-NOW project, see Jan et al. 2001 were compared with simulations from two dimensional numerical models. The results show that as the Pacific M 2 tide propagates into the East China Sea, its amplitude increases from 0.5 to 1.5 m due to shoaling topography and quarter-wave resonance (Jan et al. 2004). Part of this amplified tidal wave propagates southward along the China coast and enters the Taiwan Strait, where further amplification occurs. The southward propagating wave interacts with the part reflected from the abruptly deepened topography south of the Taiwan Strait leading to the coexistence of a progressive wave to the west and a partial standing wave in the eastern part of the strait. These investigations, however, neglected the effect of the Taiwan Banks, which are likely to be important.
The tidal waves that propagate from the Pacific Ocean into the Taiwan Strait become distorted and asymmetric due to bottom friction effects. The shallow water constituents are generally small compared to the semidiurnal and diurnal constituents but they are important to understand the details of the tides in shallow water (Pugh 1987).
The aim of this paper is to investigate the role of the shallow Taiwan Banks on the tidal dynamics with a numerical model. The shallow water constituents, which are shown to be important, largely depend on the bathymetry; the shape of coastline as well as the local tidal regime, and will also be discussed.

The ShAllOW WATer cOnSTITuenTS
Tidal waves in shallow water become distorted and asymmetric because of bottom friction effects and flow curvatures (Defant 1961;Pugh 1987). Several harmonics, with frequencies higher than the original astronomical tidal constituents, can arise. Harmonic analysis suggests that nonlinear distortion effects create compound and overtides normally referred to as shallow water tides, within the diurnal, semidiurnal, quarter-diurnal and even in higher constituent bands.
The bottom friction effect is one of the major causes of tidal energy dissipation in shallow water. Amin (1993) pointed out that bottom friction plays a major role in setting up higher harmonic constituents in shallow water as it redistributes tidal wave energy from lower frequency constituents to higher harmonic bands. The drag of the current on the seabed and the associated energy transfer can be modeled by the kinematics of a frictional boundary layers theory (Amin 1993). The stress due to bottom friction (τ b ) is: where C D is the drag coefficient, t is the sea water density and u is the current velocity. It follows that the work (W) done by the bottom stress per unit time can be expressed by: and W is assumed to balance the energy dissipation that occurs in the boundary layer. As the tidal current flows across a shoal of very small depth, the speed increases due to the conservation of mass. The energy dissipation will therefore grow as the third power of the current speed, Eq.
(2). The frictional loss of momentum per unit volume of fluid increases with decreasing water depths. Due to the quadratic nature of the frictional momentum loss, the increased loss during the ebb phase is greater than the decreased loss during the flood phase. Also, for a progressive long wave in shallow water, where the ratio of the amplitude to the water depth is not negligible, the crest of the wave (flood tide) travels faster than the trough (ebb tide) (Parker 1991). This results in a distortion of the sinusoidal wave profile. If the undisturbed sinusoid (e.g., the semi-diurnal constituent M 2 ) is subtracted from the distorted profile, energy is found in the second harmonic (M 4 ) and other even harmonics. Quadratic friction causes maximum attenuation and minimum wave propagation velocity at both maximum flood and maximum ebb with the opposite occurring at slack waters. The result of this symmetric effect is a third harmonic, M 6 (Parker 1991). When the above mentioned second order shallow water tidal constituents, η 2 , interact with the first order constituents, η 1 , third order tidal waves arise. The third order shallow water constituents, which are S 6 , M 6 , 2MS 6 , 2MS 2 , 2SM 6 , 2SM 2 , are computed from 1 Pugh 1987). Similarly, higher harmonic shallow water constituents can be computed. With this approach, 10 shallow water constituents are obtained by considering only 2 individual source constituents. There are 30 shallow water constituents that originate from the four major constituents O 1 , K 1 , M 2 , S 2 (Pugh 1987).
Shallow water tides can also be due to nonlinear momentum advection terms, through which the effects of flow curvature is significant. It is instructive here to show the effects of nonlinear momentum advection terms, for a tidal wave propagation in a one-dimensional channel with no background rotation, bottom friction and diffusion terms are discussed below [see also Parker (1991) and Bowden 1983)]. The boundary conditions are zero vertical velocity at the bottom and at the free surface. The continuity equation is: where h is the mean water depth, η is the tidal water elevation respected to the mean water depth and u is the current velocity. By applying the hydrostatic approximation, the momentum equation is obtained as: where g is the gravity acceleration. The shallow water tides originate from the terms with flow curvature, i.e., the u · ∂u / ∂x in the momentum equation and η ∂u / ∂x in the continuity equation. To the first order approximation (indicated by suffix '1'), the solution of Eq. (4) is: The ω is the frequency and the k is the wavenumber. Substituting Eq. (5) into Eqs. (3) and (4), the second order approximation solution (indicated by suffix '2') is: The second term on the right hand side is the harmonic of the original wave. Its amplitude increases linearly along the x axis and is proportional to the amplitude squared of the original tidal wave. In sum, the nonlinear terms alone induce overtides even when bottom friction is not considered. Furthermore, if the water elevation η at x = 0 is the superposition of two source constituents with different amplitudes and phases, the solution of η is: . Eq. (8) shows that the nonlinear interaction of two constituents induces overtides as well as components in higher harmonic frequency bands (Bowden 1983). From the above discussion it can be expected that the most significant shallow water constituents are generated in the regions where bottom friction is important and in the vicinity of shoals or headlands, where the flow path is curved. Fanjul (1997) used a three dimensional numerical model to investigate the dynamics of shallow water tides in the North Atlantic Ocean (20 ~ 48°N, 34°W ~ 0°). He showed that the effect of the nonlinear momentum advection terms become non-negligible for the prediction of shallow water tidal elevation. The energy transferred from lower frequency to higher frequency was also identified. The authors took M 2 as an example to demonstrate the effects of the nonlinear term ∂ηu / ∂x in the continuity equation, of the convective term u · ∂u / ∂x in the momentum equation as well as of the friction terms. Fanjul (1997) estimates that 40% of the energy of M 4 is gained from M 2 through the nonlinear effects, 10% from the convective effects, and 50% is gained due to the energy transfer by bottom friction.

model description
The Princeton Ocean Model, or POM (Mellor 2004) was used in this study. POM has been successfully used to model the tidal circulation in domains of widely different scales, such as ocean-basins, marginal seas, continental shelves, straits and harbors (http://www.aos.princeton. edu/WWWPUBLIC/htdocs.pom/Forecasting.htm) thereby demonstrating that POM is an appropriate choice for our study. POM is a full three dimensional, sigma coordinate, finite-difference primitive equation ocean circulation model. It solves three-dimensional primitive equations on an Arakawa C grid. It is based on the assumption of hydrostatic pressure and the Boussinesq approximation is also used. A Laplacian operator handles horizontal mixing, and the coefficients of vertical mixing are determined by a Mellor-Yamada turbulence closure scheme (Mellor and Yamada 1982). The computational scheme of splitting mode is applied in POM to improve the computational efficiency. The scheme integrates the barotropic (external mode) and baroclinic (internal mode) equations at their respective time steps and thus the computation time could be reduced. Since we do not intend to include a discussion of the effects of the wind driven current and the stratification, the atmospheric input was not considered. The salinity and temperature are regarded as constants and the density is assumed to be uniformly distributed along with the water depth in present study. The model equations and the numerical algorithms used to solve them have been described in detail by Mellor (2004). The function for wet/dry mopping was not used as coastal inundation was not a focus in this study. The focus here is on the results of barotropic tidal modeling.

The model domain
We used two computational domains: the first is a large scale domain that includes the Asian marginal seas and the western Pacific and the second is a smaller scale domain that covers the Taiwan Strait region.
The large scale domain (99 ~ 130°E, 5 ~ 40°N, Fig. 2) extends much beyond the Taiwan Strait region in order to reduce the effect of erroneous boundary conditions. With this set-up, the eastern boundary, which is located in the deep Pacific, is the main open boundary with external forcing, while the other three boundaries are sheltered by lands. The eastern boundary is forced by the dominant deep ocean constituents, since, given the large water depth in the Pacific (> 5000 m), the local shallow water tide generating potential can be neglected. The grid has 186 × 210 nodes, with a resolution of 10 min.
The small scale domain (116 ~ 122°E and 22 ~ 26°N) has 250 × 250 nodes with variable spatial resolutions wherein the grid generated used a SeaGrid Orthogonal Grid Maker (Denham 2006, http://woodshole.er.usgs.gov/staffpages/ cdenham/public_html/seagrid/seagrid.html). The SeaGrid allows condensing the grid points in the area of interests, thus helping to avoid an excessive computational load. Figure 3 shows the details of the grid in the Taiwan Strait region. The smaller grid spacing was set at the 119.0°E, 23.0°N (about 5.2 and 3.6 km in the zonal and meridional directions respectively). The grid spacing is then increased linearly to the boundaries. Twelve vertical layers are used.

The Bathymetry
The ETOPO2 (Earth Topography Two Minute) bathymetry data from the US National Geophysical Data Center (NGDC) (http://www.ngdc. noaa.gov/mgg/) and the 500 m high resolution TaiDBM (Taiwan Digital Bathymetry Model, Version 6) bathymetry data of National Center of Ocean Research, Taiwan (NCOR) (http://www.ncor. ntu.edu.tw/ODBS/MGG/) were used. The TaiBDM (117 ~ 125°E, 18 ~ 27°N) was compiled from echo sounding datasets from research vessels and has a better resolution than ETOPO2 in the shallow water regions of the Taiwan Strait. ETOPO2 was used to cover the areas for which the TaiBDM v.6 has no data.

The Open Boundary
The tidal forcing was applied to the open eastern boundary of the large scale domain by specifying the time dependent sea level generated by the eight dominant diurnal and semi-diurnal tidal constituents (M 2 , O 1 , K 1 , S 2 , P 1 , N 2 , K 2 , and Q 1 ), as calculated by the NAOTIDE global tidal model (Matsunoto and Takannezawa 2000). The surface elevation η(t), due to n constituents, is given by where a i is the amplitude, σ i is the frequency and ε i is the phase. NAOTIDE was developed by assimilating five years of TOPEX/POSEIDON altimeter data into a hydro-dynamical model. A comparison of three global tidal models, CSR4.0, NAOTIDE and FES, with observations in the East China Sea and the seas around Japan shows that NAOTIDE performs better when estimating M 2 , K 1 , S 2 , N 2 , K 2 , and Q 1 (Matsumoto 2000). For this study NAOTIDE prescribes the sea-level at the open boundary by adding, at each time step, the incremental adjustment due to the constituents. The Taiwan Strait model started from an initial condition interpolated from the results of the Asian marginal Seas and the western Pacific model. Along the boundaries of the Taiwan Strait model, the elevation and current velocity were obtained by spatial interpolation of the corresponding coarser model. In the present study a one-way nesting method is adopted. The boundary conditions of the nested region are provided by the results of the coarser model.

model Parameters
A splitting mode technique is used in POM. The external mode, which is barotropic, calculates the fast two dimen- sional dispersive surface waves, and an explicit numerical scheme is used to reduce the computation time. The internal mode, which is baroclinic, calculates in full the three dimensional gravity waves. The implicit scheme is adopted for the internal mode to gain higher spatial resolution and avoid computational divergence. The Asian marginal seas and the western Pacific model was run with a 60 s time step in the internal mode time-step, and 2 s time step in the external mode time-step. The horizontal eddy viscosity formulation follows the Smagorinsky Diffusivity (Mellor 2004) with the HORCON parameter (Mellor 2004) set to 0.20. In the POM model, the stress induced by the bottom friction is described by where C z is the non-dimensional coefficient, which is a function of the roughness length of the bed, Z 0 , and the water depth: where, ĸ = 0.4 is the von Karman constant, ζ is the bottom topography and σ kb -1 is the depth of the layer overlying the bottom one in sigma coordinates. In other words a logarithmic current profile is assumed and the bottom friction is determined by the current speed of the layer nearest to the bed. Numerical investigations (Kuo et al. 1996) showed that in the oscillatory boundary layer flow, the vertical velocity distribution near the bottom may deviate from the logarithmic profile. The small deviation of the current profile to the logarithmic law will result in significant errors when the profiles are used to calculate the bed roughness length and shear stress. Kuo et al. (1996) present the deviation amplitude of calculated roughness length and bottom shear stress may exceed 100% when strong flow acceleration occurs. This will further incur the physically unrealistic roughness height that varies with tidal phase. It should be noted that with the practical range of vertical grid spacing, the error in bottom stress calculated by the logarithmic profile becomes non-negligible, particularly when shallow water overtides are present. For numerical simulations of flows in shallow water region, such as estuaries and channels, the Manning formula for bottom friction is often adopted (Liu at al 2003) and C z of Eq. (11) becomes: where n is the Manning coefficient, h is the water depth and g is the acceleration of gravity. As stated by Eq. (11) through Eq. (13), bottom friction is determined from the current speed averaged over the whole water depth together with the coefficient n at the corresponding location. The velocity is therefore much less sensitive to the logarithmic profile assumption.
In numerical simulations of tidal waves propagating over Taiwan Banks, the abrupt bathymetry shoals, an undesirable spurious circulation was created when the sigma coordinate system was applied and resulted in a velocity error field. Since the sigma-layers were squeezed together in very shallow water (less than 20 m), the errors in the logarithmic current profile could be easily augmented and thus the bed stress might not properly reflected. This pressure gradient error of sigma coordinate over steep topography could be made smaller by increasing the grid resolution (Mellor 1994(Mellor , 1998 or using a horizontally averaged formula as in Eqs. (12) and (13). We found that more stable computations were obtained when the Manning formulation was adopted with the coefficient n set to 0.032.

Validation of the model
To validate the model, a year's worth (2005) of sea level data from 34 tidal stations along the coasts of Taiwan and mainland China were used for comparison (Fig. 3). It should be noted that the data from Lialuo Wan station is considered to be invalid due to the fact that the phase reveals significant discrepancies with regard to the overall trend along the coasts and the nearby (25 km) Xiamen Gang Station.
The comparisons of the amplitudes and phases of M 2 , S 2 , K 1 , and O 1 obtained from the harmonic analysis of the observations with their numerically simulated counterparts are shown in Tables 1 through 4, respectively [see Pawlowicz (2002) for a description of the methodology used here]. The error indexes used to quantify the comparison are the amplitude error (APE) and the Phase lags (∆g). The root mean square error (RMSE) of the amplitudes is used to express the overall performance, where A is the amplitude, φ is the phase and subscript ob and si denote the values from observation and numerical simulation, respectively. As illustrated in Figs. 4, 5 and Tables 1, 2, the simulated amplitudes of the semi-diurnal constituents agree with the observations quite well. Overall, the simulated M 2 amplitudes are about 0.041 m larger than the observations and 0.001 m smaller for S 2 . The RMSE of the M 2 and S 2 amplitude simulations are 0.0848 and 0.0415 m, respectively. The maximum APE of the model occurs at the northeast boundary of the Taiwan Strait. The maximum phase errors of semi-diurnal constituents occur along the southwest coast of Taiwan (Kao-hsiung Kang and Tung-Kang). The averaged phase error is about 5 degrees. The overall RMSE of amplitude differences for K 1 is 0.0251 m. The model amplitude is under-estimated by a few centimeters on the northern coast of Taiwan and the same magnitude over-estimation along the Taiwan coast from the middle part of the island to the south (Figs. 6 and 7). As for O 1 , the amplitude RMSE is less than 0.0238 m. The phase discrepancies for K 1 and O 1 are always less than 27° (Tables 3 and 4). The comparison between modeled and observed data for M 2 , S 2 , K 1 , O 1 shows that the RMSE error indexes are smaller than the ones obtained from recent simulations in the region (Lin et al. 2000;.

SImulATIOn reSulTS And dIScuSSIOn
The tidal elevation was simulated from 2005/01/01 to 2005/12/31 and the spatial characteristics of the dominating M 2 and other shallow water constituents induced by nonlinear effects are described and discussed in this section. In order to verify and validate the model results, data from two surface velocity program (SVP) drifters (Niiler 2001) were used for comparison. On 7 May 2005, 25 SVP drifters were towed at a 15 m depth in the South China Sea (SCS) west of the Luzon Strait to measure the surface expression of the internal tide and non-linear internal waves that reach their maximum amplitude during spring tides (Ramp et al. 2004). See Centurioni et al. (2007) for a description of this experiment. Two drifters marked 56440 and 56420 exited the SCS through the Taiwan Strait after several months of deployment and the trajectories and velocities of the devices that drifted in the vicinity of the Taiwan Banks are used here. Their trajectories are illustrated in Figs. 8 and 9.
To compute the spectrum of the tidal current at each point of the drifter tracks, the wavelet based short-term rotary spectral analysis was applied (Chien 2005). This methodology allows estimation of the frequency-temporal structure of the drifter records at each time/location. The wavelet based rotary spectrum RWT is, The RWT is a density function of a and b, which are the parameters representing the frequency and time of the signal, respectively. w(t) is the velocity function of drifter, w(t) = u(t) + iv(t), and the ψ is the wavelet function. We adopted the Morlet wavelet function in present study, In the calculation, the Morlet wavelet functions were first re-scaled and shifted by alternating the variable x by given values of a and b to form the wavelet function sets. The integral of the inner product of the signals of drifter velocity and the wavelet function sets was then taken. The major advantage of using the wavelet-based rotary spectrum rather than Fourier transform based rotary spectrum is that the wavelet-based spectral analysis is capable of dealing with nonstationary signals. It gives the spectral structure of a signal at a specific time and location with reasonable frequency  and temporal resolutions. The Fourier transform gives the information of mean spectrum over the whole duration of the signal. This fact is crucial as the drifters' velocity signal in the present study varied with time/location. The results of the wavelet-based rotary spectra of the signals retrieved from Argos drifters 56420 and 56440, when they were approaching the Taiwan Banks in July and September of 2005, are illustrated in Figs. 10 -13. Figures 10 and 12 are the current velocity oscillation of each constituent, and Figs. 11 and 13 are the temporal variation of velocity amplitudes of SVP 56420 and 56440, respectively.

Tidal Propagation
The amplitude of M 2 (Fig. 4), which is the dominating constituent in the Taiwan Strait, is in excess of 1.65 and 2.05 m in the middle of the Taiwanese coast (Tai-Chung station) and of China's coast (Wuqiu Yu station), respectively. The phase lines also show that the M 2 wave propagates from the north into the Taiwan Strait. This wave is essentially a progressive barotropic Kelvin wave that moves southward along the western boundary of the Taiwan Strait. These simulated results are consistent with the observations, for which the phase increases to the south and along the western boundary of the Taiwan Strait (Lin 2001). The pattern of the phase lines of S 2 (Fig. 5) is similar to the M 2 , the amplitude is less than 1/3 of M 2 . The co-tidal contours of M 2 in the vi-cinity of Taiwan Banks and the co-phase lines are relatively dense in the region of Taiwan Banks and southwest of the Penghu Channel since the celerity of the M 2 wave is smaller in the shallow water. The mean rate of energy dissipation per unit area due to bottom friction could be estimated using the following equation:  To estimate the tidal energy dissipation, the integration in Eq. (19) was taken over a one year run using the 30 minute depth-averaged current data. The rate of energy loss due to bottom friction could be yielded by dividing the total current energy per unit time per unit area by the value obtained from Eq. (19). The result is shown in Fig. 14. It could be found that the maximum of the energy dissipation is about 90% over the Taiwan Banks area.
On the western part of the lower domain, the southward M 2 tide west of the Taiwan Banks propagates along the coast of mainland China as the Kelvin Waves, and part of it then diffracts around the point 118°E -22.5°N, from the southwest to the south, due to the topographic slope in the south of the Taiwan Banks (Fig. 4). In the middle part of the lower domain, where the Taiwan Banks are located, more than 90% of the southward M 2 tide is dissipated by bottom friction (Fig. 14) and the amplitude decreases from 0.4 to 0.05 m while the celerity also decreases rapidly when the wave reaches the shallowest area.
Finally on the eastern side of the lower domain, the southward M 2 tide becomes weak after it propagates around the Penghu Islands where it loses most of its energy (Fig. 14). On the other hand, as marked with the alphabets from a to e in Figs. 4 and 5, five tidal gauge stations, i.e., the Chiangchun Kang, Kuo-sheng Kang, An-ping Kang, Kao-hsiung Kang and Tung Kang stations, aligned along the south-west coast of Taiwan, which is on the eastern side of the Penghu Channel entrance. "Kang" is the Mandarin translation of harbor. These tidal gauges are located at the harbor mouths, which could be regarded as in open waters. For other gauges with similar names where "Kang" is omitted, e.g., Chiang-chun, An-Ping and Kao-Hsiung, they were setup in the inner harbor or estuaries which are closer to urban areas. As the tidal signals measured from those stations feature phase lag to those obtained from the harbor entrances, we only focus on the tidal signal obtained from open waters here. According to the measured data obtained from the five stations, the decreasing phase along the coastline section from the north to the south can be identified. This indicates that the M 2 tide propagates from the south to the north as another Kelvin Wave along the Penghu Channel.
The simulated phases of northward propagating and southward propagating M 2 tide are gridded and the corresponding co-phase contour lines are drawn. From Figs. 4 and 5, it seems that 119.5°E, 22.8°N, which is at the western bank of the Penghu Channel is the point that separates the northward and southward tidal waves.
In the above scenario, the northeast-southwest Taiwan Banks act as the divide that separate the southward M 2 tide and guide the directions of the western branch as the Kelvin waves along the coast of Mainland China. The eastern branch M 2 wave then diminished due to the bottom friction. Along the coastline of Taiwan, the southward tidal waves are weakened by the Penghu Islands, while the northward M 2 wave dominates the tides along the Penghu Channel. From the numerical simulations, it could be seen that the 119.5°E, 22.8°N is the critical point that separates these two waves. At this point, the rotation-like current that was incurred by the tidal residual could be formed. The track of a SVP drifter marked 56440 in Fig. 9 shows a counter-clockwise rotation around the above-mentioned point, beginning in 7 September 2005.
The drifter's tracks of 56420 and 56440 were denoted by solid lines and the locations of the drifters at 0000 UTC each day were marked by triangles respectively in Figs. 8 and 9. The drifters, which were equipped with drogues, were designed to move in the same velocity as water particles of the upper ocean mix layer. The orbital tracks of the tidal oscillations were not closed due to the existence of the mass transport of the upper ocean layer. This surface current could be regarded as the superposition contributed from two different mechanisms, i.e., the summer monsoon associated Ekman current in Taiwan Strait (Centurioni et al. 2004(Centurioni et al. , 2007 and the shallow water, Stoke's drifting current. The summer monsoon currents a larger, more pronounced circulation scale and have much lower oscillation frequencies than astronomy tides. Thus, they feature comparatively Fig. 14. The spatial distribution of the average energy dissipation (W m -2 ) due to bottom friction.
steadier than the tidal current as well as the Stoke's drifting current induced by shallow water tidal residual. The mean speed of northeastward surface current measured by drifter 56440 was about 0.25 m s -1 , which was only half of the magnitude as the mean speed retrieved from drifter 56420. In present study, the shallow water tidal oscillations with diurnal or higher frequencies are focused, thus the components below 0.7 cycles per day are omitted as shown in Figs. 10 -13.
As the drifter 56440 moved with a comparatively slower monsoon surface current, it was captured by a tidal residual current at 119.5°E, 22.8°N. The current direction is paralleled to the tidal wave propagation direction, which is perpendicular to its co-phase lines, and as a result drifter 56440 rotates counter-clockwise around the point. It took 120 hours to drift around the point (Fig. 9). The behavior of the drifter is consistent with the simulation by present numerical model. Hwung et al. (1986) suggested that the northward tidal waves that enter the Taiwan Strait come from the Luzon Strait and interact with the southward wave. Recent numerical experiments (Jan et. al. , 2004 suggest that the northward M 2 is generated by the reflection of the incoming wave due to abrupt step-like topography. Jan et. al. (2004) stated that the step-like bathymetry in the Penghu Channel causes the reflection of the M 2 tide, and the energy flux of M 2 tide across this boundary from Luzon Strait can be neglected. In our study, we showed that the simulation agrees with the observational evidence from SVP drifter. The area where the northward waves predominate is demonstrated. Still the source of the northward M 2 tide needs further investigation.

Shallow-Water harmonic Generation
From the corresponding wavelet-based rotary spectra of drifter 56420 and 56440 (Figs. 10 -13), periodical oscillations of drifters' velocity, which were dominated by semidiurnal tides, can be obviously identified. As the drifter 56420 and 56440 approached the Taiwan Banks on 2 September and 12 July 2005 respectively, the strength of high frequency velocity oscillations was intensified as shown in Figs. 10 and 12. As the drifters reached the shortest distance from the Taiwan Banks, the energy of 1/5 diurnal tides increases significantly (Figs. 11 and 13). These higher frequency tides decay rapidly as the drifter moved away from the Taiwan Banks. The difference between the characteristics of the behavior of drifter 56420 and 56440 is that the 56440 trajectory featured higher wave nonlinearity as there were more significant components in the higher frequency bands as shown in Figs. 11 and 13. This demonstrates the effects of nonlinearity of local shallow water tidal waves in the period when the drifter 56440 passed along the southern border of Taiwan Banks.
Higher-order tidal harmonics are generated in shallow water regions where bottom friction and non-linear interactions are important. A harmonic analysis (Pawlowicz 2002) of the one year simulated water elevation data was performed (Figs. 15 to 18). MO 3 , M 4 , 2MK 5 , and M 6 are used to represent 1/3, 1/4, 1/5, and 1/6 diurnal shallow water tidal constituents, respectively. The largest of these higher harmonic tidal amplitudes are generated in the shallow regions in the immediate vicinity of Taiwan Banks. It could be found that, in the area of Taiwan Banks, the modeled 1/4 diurnal tide has the largest amplitude (Fig. 16) and the 1/5 diurnal tide (Fig. 17) has a smaller amplitude than the 4 shallow water constituents. It should be noted however, the observed velocity amplitude of the 1/5 and 1/6 diurnal constituents were even larger than the 1/3 and 1/4 diurnal tides for both drifters (Figs. 11 and 13).
From the results of harmonic analysis (Pawlowicz 2002), the amplitude of the shallow water harmonics range from 2% to 10% of the amplitude of M 2 in the Taiwan Banks region and decrease toward the south and the north ends of the Taiwan Strait. The simulated phase of MO 3 shows that the waves generated by the Taiwan Banks propagate both southeastward and northward to the South China Sea and East China Sea, as M 6 were generated by the bathymetric change at the East China Sea located at the northeast region of the Taiwan Strait and then being enhanced in the Taiwan Banks region. 2MK 5 is generated near Pingtan, Fujian on the western boundary of the Taiwan Strait and then propagates both to the north and the south. These phenomena are compared and confirmed with the observed harmonic constants along the coast of Taiwan (Lin et al. 2001). The phase pattern of M 4 was more complicated. There exist nodal points in the Taiwan Strait.
The amplitude of M 4 is larger than 2MK 5 and M 6 near the Taiwan Banks and decays rapidly with increasing water depth. It could be concluded from the above discussion of amplitude and co-phase patterns that the Taiwan Banks is not the source of the generation of shallow water constituents in the Taiwan Strait, but it significantly intensified their strength. Our wavelet-based rotary spectral analysis of drifters records also demonstrates that there exist background higher harmonic constituents in the Taiwan Strait, and these components are magnified in the vicinity of the Taiwan Banks.

Tidal Wave resonance
To investigate the amplification of semidiurnal constituents in the mid section of the Taiwan Strait, Lin et al. (2000Lin et al. ( , 2001 analyzed the harmonic constants for more than 30 tidal stations along the western coast of Taiwan. It was found that the maximum amplification factor was 21.19 (10.72 for S 2 , 14.60 for M 2 , 17.49 for N 2 , 20.53 for O 2 , and 21.19 for 2NS 2 ), which is at least twice the magni-    tude that could be generated by the interaction of solitons in the shallow water (Soomere 2004). This result implies that the resonance effect would be the only possible cause of this amplification.
The resonance is caused by the existence of standing waves in the channel. Standing waves have zero celerity and most of the wave energy is therefore conserved. With increasing energy input from incident waves the amplitude of the standing waves is amplified. The maximum amplitude occurs when energy input is balanced by energy dissipation from bottom friction, viscosity and energy leakage. For a channel with only one open end the co-oscillation standing wave occurs when the following quarter-wave resonance condition is met (Mellor 1998 where L is the length of the continental shelf or of the basin, ω and λ are the wave angular frequency and wave length, respectively and c gh = is the phase speed. The frequencies of the normal modes for this quarter-wave resonance can be written as: Since the averaged phase speed of tidal waves in the Taiwan Strait is c = 23 m s -1 and the frequency of the semidiurnal tide M 2 is ω = 1.41 × 10 -4 rad s -1 , we now examine the boundary condition at the ends of the Taiwan Strait based on the above theories. Since the actual length of the Taiwan Strait is obscure, we could either measure the length of it from 300 km (north entrance: 120.5°E, 25.5°N ~ south entrance: 119°E, 23.25°N) to 500 km long (north entrance 121°E, 26°N ~ south entrance: 118.25°E, 22.25°N).
For the case where both ends of the channel being open, according to Eq. (26), for the first mode (n = 1) and the second mode (n = 2), the resonance of M 2 and M 4 tides will occur if the length of channel L = 500 km is presumed. For the case of only one end open, according to the quarterwave resonance condition, the resonance of M 2 and 2MK 5 will occur when Eq. (24) is satisfied for n = 1 and n = 3 as we take L = 300 km. It is noted that the corresponding wave period of n = 3 is about 4.8 hr, which the periods of the 1/5 diurnal constituent.
As the resonance of M 2 could occur in both cases just by presuming different lengths of the Taiwan Strait, we now focus on the patterns of higher harmonic constituents. The model results show that two nodal points in the M 4 simulation could be identified as denoted by arrows in Fig. 16, while there are no nodal points found in the 2MK 5 and M 6 co-tidal charts (Figs. 17 and 18). Since the nodal points are the signatures of standing waves we can hypothesize that the tidal co-oscillation in the Taiwan Strait is primarily due to the condition of two ends open channel. From the previous discussion, the whole south end of the Taiwan Strait should be considered as open for transmitting tidal waves.

cOncludInG remArkS
The tidal dynamics of the Taiwan Strait was investigated with the POM model. A good agreement with the available field data of 34 water elevation stations was shown. The error indexes of the numerical simulation were smaller than those of previous studies in this region. The RMSE of M 2 , S 2 , O 1 , and K 1 are 0.0848, 0.0415, 0.0251, and 0.0238 m, respectively.
In order to obtain a better understanding of the effects of the Taiwan Banks on the tidal waves, high resolution numerical simulations were carried out. The effects of the Taiwan Banks on the tides mechanisms can be summarized as follows: first, it divides, separates and guides the southward propagating branch of M 2 tide along the western boundary of the Taiwan Strait as Kelvin Waves. The eastern branch, which is dissipated the energy by the bottom friction around the Taiwan Banks and the Penghu Islands gradually diminished along the south-eastern boundary of the Taiwan Strait. The northward propagating M 2 dominates the tides in the waters east of the point 119.5°E, 22.8°N along the Penghu Channel instead. It contributes to the resonance mechanism and induces anomaly amplification of semidiurnal tide in the mid-section of Taiwan Strait. Moreover, as the significant patterns of the standing waves of the 1/4 diurnal constituents are obvious, it could be hypothesized that the amplification of semi-diurnal tides in the Taiwan Strait are primarily due to the condition of open channel co-oscillations at the two ends.
Second, the Taiwan Banks enhance the strength of shallow water constituents in the Taiwan Strait. The shallow water constituents are amplified, propagate southward to the South China Sea and diminished rapidly. From the SVP drifter data short-time rotary spectral analysis, it is also interesting to note that the 1/5 diurnal constituent becomes significant compared to other shallow water constituents in the local vicinity of the Taiwan Banks. This phenomenon was not well resolved in the numerical models.