滕 飛 方國洪, 魏澤勛 徐曉慶 崔欣梅 吳 頔

(1. 中國海洋大學(xué) 青島 266100; 2. 國家海洋局第一海洋研究所 青島 266061)

為了研究底摩擦對潮汐的影響和提高潮汐數(shù)值模擬的精確度, 我們分別應(yīng)用 Chezy型和廣義Manning型摩擦關(guān)系模擬渤、黃、東海的潮汐。

Chezy于 1775年首先提出定常流動水渠中摩擦力與流速的關(guān)系, 后來被稱為 Chezy公式(Dronkers,1964)。在潮汐動力學(xué)中可寫成

其中Fb為摩擦力, ρ為流體密度, g為重力加速度, C為 Chezy系數(shù), u為流速。如果 C為常數(shù), 式(1)表明摩擦力與流速的平方成比例, 為方便常常用系數(shù)r代替Chezy系數(shù), 取

r(亦即C)取常數(shù)的情形常被稱為平方摩擦, 也可稱為 Chezy 型摩擦(Lefèvre et al, 2000)。Taylor首次研究了愛爾蘭海的潮汐摩擦問題(Taylor, 1920), 證明采用r=0.0016至0.002所得到的潮能消耗與依據(jù)潮汐和潮流觀測值所算出的潮能消耗基本一致。Proudman(1953)按照Taylor的方法進(jìn)一步給出r的具體數(shù)值為0.0026。此后, 大多數(shù)潮汐數(shù)值模擬中均取r=0.002至0.003的數(shù)值。

Manning于 1890年通過實驗, 提出 Chezy系數(shù)并不是常數(shù),而是與水深 h有關(guān)(Dronkers, 1964), 他給出的關(guān)系式,

稱為Manning公式, 其中n稱為Manning粗糙度系數(shù),將式(3)代入式(2), 可得

式(4)表明, 水深越淺, r值越大。

對于中國近海, Fang等(1987)用三種不同的方法計算了杭州灣潮流的摩擦系數(shù),得到的 r值介于0.0005至0.00067之間, 顯著小于Proudman(1953)和后來的數(shù)模研究所采用的數(shù)值。值得注意的是,杭州灣水深很淺, 大約10m左右。周朦等(1987)估計了渤海的摩擦系數(shù), 得到r=0.0013, 也比 Proudman(1953)的數(shù)值小很多, 而渤海的水深大約在20m左右。這些結(jié)果表明, 與其他模擬研究相比,渤、黃、東海海底摩擦系數(shù)較小, 而且與水深的關(guān)系可能與Manning公式不一致。我們稱平方摩擦系數(shù)為常數(shù)的形式為Chezy型, 稱平方摩擦系數(shù)與水深的冪函數(shù)成比例形式為廣義Manning型。這里“廣義”指冪值可以是任意實數(shù),不局限于Manning公式中的-1/3。本研究嘗試在渤、黃、東海潮汐模擬中分別應(yīng)用 Chezy型以及廣義Manning型摩擦關(guān)系,并通過與實測對比來確定其最優(yōu)的參數(shù)值。

1 數(shù)值模式和實測數(shù)據(jù)

1.1 基本方程

本研究采用有限體積海洋模式 FVCOM, 該模式采用σ垂向坐標(biāo)來模擬不規(guī)則的海底地形, 坐標(biāo)變換公式如下(Chenet al, 2006)

由于計算范圍大, 考慮緯度不同所引起的科氏力變化, 采用球坐標(biāo)連續(xù)方程, 球坐標(biāo)系下的σ坐標(biāo)的連續(xù)方程和動量方程如下

式中,R是地球半徑,γ代表計算點經(jīng)度,Φ代表計算點緯度,D代表未受擾動水深,ζ為相對于未受擾動海面的高度,u、v、w分別代表東向、北向和垂向速度,ρ代表海水密度,f代表科氏參數(shù),g代表地球重力加速度,Km代表垂直渦動黏性系數(shù)。本研究采用它的二維模式。

在大范圍的區(qū)域中, 引潮勢產(chǎn)生的平衡潮是不能忽略的, 因此模擬過程必須包括引潮力。對半日潮,平衡潮為

對全日潮, 平衡潮為

其中,Ai和ωi代表第i個分潮的振幅和頻率,Nsemi和Ndiurnal代表模式中參與計算的半日潮和全日潮的個數(shù),βi為地球彈性所產(chǎn)生的訂正因子, 對各個分潮所采用的數(shù)值見表1。

表1 平衡潮訂正因子βiTab.1 Correction factor for equilibrium tides βi

1.2 模式水深

在潮汐計算中, 水深是重要的因素, 且本文實驗方法中廣義Manning型底摩擦系數(shù)與水深關(guān)系密切。對于陸架海域, 直接從海圖上讀取的水深值一般要優(yōu)于國外現(xiàn)有數(shù)據(jù)集的數(shù)值(例 如吳 頔 等, 2015), 因此為了提高計算準(zhǔn)確度, 本文采用的水深數(shù)據(jù)來自兩部分。第一部分是林美華等(1991)水深數(shù)據(jù)集, 該數(shù)據(jù)集的水深系由海圖水深讀取得到, 分辨率為 5′;第二部分是 Etopo5水深數(shù)據(jù), 對于第一部分?jǐn)?shù)據(jù)中沒有涵蓋的區(qū)域從 Etopo5水深數(shù)據(jù)中選取, 最后得到研究區(qū)域完整的水深數(shù)據(jù), 其分布見圖1。

1.3 計算海區(qū)、邊界條件的選取與計算過程

本文研究海區(qū)為 24°—41°N, 116°—128°E 的渤海、黃海和東海陸架海域。模擬區(qū)域中不包括東海的深海區(qū), 是由于深水區(qū)的潮能耗散主要來自內(nèi)潮效應(yīng), 需要另作研究。本研究采用的計算網(wǎng)格水平分辨率在島嶼以及大陸邊界附近約為 1/20°, 其余部分約為1/10°, 共包含26851個網(wǎng)格節(jié)點, 51486個三角單元(圖2), 垂向分12層。FVCOM模型采用內(nèi)外模分離的方式求解。二維外模數(shù)值格式是基于三角形網(wǎng)格的有限體積法, 將連續(xù)方程和動量方程在三角形單元內(nèi)積分后, 通過改進(jìn)的四階龍格-庫塔方法求解。三維內(nèi)模的動量方程采用簡單的顯式和隱式相結(jié)合的差分格式求解。外模時間步長為6s, 內(nèi)外模時間步長比率為10。

圖1 渤、黃、東海水深分布圖(單位： m)Fig.1 Bathymetry of the Bohai, Yellow, and East China Seas

圖2 計算網(wǎng)格Fig.2 The computation grid

計算的初始條件假設(shè)海洋是靜止的, 海表面的擾動水位初始值全部為0, 所有三角形網(wǎng)格中心點的流速u=v=0。本研究采用正壓模式, 取整個海域內(nèi)的溫度和鹽度均為常數(shù), 溫度為18°C, 鹽度為33。在臺灣海峽、東海陸坡、朝鮮海峽共設(shè)有 3個開邊界,開邊界上每個點的水位給定,由 T_tide的預(yù)報程序給出(Pawlowiczetal, 2002),其中 M2、S2、K1和 O1調(diào)和常數(shù)來自DTU10(Chengetal, 2011)數(shù)據(jù)。

模式運行30d。對后15d時間序列進(jìn)行調(diào)和分析,得到網(wǎng)格點上的水位的調(diào)和常數(shù)、每個三角形中心點的潮流橢圓要素以及速度各分量的調(diào)和常數(shù)。

為了對模式結(jié)果和觀測值進(jìn)行比較并優(yōu)化模式參數(shù), 本文在研究區(qū)域內(nèi)一共選取了43個實測站, 其分布圖見圖3。在這43個站的數(shù)據(jù)中, 26個來自國際水文組織IHO潮汐調(diào)和常數(shù)數(shù)據(jù)集, 11個來自TOPEX/Poseidon(T/P)衛(wèi)星軌道交叉點調(diào)和常數(shù)數(shù)據(jù)(Wanget al, 2012), 2個站調(diào)和常數(shù)來自Jan等(2002), 4個根據(jù)TOGA計劃的實測逐時水位資料經(jīng)調(diào)和分析得到。

圖3 實測站位置分布圖Fig.3 Deployment of observational stations

1.4 評價方法

為了評估模式結(jié)果的準(zhǔn)確性和優(yōu)化模式參數(shù),采用模式結(jié)果和實測調(diào)和常數(shù)的向量均方根(RMS)偏差F作為價格函數(shù), 其計算公式如下

其中Fi為各個分潮的均方根偏差

上面k代表實測站位編號, i代表分潮, H代表振幅, G代表遲角, Com代表模擬結(jié)果, Mea代表實測結(jié)果。本文中對兩種底摩擦類型的結(jié)果評價和參數(shù)優(yōu)化都是以價格函數(shù)F值大小為標(biāo)準(zhǔn), F值越小代表模擬結(jié)果越好, 反之越差。

2 實驗方法

2.1 Chezy型摩擦

采用Chezy型摩擦?xí)r, 底摩擦系數(shù)r為常數(shù)。采用不同的r值可以得到不同的F值, 最優(yōu)的r值對應(yīng)最小的F值。為了找出最優(yōu)的r值, 一種較先進(jìn)的方法是通過伴隨方程進(jìn)行反演, 但是這會涉及難度較大的伴隨模式的研發(fā)。對于當(dāng)前的需要優(yōu)化的參數(shù)很少的情況, 我們可采用多次直接模擬的方法。我們先用范圍較大、分辨率較低(即間隔較大)的不同r值進(jìn)行模擬。首先, 取r=(0.0005, 0.0010, 0.0015, 0.0020,0.0025, 0.0030, 0.0035, 0.0040), 結(jié)果得到其中r=0.0010時 F值最小; 然后將分辨率增加一倍, 取r=(0.00075, 0.0010, 0.00125)再進(jìn)行模擬(其中對 r=0.0010不需要重新模擬, 但參與比較), 結(jié)果得到其中r=0.0010時F值最小; 最后再將分辨率增加一倍,取 r=(0.0008725, 0.0010, 0.001125)進(jìn)行模擬(和前面一樣, 其中對 r=0.0010不需要重新模擬, 但參與比較), 結(jié)果得到其中 r=0.001125時, F最小, 其值為11.04cm。本次實驗F值分布曲線見圖4 (其中黑色圓點代表模擬過程中各個r和對應(yīng)的 F值)。由于最小的幾個F值實際上差別很小, 故不作更高分辨率的進(jìn)一步模擬, 而認(rèn)為0.001125為最優(yōu)的r值。

圖4 采用Chezy型摩擦模擬所得價格函數(shù)F值與摩擦系數(shù)rFig.4 The dependence of cost function F on friction coefficient r in Chezy-type relationship

2.2 廣義Manning型底摩擦

為了避免當(dāng)h→0和h→∞時r→0,將廣義Manning型摩擦系數(shù)r取作如下形式

若 h≤h1, 取 r=r1;

若 h≥h2, 取 r=r2; (14)

為滿足作為h之函數(shù)的r的連續(xù)性, a和m需分別等于

本研究中取1h=10m,2h=200m。數(shù)值模擬結(jié)果與實測值的偏差與r1和r2有關(guān)。因而價格函數(shù)F可以寫成

對應(yīng)著最小F值的(r1, r2)組合便是最優(yōu)的組合。

為了尋找最優(yōu)的(r1, r2)組合, 先取范圍較大、分辨率較低的(r1, r2)組合： (r1=0.0002, 0.0004, …, 0.0014;r2=0.0005, 0.0010, …, 0.0040)這 56 組(r1, r2)值進(jìn)行模擬, 并計算出對應(yīng)的 F值。結(jié)果得到(r1, r2)=(0.0008,0.0015)時 F值最小。然后以該點為中心, 縮小(r1, r2)的取值范圍并提高一倍分辨率繼續(xù)進(jìn)行數(shù)值模擬, 并計算對應(yīng)的F值, 結(jié)果得到(r1, r2)=(0.0009, 0.0015)時F值最小。最后再以該點為中心, 進(jìn)一步縮小(r1,r2)的取值范圍并再提高一倍分辨率繼續(xù)進(jìn)行數(shù)值模擬, 并計算對應(yīng)的F值, 結(jié)果得到(r1,r2)=(0.0009, 0.001375)時F值最小, 等于 10.88cm。由于在這個組合附近F值的變化已很小, 我們認(rèn)為它已經(jīng)是最優(yōu)的參數(shù)組合,不再作更高分辨率組合的模擬。由式(15)可得a=0.001375,m=0.14。與原始的Manning公式不同, 該公式的冪值為負(fù)數(shù), 而本研究得到的冪值為正。這表明在渤、黃、東海陸架區(qū), 水深越淺摩擦系數(shù)r值越小。

全部實驗結(jié)果以F值等值線的形式示于圖5, 其中黑色圓點為(r1,r2)第一次實驗取值點, 正方形點為第二次實驗取值點, 菱形點為第三次實驗取值點。為了更好顯示第二和第三次實驗結(jié)果, 我們將圖5方框區(qū)域予以放大, 示于圖6。

圖5 采用廣義Manning型摩擦模擬所得價格函數(shù)F值與參數(shù)組合(r1, r2)的關(guān)系： 全部實驗結(jié)果Fig.5 The dependence of cost function F on parameter combination (r1, r2) in Manning-type relationship： all results

圖6 采用廣義Manning型摩擦模擬所得價格函數(shù)F值與參數(shù)組合(r1, r2)的關(guān)系： 圖5方框區(qū)域放大圖Fig.6 The dependence of cost function F on parameter combination (r1, r2) in Manning-type relationship： the enlarged view of the rectangular area in Fig.5

通過比較F值可以發(fā)現(xiàn), 采用廣義Manning型摩擦關(guān)系得到的結(jié)果優(yōu)于Chezy型摩擦關(guān)系。對比用兩種摩擦關(guān)系所得的同潮圖也可以看出, 采用廣義Manning型摩擦關(guān)系結(jié)果更合理。例如對M2分潮, 由于Chezy型關(guān)系給出的摩擦系數(shù)在深海區(qū)偏小, 計算結(jié)果與驗潮站實測值比較在海域東南部, 特別是濟(jì)州島附近和朝鮮半島西南沿岸顯著偏大; 在渤海由于Chezy型關(guān)系給出的摩擦系數(shù)偏大, 模擬得出的兩個無潮點可能離岸太近甚至消失(例如見Lefèvreet al,2000; 朱學(xué)明等, 2012)。因此下文的分析均采用由最優(yōu)參數(shù)[即(r1,r2)=(0.0009, 0.001375)]廣義Manning型摩擦關(guān)系計算所得的結(jié)果。

3 模擬結(jié)果驗證

本文選取的43個代表性實測站大體均勻地分布在研究區(qū)域內(nèi), 我們在計算網(wǎng)格中找到與實測站距離最近的節(jié)點, 比較計算和實測所得的潮汐調(diào)和常數(shù), 見表2(表中H代表振幅, 單位為cm;G代表格林尼治遲角, 單位為度(°)。由表可以得出, M2、S2、K1和 O1分潮振幅/遲角偏差絕對值的平均值分別為11.4cm/6.4°, 4.3cm/6.1°, 4.3cm/6.1°和 1.1cm/3.7°?？傮w上模式和實測結(jié)果具有較好的一致性。

表2 四個主要分潮計算和實測結(jié)果對比表(H: 振幅, G: 遲角)Tab.2 Results of computation and observation on harmonic constants for four principal constituents (H： amplitude, G： phase-lag)

4 渤、黃、東海海域的潮波分布

4.1 潮汐分布

對數(shù)值模擬的水位場進(jìn)行調(diào)和分析, 得到 M2、S2、K1和 O1四個分潮的調(diào)和常數(shù)。根據(jù)這些調(diào)和常數(shù)繪制的同潮圖分別示于圖7、圖8、圖9和圖10。對于渤、黃、東海潮汐、潮流的觀測和研究工作已有很多, Fang(1986)根據(jù)大量的實測資料以及數(shù)值模擬結(jié)果繪制了一份較完整的潮汐、潮流圖; 王凱等(1999)曾經(jīng)用三維數(shù)值模式模擬了該區(qū)域的 M2分潮的潮汐、潮流(王凱等, 1999); Fang等(2004)基于10年的衛(wèi)星高度計資料給出了該區(qū)域的潮汐同潮圖(Fang et al,2004); 王永剛等(2004)用驗潮站資料做了渤、黃、東海的同化數(shù)值模式; 朱學(xué)明等(2012)用 FVCOM 海洋數(shù)值模式模擬了該區(qū)域的潮汐、潮流。

圖7 M2分潮振幅(cm)和格林威治遲角(°)Fig.7 The amplitude (cm) and Greenwich phase-lag(°) of M2

由M2分潮的同潮圖(圖7)可以看到, 本文得到的黃河口外海無潮點的位置上與 Lefèvre等(2000)及朱學(xué)明等(2012)的結(jié)果有區(qū)別。后兩項研究給出的無潮點已退化到了陸地上, 而本文的結(jié)果顯示無潮點不但存在, 而且離岸有一定距離。這主要是由于黃河口岸形的改變引起的, 與王永剛等(2014)模擬得出的新近情況一致。其他三個無潮點(渤海西北部、山東半島東側(cè)海域和黃海南部)的位置與前人結(jié)果無太大差別。此外, 在臺灣島東北有一個明顯的退化的無潮點,這與前人結(jié)果也一致。

圖8是模擬得到的S2分潮同潮圖, 與M2分潮一樣, 一共有四個無潮點。與Fang等(2004)的結(jié)果相比,在本文中渤海內(nèi)部的兩個無潮點完全形成, 而后者基本在海岸線上, 其他地方等遲角線差別不大。關(guān)于S2分潮的等振幅線在蘇北外海的形狀, Fang等(2004)的結(jié)果中等振幅線有一個明顯的北向凸起, 本文則沒有。但是本文的結(jié)果與王永剛等(2004)和朱學(xué)明等(2012)的結(jié)果非常接近, 這可能是由于Fang等(2004)的研究中所用的衛(wèi)星高度計地面軌道較稀疏, 內(nèi)插引入的誤差造成的。從半日潮來看, 本文的結(jié)果與Fang等(2004)實測的結(jié)果相比, 除在渤海海域差別較大, 其他地方基本一致。

圖8 S2分潮振幅(cm)和格林威治遲角(°)Fig.8 The amplitude (cm) and Greenwich phase-lag (°) of S2

圖9和圖10給出了K1和O1分潮的振幅和格林威治遲角分布。它們與Fang等(2004)結(jié)果非常接近?？傮w來說, 與半日潮相比, 全日潮在該海域結(jié)構(gòu)較簡單, 振幅較小, K1分潮最大 40cm 左右, O1最大在30cm左右, 因此模擬結(jié)果較好。

圖9 K1分潮振幅(cm)和格林威治遲角(°)分布圖Fig.9 The amplitude (cm) and Greenwich phase-lag (°) of K1

圖10 O1分潮振幅(cm)和格林威治遲角(°)Fig.10 The amplitude (cm) and Greenwich phase-lag (°) of O1

4.2 潮能通量和潮汐傳播

一個潮周期內(nèi)單位時間通過自海底至海面單位寬度斷面的潮能通量叫做能通量密度, 其算式為

其中Φx和Φy分別為潮能通量密度的東向和北向分量,H和G分別為水位的調(diào)和常數(shù), U、ξ和V、η分別為潮流東向和北向分量的調(diào)和常數(shù)。

本文根據(jù)式(18)和式(19)計算了四個分潮的潮能通量, 其中最大的半日分潮M2和最大的全日分潮K1的能通量分布見圖11和12。為了顯示潮能通量區(qū)域的特征, 圖中矢量的長度與能通量的平方根成正比例, 其大小比例尺見圖的左下方。

從圖11的M2分潮來看, 潮能輸送路徑主要分為三支： 第一支是通過東海陸架邊緣南側(cè)進(jìn)入東海后反轉(zhuǎn), 由臺灣海峽輸入南海; 第二支是通過東海陸架邊緣北側(cè)經(jīng)東海進(jìn)而到達(dá)黃海, 在蘇北外海和山東半島南部海域形成一個逆時針旋轉(zhuǎn)的潮波系統(tǒng)進(jìn)而耗散; 第三支是沿著朝鮮半島西海岸一直向北進(jìn)入渤海耗散。從圖可以看出, 大部分潮能是在黃海消耗的, 進(jìn)入渤海的潮能較小; 同時在渤海海峽處, 半日潮能通量均向西, 即使在海峽南部也看不到向東的反射波能通量, 這與全日潮波有明顯不同(見下)。S2分潮能通量分布特征與 M2基本相同; 量值較小, 只有M2的10%—20%。

與半日潮相比, 全日潮能量較小(其中 K1分潮見圖12)。全日潮的輸送路徑分三支： 一支由臺灣島東側(cè)由東海反轉(zhuǎn)進(jìn)入臺灣海峽; 另一支在黃海中部向西北方向前進(jìn), 沿山東半島南部沿岸和蘇北沿岸向南發(fā)生耗散; 還有一支在濟(jì)州島南部沿著朝鮮半島西岸向北前進(jìn), 通過渤海海峽北部進(jìn)入渤海耗散大部分, 其余部分通過渤海海峽南部回到黃海, 并在山東半島北側(cè)外海消耗殆盡。全日潮波南向能通量基本上可以到達(dá)長江口北岸, 而半日潮波南向能通量只能到達(dá)蘇北外海中部。這說明從反射波相對于入射波的能量來看, 全日潮要比半日潮強(qiáng)。

圖11 M2分潮能通量密度分布圖Fig.11 The energy flux density vectors of M2

圖12 K1分潮能通量密度分布圖Fig.12 The energy flux density vectors of K1

5 結(jié)語

基于FVCOM海洋數(shù)值模式, 建立了渤、黃、東海陸架海域的潮汐數(shù)值模型, 分別采用Chezy型和廣義Manning型底摩擦系數(shù)模擬了潮汐、潮流, 并通過與實測數(shù)據(jù)對比確定了其最優(yōu)參數(shù)。將計算結(jié)果與實測數(shù)據(jù)進(jìn)行比較, 結(jié)果良好。最后, 用廣義 Manning型底摩擦系數(shù)最優(yōu)結(jié)果對該海域的潮汐、潮流和潮能進(jìn)行了分析。主要結(jié)論有：

(1) 采用廣義Manning型比采用Chezy型底摩擦系數(shù)得到的結(jié)果更好, 因此建議在模擬渤、黃、東海潮汐時最好采用前者。

(2) 本研究得到的底摩擦系數(shù)在0.009至0.0014之間, 顯著低于 Proudman(1953)得到的 0.0026, 這也表明渤、黃、東海海底摩擦系數(shù)要比大多數(shù)模擬研究中采用的數(shù)值要小。

(3) 原始的Manning公式的冪值為負(fù)數(shù), 而本研究得到的冪值為正。這表明在渤、黃、東海陸架, 總體上水越淺摩擦系數(shù) r值越小, 也表明原始的Manning公式必須推廣才能適用于渤、黃、東海大范圍潮汐模擬。

王 凱, 方國洪, 筰馮士, 1999. 渤海、黃海、東海M2潮汐潮流的三維數(shù)值模擬. 海洋學(xué)報, 21(4)： 1—13

王永剛, 方國洪, 曹德明等, 2004. 渤、黃、東海潮汐的一種驗潮站資料同化數(shù)值模式. 海洋科學(xué)進(jìn)展, 22(3)： 253—274

朱學(xué)明, 鮑獻(xiàn)文, 宋德海等, 2012. 渤、黃、東海潮汐、潮流的數(shù)值模擬與研究. 海洋與湖沼, 43(6)： 1103—1113

吳 頔, 方國洪, 崔欣梅等, 2015. 泰國灣及鄰近海域潮汐潮流的數(shù)值模擬. 海洋學(xué)報, 37(1)： 11—20

林美華, 方國洪, 1991. 中國海標(biāo)準(zhǔn)經(jīng)緯度水深和基準(zhǔn)面數(shù)據(jù)表. 青島： 中國科學(xué)院海洋研究所, 1—144

周 朦, 方國洪, 1987. |U|U的Fourier展開和渤海海底拖曳系數(shù) CD. 海洋與湖沼, 18(1)： 1—11

Chen C S, Beardsley R C, Cowles G, 2006. An Unstructured Grid,Finite-Volume Coastal Ocean Model-FVCOM User Manual,Technical Report SMAST/UMASSD-06-0602. 2nd edn. New Bedford： School for Marine Science and Technology,University of Massachusetts Dartmouth, 6—8

Cheng Y C, Andersen O B, 2011. Multimission empirical ocean tide modeling for shallow waters and polar seas. Journal of Geophysical Research Oceans, 116(C11), doi： 10.1029/2011JC007172

Dronkers J J, 1964. Tidal Computations in Rivers and Coastal Waters. Amsterdam, North-Holland： Interscience Publishers,518

Fang G H, 1986. Tide and tidal current charts for the marginal seas adjacent to China. Chinese Journal of Oceanology and Limnology, 4(1)： 1—16

Fang G H, Li S D, Cao D M, 1987. Three methods for estimating turbulent stress and drag coefficient in tidal currents of the Hangzhou Bay. Science Bulletin, 32(4)： 252—257

Fang G H, Wang Y G, Wei Z X et al, 2004. Empirical cotidal charts of the Bohai, Yellow, and East China Seas from 10 years of TOPEX/Poseidon altimetry. Journal of Geophysical Research： Oceans, 109(C11), doi： 10.1029/2004JC002484

Jan S, Chern C S, Wang J, 2002. Transition of tidal waves from the East to South China Seas over the Taiwan Strait：Influence of the abrupt step in the topography. Journal of Oceanography, 58(6)： 837—850

Lefèvre F, Le Provost C, Lyard F H, 2000. How can we improve a global ocean tide model at a regional scale? A test on the Yellow Sea and the East China Sea. Journal of Geophysical Research： Oceans, 105(C4)： 8707—8725

Pawlowicz R, Beardsley B, Lentz S, 2002. Classical tidal harmonic analysis including error estimates in MATLAB using T_TIDE. Computers & Geosciences, 28(8)： 929—937 Proudman J, 1953. Dynamical Oceanography, London： Methuen and Co, 409

Taylor G I. 1920. Tidal friction in the Irish Sea. Philosophical Transactions of the Royal Society of London. Series A, 220：1—33

Wang Y H, Fang G H, Wei Z X et al, 2012. Cotidal charts and tidal power input atlases of the global ocean from TOPEX/Poseidon and JASON-1 altimetry. Acta Oceanologica Sinica,31(4)： 11—23