印度尼西亞近海潮汐潮流的數值模擬

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

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

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

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

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

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

1.1 基本方程組

FVCOM模式采用σ垂向坐標來模擬不規則的底部地形，坐標變換公式如下：

(1)

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

(2)

(3)

(4)

(5)

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

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

βiAicos2φcos(ωit+2λ)

(6)

對于全日潮，引潮勢為

βiAisin2φcos(ωit+λ)

(7)

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

βi=1+klove-hlove

(8)

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

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

1.2 模式水深的確定

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

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

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

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

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

圖2 計算網格Fig.2 Computation grids

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

2 計算結果與實測結果的比較

2.1 評價方法

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

(9)

(10)

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

(11)

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

(12)

相對偏差：

k=RMSd/RMSh

(13)

2.2 與T/P衛星數據的比較

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

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

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

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

圖4 M2分潮模式結果和衛星數據的相關性Fig.4 Correlation between the model results and the T/P data for M2

圖5 K1分潮模式結果和衛星數據的相關性Fig.5 Correlation between the model results and the T/P data for K1

2.3 與驗潮站數據比較

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

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

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

圖7 M2分潮模式結果和驗潮站數據的相關系數Fig.7 Correlation between the model results and the tidal gauge data for M2

圖8 K1分潮模式結果和驗潮站數據的相關系數Fig.8 Correlation between the model results and the tidal gauge data for K1

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

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

3.1 潮汐分布

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

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

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

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

O1分潮的分布特征與K1相近，其振幅大約為K1的2/3左右。卡奔塔利亞灣口也有一個無潮點，西爪哇海也存在一個波節帶，這里是否形成無潮點還難以確定。

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

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

3.2 潮流分布

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

對于M2分潮來說，在連接太平洋和爪哇海東北部的望加錫海峽、馬魯古海、哈馬黑拉海以及斯蘭海潮流呈現往復流，流速比較大，其中在馬魯古海與班達海之間的個別海峽中最大流速可以超過100 cm/s。而在班達海和蘇拉威西海中部的深水區域潮流非常小，都小于11 cm/s。在爪哇海的內部水深較淺，流速較強，其中在爪哇海東南部最大可以達到51 cm/s。卡奔塔利亞灣水深較淺，在該灣東北部流速最大可達80 cm/s，這個區域的S2分潮與M2相比除流速偏小以外，分布特征基本一樣。

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

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

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

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

圖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，在卡奔塔利亞灣中潮余流呈現非常規律的順時針環流。在澳大利亞西北部潮余流很強，達到將近4 cm/s。在納土納海西部沿岸區域存在比較弱的由北向南的余流。

圖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為當地水深。式(14)對各分潮求時間平均，可得[12]：

ρghHUcos(G-ξ)

(15)

(16)

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

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

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

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

5 結 語

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

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

參考文獻(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.