印度尼西亞近海潮汐潮流的數(shù)值模擬

滕 飛，方國洪,王新怡，魏澤勛，王永剛

(1. 國家海洋局 第一海洋研究所，山東 青島 266061;2. 海洋環(huán)境科學(xué)和數(shù)值模擬國家海洋局重點實驗室，山東 青島 266061)

潮汐是海洋中最重要的現(xiàn)象之一，而印度尼西亞的潮汐現(xiàn)象在海洋中是最復(fù)雜的。復(fù)雜的沿岸地形、狹窄的海峽和無數(shù)的小島嶼、廣闊的大陸架以及迅速變化的水深，再加上從印度洋和太平洋傳入的巨大的潮汐能量，這些組成了一個復(fù)雜的潮汐系統(tǒng)。印尼海中有多個無潮點，海峽的潮流流速大。潮汐對印尼海的海水混合和環(huán)流有很大影響。

早在1981年，Webb[1]模擬了卡奔塔利亞灣和阿拉弗拉海的潮汐；1993年，Wolanski[2]模擬了卡奔塔利亞灣的潮汐；1996年，Hatayama等[3]計算了整個印尼海的M2和K1分潮。2002年，Egbert等[4]將衛(wèi)星高度計同化到動力學(xué)模式中，給出了印尼海的正壓潮。根據(jù)這個結(jié)果，于2005年，Ray[5]對印尼海域潮汐作了詳細的闡述。2008年，Koropitan和Ikeda[6]模擬了爪哇海的潮汐環(huán)流及混合。

前人對于印尼海及其周邊海域的水動力模擬研究均采用矩形網(wǎng)格進行計算。采用目前較為流行的有限體積近岸海洋模式FVCOM[7]，該模式采用非結(jié)構(gòu)三角形網(wǎng)格，對印尼海主要的島嶼、海峽等地形變化劇烈地方進行加密處理，大大提高了海洋模式對復(fù)雜海岸線的擬合程度。模式采用干-濕網(wǎng)格判別技術(shù)較好地處理了模型的固體邊界，考慮地轉(zhuǎn)科氏力、非線性效應(yīng)、水平對流擴散等，建立了分辨率比較高的印尼海及其周邊海域的數(shù)值模式。模擬結(jié)果與前人的結(jié)果相比較有所改進。根據(jù)模式結(jié)果，本文給出了印尼海主要的4個分潮(M2，S2，K1，O1)的同潮圖，并分析了印尼海及其周邊海域的潮汐、潮流的分布特征和潮波的傳播特征。

1 計算方法與計算過程基本方程組

1.1 基本方程組

FVCOM模式采用σ垂向坐標(biāo)來模擬不規(guī)則的底部地形，坐標(biāo)變換公式如下：

(1)

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

(2)

(3)

(4)

(5)

式中，t為時間；r為地球半徑；λ為計算點經(jīng)度；φ為計算點緯度；σ為垂向分層(海表面為0，海底為-1)；H為未擾動水深；ζ為相對于未擾動海面的潮汐高度；D為擾動水深；u，v，ω分別為東向，北向和垂直速度；ρ為海水密度；f為科氏參數(shù)；g為地球重力加速度；Km為垂直混合參數(shù)；Fu和Fv為水平動量的擴散項。

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

βiAicos2φcos(ωit+2λ)

(6)

對于全日潮，引潮勢為

βiAisin2φcos(ωit+λ)

(7)

式中，Ai和ωi為第i個分潮的振幅和頻率；Nsemi和Ndiurnal為模式中參與計算的半日潮和全日潮的個數(shù)；βi可表示成：

βi=1+klove-hlove

(8)

式中，k和h為love數(shù)。在FVCOM中，各個分潮的klove和hlove所采用的數(shù)值見表1。

表1 模式中各分潮love數(shù)klove和hloveTable 1 Love numbers klove and hlove for the tidal constituents in model

1.2 模式水深的確定

在潮汐模擬中，水深是至關(guān)重要的參數(shù)。但是我們發(fā)現(xiàn)，現(xiàn)有的格點化水深數(shù)據(jù)etopo1(1′分辨率)和etopo5(5′分辨率)在印尼近海深水區(qū)一致性較好，而在淺海區(qū)則差別較大，而且兩者與海圖資料都不相符。我們通過多次試驗，最后確定了本研究采用的水深數(shù)據(jù)。在深海等大部分海域采用etopo1的數(shù)據(jù)，在納土納海和爪哇海采用從海圖讀取的水深數(shù)據(jù)，在阿拉弗拉海從131°E到135°E用etopo1和etopo5的數(shù)據(jù)進行過渡，在卡奔塔利亞灣中采用etopo5的數(shù)據(jù)，模式中采用的水深分布見圖1。

MK-望加錫海峽，ML-馬魯古海，HH-哈馬黑拉海，F(xiàn)L-弗洛勒斯海，SW-薩武海，TM-帝汶海，AF-阿拉弗拉海，CP-卡奔塔利亞灣

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

本研究海區(qū)為從16°S～8°N，100°～145°E的印尼近海(圖1)，但數(shù)值計算區(qū)域較大，包括印度尼西亞內(nèi)部和周邊海域、南海大部(包括泰國灣)、印度洋東北部、太平洋西南部以及澳大利亞北部的卡奔塔利亞灣(圖2)。計算網(wǎng)格水平分辨率在印尼海域島嶼平均為1/12度，在大陸邊界平均為1/5度，在開邊界平均為1/2度。共包含80 096個網(wǎng)格節(jié)點，155 977個三角單元(圖2)；垂向分6層。FVCOM模型采用內(nèi)外模分離的方式求解，二維外模數(shù)值格式是基于三角形網(wǎng)格的有限體積法，將連續(xù)方程和動量方程在三角形單元內(nèi)積分后，通過改進的四階龍格庫塔方法求解。三維內(nèi)模的動量方程采用簡單的顯式和隱式相結(jié)合的差分格式求解，時間步長外模為6 s，內(nèi)外模時間步長比率為10。

計算的初始條件假設(shè)海洋是靜止的，海表面的擾動水位初始全部為0，所有三角形網(wǎng)格中心點的流速u=v=0。本文采用正壓模式，取整個海域內(nèi)的溫度和鹽度均為常數(shù)，溫度為18 ℃，鹽度為35.0。在南海、印度洋、太平洋共設(shè)有4個開邊界，開邊界上每個點的水位給定，由T_tide的預(yù)報程序給出[8]，其中M2，S2，K1和O1調(diào)和常數(shù)來自DTU10[9]數(shù)據(jù)。

圖2 計算網(wǎng)格Fig.2 Computation grids

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

2 計算結(jié)果與實測結(jié)果的比較

2.1 評價方法

為了比較模式的結(jié)果，我們用TOPEX/Poseidon(T/P)衛(wèi)星數(shù)據(jù)和地面驗潮站數(shù)據(jù)分別與模式結(jié)果進行比較驗證。下面介紹幾種比較常用的計算值與觀測值之偏差不同的度量方法[10]。首先是一種較直觀的方法，分別是振幅之間和遲角之間的均方偏差：

(9)

(10)

式中,H和G分別為分潮的振幅和遲角，下標(biāo)a和b分別為模式結(jié)果和觀測結(jié)果；K為比較的站位總個數(shù)。這種方法比較直觀，但是當(dāng)H較小時，G的計算和觀測值都不穩(wěn)定，同時,H比較小時,G的誤差對潮高誤差影響也較小，因此在無潮點附近遲角很容易出現(xiàn)較大的誤差，因此不考慮H的差別而將所有遲角差值的平方進行平均具有相當(dāng)大的偶然性。另一種方法是考察潮高之差的均方根。觀測站潮高的均方根值為

(11)

觀測潮高和模式潮高之差的均方根為

(12)

相對偏差：

k=RMSd/RMSh

(13)

2.2 與T/P衛(wèi)星數(shù)據(jù)的比較

為對模式結(jié)果的準確度進行評估，我們選用T/P衛(wèi)星軌道交叉點的調(diào)和常數(shù)對模式結(jié)果進行驗證。由于在陸地邊緣的結(jié)果有誤差，我們對落在模式計算范圍內(nèi)的所有T/P交叉點數(shù)據(jù)中剔除了幾個島嶼附近的數(shù)據(jù)。實際采用的比較站位共104個，其分布圖見圖3。表2為模式數(shù)據(jù)與T/P衛(wèi)星數(shù)據(jù)的振幅之間和遲角之間的均方偏差，潮高之差的均方根以及相對偏差。

圖4和圖5分別給出在T/P交叉點上計算和T/P觀測所得的M2和K1分潮振幅和遲角的相關(guān)關(guān)系圖。圖中RMS值與表2相同。r為模式值和觀測值的相關(guān)系數(shù)。

圖3 T/P軌道交叉點站位Fig.3 Positions of the TOPEX/Poseidon crossover points

表2 模式結(jié)果與T/P衛(wèi)星軌道交叉點觀測值比較Table 2 Comparison of the model results to the observations at TOPEX/Poseidon crossover points

圖4 M2分潮模式結(jié)果和衛(wèi)星數(shù)據(jù)的相關(guān)性Fig.4 Correlation between the model results and the T/P data for M2

圖5 K1分潮模式結(jié)果和衛(wèi)星數(shù)據(jù)的相關(guān)性Fig.5 Correlation between the model results and the T/P data for K1

2.3 與驗潮站數(shù)據(jù)比較

本文的驗潮站取自國際水文組織IHO潮汐調(diào)和常數(shù)數(shù)據(jù)集，共選取其中79個站位(圖6)。挑選的原則是用于分析調(diào)和常數(shù)的觀測時間大于25 d，并且較均勻分布在研究區(qū)域內(nèi)。用驗潮站數(shù)據(jù)與模式結(jié)果進行了對比，比較結(jié)果見圖7、圖8。模式數(shù)據(jù)與驗潮站數(shù)據(jù)的振幅之間和遲角之間的均方偏差，潮高之差的均方根以及相對偏差見表3。

圖6 驗潮站站位Fig.6 The positions of the tide gauge stations

表3 模式數(shù)據(jù)與驗潮站數(shù)據(jù)比較Table 3 Comparison of the model results to the data from tide gauge

圖7 M2分潮模式結(jié)果和驗潮站數(shù)據(jù)的相關(guān)系數(shù)Fig.7 Correlation between the model results and the tidal gauge data for M2

圖8 K1分潮模式結(jié)果和驗潮站數(shù)據(jù)的相關(guān)系數(shù)Fig.8 Correlation between the model results and the tidal gauge data for K1

從表2和表3來看，模式結(jié)果與T/P結(jié)果比較一致，振幅均方根差大約為2～5 cm，遲角均方根偏差大約為5°～8°，潮高的均方根偏差大約為2～4 cm，相對偏差為0.11～0.19。模式結(jié)果與驗潮站結(jié)果的偏差大一點，振幅均方根偏差為4～9 cm，遲角均方根偏差為9°～12°，潮高的均方根偏差為4～8 cm，相對偏差為0.2～0.3之間。模式和驗潮站觀測結(jié)果偏差較大是由于驗潮站大多位于港口內(nèi)部，局地效應(yīng)不能很好地在大范圍數(shù)值模式中考慮。但總的來看模式和觀測值有較好的一致性。特別與Hatayama等[3]的圖4和圖14比較，可以看到本文結(jié)果更符合觀測值。

3 印尼海及其周邊海域的潮波分布

3.1 潮汐分布

模式計算所得到的4個分潮的同潮圖見圖9。M2分潮最大振幅出現(xiàn)在澳大利亞西北海域，其值可超過3 m。次大值出現(xiàn)在加里曼丹島西北岸，其值超過1.5 m，與Fang等[11]的南海潮汐數(shù)值模擬結(jié)果一致。印尼海大部海域M2分潮振幅都不大，在0.6 m以下。在澳大利亞的西北海域約瑟夫-波拿巴灣西北口有一個順時針的無潮點，在阿拉弗拉海的北部有一逆時針無潮點。在南納土納海、爪哇海到弗洛勒斯海一帶M2存在多個無潮點，振幅非常小。

……為振幅/cm；——為遲角/°

S2分潮的分布特征與M2相近(圖9a,b)。在南納土納海、爪哇海至弗洛勒斯海一帶也與M2一樣存在多個無潮點，振幅很小，但是無潮點的位置與M2有些不同。另外，值得特別注意的是，在望加錫海峽，M2的遲角從南向北推遲，而S2的遲角反過來從北向南推遲。

K1分潮振幅較大值出現(xiàn)在南海南部、納土納海南部、爪哇海東北部、阿弗勒斯海北部和約瑟夫·波拿巴灣，其中尤以阿弗勒斯海北部的新幾內(nèi)亞西南沿岸為最大，可超過1.2 m(圖9c)。卡奔塔利亞灣有一個順時針旋轉(zhuǎn)的無潮點，西爪哇海有一個同潮時線密集的波節(jié)帶，但似乎未能形成無潮點。

O1分潮的分布特征與K1相近，其振幅大約為K1的2/3左右?？ū妓麃啚晨谝灿幸粋€無潮點，西爪哇海也存在一個波節(jié)帶，這里是否形成無潮點還難以確定。

圖10為潮汐類型圖。由于印尼海域半日潮波結(jié)構(gòu)復(fù)雜，一些區(qū)域(特別是M2無潮點附近)S2分潮振幅可能大于M2，本文采用F=(HO1+HK1)/(HM2+HS2)的計算結(jié)果作為潮汐類型的判斷參數(shù)。按照國外的習(xí)慣，F(xiàn)=0～0.25為半日潮類型，F(xiàn)=0.25～1.5為半日為主的混合潮汐類型，F(xiàn)=1.5～3為全日潮為主的混合潮類型，F(xiàn)>3.0為全日潮類型。從圖中可以看出,在馬來半島以東海域、納土納海南部(包括卡里馬塔海峽)、爪哇海東部、新幾內(nèi)亞島西南外海及卡奔塔利亞灣東南部為規(guī)則全日潮。其它海域大部分為混合潮。半日潮類型只出現(xiàn)在澳大利亞西北近岸海域。

圖10 潮汐類型分布圖Fig.10 Distributions of tide types

3.2 潮流分布

為了研究印尼海區(qū)域的潮流分布情況，分別計算了4個分潮的潮流橢圓要素(長軸、短軸、傾斜角、格林威治遲角)，其中M2和K1分潮的潮流橢圓分別見圖11和圖12，陰影海域部分代表潮流以逆時針方向旋轉(zhuǎn)，無陰影海域部分為順時針旋轉(zhuǎn)。

對于M2分潮來說，在連接太平洋和爪哇海東北部的望加錫海峽、馬魯古海、哈馬黑拉海以及斯蘭海潮流呈現(xiàn)往復(fù)流，流速比較大，其中在馬魯古海與班達海之間的個別海峽中最大流速可以超過100 cm/s。而在班達海和蘇拉威西海中部的深水區(qū)域潮流非常小，都小于11 cm/s。在爪哇海的內(nèi)部水深較淺，流速較強，其中在爪哇海東南部最大可以達到51 cm/s?？ū妓麃啚乘钶^淺，在該灣東北部流速最大可達80 cm/s，這個區(qū)域的S2分潮與M2相比除流速偏小以外，分布特征基本一樣。

圖11 M2分潮潮流橢圓Fig.11 M2 current ellipses

對于K1分潮而言，在納土納海南部及卡里馬塔海峽附近的流速呈現(xiàn)往復(fù)流，流速較大，最大流速出現(xiàn)在卡里馬塔海峽西部海峽中。在爪哇海呈現(xiàn)東西的往復(fù)流，流速大約在30 cm/s左右，最大流速出現(xiàn)在爪哇海的西北部。在班達海和蘇拉威西海水深較深地方的流速較小。在望加錫海峽、馬魯古海、哈馬黑拉海以及斯蘭海中流速比M2分潮的流速小。在卡奔塔利亞灣的K1分潮流基本上都按順時針方向旋轉(zhuǎn)，特別在灣口，旋轉(zhuǎn)特性明顯。O1分潮的分布情況與K1分潮相似，只是流速較小。

圖12 K1分潮潮流橢圓Fig.12 K1 current ellipses

為了對印尼近海潮流的大小有一個總體的了解，圖13給出了4個主要分潮最大流速之和的分布。由圖可以看出，強的潮流主要出現(xiàn)在陸架淺水區(qū)和狹窄的海峽中。

圖13 M2、S2、K1和O1分潮最大流速之和(單位：cm/s)Fig.13 Sum of the maximum current speeds of M2, S2, K1 and O1 (unit: cm/s)

3.3 潮余流分布

由4個分潮共同引起的歐拉潮余流分布圖見圖14，在卡奔塔利亞灣中潮余流呈現(xiàn)非常規(guī)律的順時針環(huán)流。在澳大利亞西北部潮余流很強，達到將近4 cm/s。在納土納海西部沿岸區(qū)域存在比較弱的由北向南的余流。

圖14 潮余流分布圖Fig.14 Distributions of tidal residual currents

圖15 M2，S2，K1和O1分潮能通量密度分布圖Fig.15 The energy flux density vectors of tidal constituents M2, S2, K1 and O1

4 潮能通量和潮波傳播

潮能通量強度又叫做能通量密度，即單位時間通過自海底至海面單位寬度斷面的潮能通量，其計算公式為

(14)

式中，v和?分別為流速和能通強度向量；T為潮波周期；ρ為海水密度，本文取1 025 kg/m3；h為當(dāng)?shù)厮?。?14)對各分潮求時間平均，可得[12]：

ρghHUcos(G-ξ)

(15)

(16)

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

Ray[5]曾給出印尼海M2和K1分潮潮能通量的輸送特征。由于他采用了統(tǒng)一的線性比例尺，許多潮能通量較弱的海區(qū)如納土納海和爪哇海都未能給出。本文根據(jù)式(15)和(16)計算了M2，S2，K1和O1分潮的潮能通量，見圖15。為了顯示潮能通量區(qū)域的特征，圖中矢量的長度取與能通量的平方根成正比例，其大小比例尺見圖的左下方。

從半日潮M2分潮來看，潮能主要是由印度洋進入印尼海和卡奔塔利亞灣，其中部分能量穿過望加錫海峽、馬魯古海和哈馬黑拉海后進入太平洋。南海有能量向南傳輸，但只能影響到納土納海，爪哇海的M2分潮主要還是依靠印度洋來的能量來維持。S2分潮能量在大部分區(qū)域與M2分潮相同，也是由南向北輸向太平洋。與M2分潮不同的是在望加錫海峽潮能由北向南進入印尼海域，而M2分潮的潮能輸送方向是由南向北進入蘇拉威西海，進而進入太平洋。

在印尼海域全日潮能量傳輸與半日潮正好相反，K1和O1分潮的潮能是由太平洋通過印尼各海域輸向印度洋和卡奔塔利亞灣。來自南海的潮能基本上只到達卡里馬塔海峽附近，爪哇海的全日潮運動主要依靠來自太平洋，通過望加錫海峽的潮能來維持。

5 結(jié) 語

基于FVCOM海洋數(shù)值模式，本研究采用非結(jié)構(gòu)三角形網(wǎng)格對計算海域進行分辨率的調(diào)整，應(yīng)用有限體積法對原始方程進行離散求解，建立了印尼海及其周邊海域的潮汐、潮流數(shù)值模型，并與T/P衛(wèi)星數(shù)據(jù)和驗潮站數(shù)據(jù)進行比較，結(jié)果良好。用該模式結(jié)果進行了印尼海域的潮汐類型，潮流和潮能的分析。

與Hatayama等[3]的結(jié)果相比，本文結(jié)果與實測符合程度有顯著提高。但是由于印尼近海島嶼眾多，潮波結(jié)構(gòu)復(fù)雜，本研究尚有待進一步改進。一個有效的方法是將印尼近海劃分為若干個海區(qū)分別進行模擬。

參考文獻(References):

[1] WEBB D J. Numerical model of the tides in the Gulf of Carpentaria and the Arafura Sea[J]. Marine & Freshwater Research, 1981, 32(1): 31-44.

[2] WOLANSKI E. Water circulation in the Gulf of Carpentaria [J]. Journal of Marine Systems, 1993, 4(5): 401-420.

[3] HATAYAMA T, AWAJI T, AKITOMO K. Tidal currents in the Indonesian Seas and their effect on transport and mixing [J]. Journal of Geophysical Research, Ocean,1996, 101(C5):12353-12373

[4] EGBERT G D, EROFEEVA S Y. Efficient inverse modeling of barotropic oceantides [J]. Journal of Atmospheric and Oceanic Technology, 2002, 19(2):183-204.

[5] RAY R D, EGBERT G D, EROFEEVA S Y.A brief overview of tides in the Indonesian Sea [J]. Oceanography, 2005,18(4):74-79.

[6] KOROPITAN A F, IKEDA M. Three-dimensional modeling of tidal circulation and mixing over the Java Sea [J]. Journal of Oceanography, 2008, 64(1): 61-80.

[7] CHEN C S, BEARDSLEY R C, COWLES G. An unstructured grid, finite-volume coastal ocean model FVCOM user manual [R]. North Dartmouth: UMASSD, 2006:6-8.

[8] PAWLOWICZ R,BEARDSLEY B,LENTZ S. Classical tidal harmonic analysis including error estimates in MATLAB using T_TIDE [J]. Computers and Geosciences, 2002, 28(8):929-937.

[9] CHENG Y C, ANDERSEN O B.Multimission empiricalocean tide modeling for shallow waters and polar seas [J]. Journal of Geophysical Research, Ocean(1978-2012),2011, 116(C1):1001-1011.

[10] WANG Y H, FANG G H, WEI Z X, et al. Accuracy assessment of global ocean tide models base on satellite altimetry [J]. Advances in Earth Science, 2010, 25(4):353-359.

[11] FANG G H, KWOK Y K, YU K, et al. Numerical simulation of principal tidal constituents in the South China Sea, Gulf of Tonkin and Gulf of Thailand [J]. Continental Shelf Research, 1999, 19(7):845-869.

[12] FANG G H. A Two-Dimensional Numerical Model for Tidal Motioninthe Taiwan Strait [J]. Marine Geophysical Research, 1984, 7(1-2): 267-276.