徐振華, 雷方輝, 婁安剛, 曹圣山

(1. 中國海洋大學(xué) 數(shù)學(xué)系, 山東 青島 266100; 2. 中海石油研究中心, 北京 100027; 3. 中國海洋大學(xué) 環(huán)境科學(xué)與工程學(xué)院, 山東 青島 266100; 4. 中國海洋大學(xué) 海洋環(huán)境與生態(tài)教育部重點(diǎn)實(shí)驗(yàn)室 山東 青島266003)

北部灣潮汐潮流的數(shù)值模擬

徐振華1, 雷方輝2, 婁安剛3,4, 曹圣山1

(1. 中國海洋大學(xué) 數(shù)學(xué)系, 山東 青島 266100; 2. 中海石油研究中心, 北京 100027; 3. 中國海洋大學(xué) 環(huán)境科學(xué)與工程學(xué)院, 山東 青島 266100; 4. 中國海洋大學(xué) 海洋環(huán)境與生態(tài)教育部重點(diǎn)實(shí)驗(yàn)室 山東 青島266003)

基于正交曲線坐標(biāo)的ECOMSED三維水動(dòng)力模式, 對(duì)北部灣的潮汐潮流進(jìn)行數(shù)值模擬。選用不同的海底摩擦系數(shù)、海底粗糙度系數(shù)以及水平湍流摩擦系數(shù)進(jìn)行數(shù)值試驗(yàn)。試驗(yàn)結(jié)果表明, 當(dāng)海底摩擦系數(shù)取 2×10-3～3×10-3, 海底粗糙度系數(shù)取 1×10-3～2×10-3m, 而水平湍流摩擦系數(shù)取 1×102～5×103m2/s時(shí), 模擬所得潮汐、潮流結(jié)果與實(shí)測數(shù)據(jù)吻合較好; 并由此對(duì)北部灣的潮汐、潮流和潮余流等特征進(jìn)行了分析。還給出了垂直湍流黏滯系數(shù)呈拋物型的分布特征。對(duì)該海域的水動(dòng)力狀況有了進(jìn)一步的了解, 為該海域的環(huán)境保護(hù)規(guī)劃提供了科學(xué)依據(jù)。

北部灣; 潮汐; 潮流; 數(shù)值模擬

北部灣是南海的主要海灣。夏華永等[1]在 1997年依經(jīng) Sigma坐標(biāo)變換后具有自由表面的三維非線性 Navier-Stokes方程, 用分裂算子法對(duì)北部灣 M2和 K1分潮進(jìn)行了三維潮流數(shù)值模擬, 孫洪亮等[2]在2001年利用 POM 三維水動(dòng)力模式模擬了北部灣的潮汐、潮流特征, 殷忠斌等[3]在 1996年采用不同的參數(shù)對(duì)北部灣的 K1分潮進(jìn)行了數(shù)值試驗(yàn), 雖然北部灣的數(shù)值模擬工作開展得較為廣泛, 但絕大部分工作都采用矩形網(wǎng)格, 并且網(wǎng)格都是均勻分布, 導(dǎo)致岸邊界附近的流向很難與實(shí)際相符。作者基于二階湍封閉模型的 ECOMSED三維水動(dòng)力模式, 在北部灣建立了一個(gè)較為細(xì)致的正交曲線坐標(biāo)系統(tǒng), 并且對(duì)灣口處的網(wǎng)格進(jìn)行加密, 垂向分為21層,可以更好地分析出各個(gè)物理量的垂直分布特征;采用不同的海底摩擦系數(shù)、海底粗糙度系數(shù)以及水平湍流摩擦系數(shù)進(jìn)行數(shù)值試驗(yàn), 并在此基礎(chǔ)上, 分析出該海區(qū)潮汐、潮流和潮余流等特征, 以期對(duì)北部灣的水動(dòng)力狀況有進(jìn)一步的了解, 并為制定該海域環(huán)境保護(hù)規(guī)劃提供科學(xué)依據(jù)。

1 動(dòng)力模式及其配置

作者所用的ECOMSED模式是由Blumberg等人在美國普林斯頓大學(xué)的三維海洋模式(POM)及其后來發(fā)展的河口、陸架和海洋模式(ECOM)的基礎(chǔ)上發(fā)展而來的, 是一個(gè)較為成熟的集海浪和沉積輸運(yùn)為一體的淺海三維水動(dòng)力學(xué)模式, 其中三維水動(dòng)力模塊從原始三維方程出發(fā), 以自由水位、三方向速度分量、溫度、鹽度、密度以及代表湍流的兩個(gè)特征量：湍動(dòng)能和湍宏觀尺度作為預(yù)報(bào)變量。關(guān)于此模塊的特點(diǎn)和有關(guān)方程及其差分格式在Blumberg[4]2002年的有關(guān)文章中有詳細(xì)敘述。

本文的計(jì)算范圍為105.5°～111.5°E,15.5°～22.0°N。模式水平方向采用正交曲線網(wǎng)格, 網(wǎng)格空間步長最大為5 km; 垂向不等距分為21層, 采用σ坐標(biāo)變換, 第一層取 0.007, 從而各物理量可以更好地逼近海表的物理場, 其余各層均勻分布。在模式中只考慮 M2和 K1分潮作為開邊界的潮強(qiáng)迫條件, 其調(diào)和參數(shù)從海洋圖集南海分冊的同潮圖內(nèi)插得到; 模擬采用零初始條件, 計(jì)算選取的內(nèi)模時(shí)間步長為134.631 s。為了保證計(jì)算的穩(wěn)定性, 強(qiáng)迫的邊界潮位從零開始逐步增加, 經(jīng)過 3個(gè)潮周期后達(dá)到正常變化, 第四個(gè)潮周期后形成穩(wěn)定的潮波。計(jì)算海區(qū)的水深分布及正交曲線計(jì)算網(wǎng)格分別見圖1和圖2。采用不同的海底摩擦系數(shù)、海底粗糙度系數(shù)以及水平湍流摩擦系數(shù)對(duì)北部灣的潮汐潮流進(jìn)行數(shù)值試驗(yàn)。

圖1 北部灣海底地形(單位：m)Fig. 1 The sea-floor topographic map in Beibu Bay(unit：m)

圖2 北部灣正交曲線計(jì)算網(wǎng)格Fig. 2 Orthogonal curvilinear computation grids in Beibu Bay

2 有關(guān)參數(shù)的數(shù)值試驗(yàn)結(jié)果

作者用K代表海底摩擦系數(shù),E代表水平湍流摩擦系數(shù),C代表海底粗糙度系數(shù)。對(duì)于這3個(gè)參數(shù)分別進(jìn)行試驗(yàn)。

對(duì)于K值, 以 1.25×10–3為間隔從 1×10–3到6×10–3取值。對(duì)于E值, 具體取值分別為：10, 1×102,5×102, 1×103, 5×103, 1×104, 5×104m2/s。表1 和表2是分別對(duì)應(yīng)不同K值和E值的計(jì)算方案, 北部灣中幾個(gè)計(jì)算點(diǎn) K1分潮潮振幅和位相的變化, 點(diǎn)的位置見圖1。結(jié)果分析發(fā)現(xiàn), 當(dāng)K取2×10-3～3×10-3,E取1×102～5×103m2/s時(shí)計(jì)算結(jié)果與實(shí)測結(jié)果比較符合,這與文獻(xiàn)[3]基本一致。

表1 不同K值K1潮汐要素的變化Tab. 1 Variation of the K1 tide at different values of K

表2 不同E值K1潮汐要素的變化Tab. 2 Variation of the K1 tide at different values of E

對(duì)于C值, 以 5×10–4為間隔從 0 m 到 6×10–3m取值,K取為 2×10–3,E取為 1×102m2/s。分析發(fā)現(xiàn), 當(dāng)C取 1×10–3～2×10–3m 時(shí), 計(jì)算結(jié)果與實(shí)測結(jié)果符合較好。當(dāng)C為1×10–3m時(shí)計(jì)算和實(shí)測潮汐調(diào)和常數(shù)的比較結(jié)果見表3, 其中ΔH和Δg分別代表計(jì)算和實(shí)測分潮的振幅和遲角的差值, 并由此計(jì)算出K1和M2分潮振幅和遲角的絕對(duì)平均誤差, 其中K1分潮振幅的絕對(duì)平均誤差為 5.4 cm, 遲角的絕對(duì)平均誤差為5.4°, M2分潮振幅的絕對(duì)平均誤差為3.5 cm, 遲角的絕對(duì)平均誤差為7.5°。

3 潮汐潮流的基本特征

3.1 潮汐

圖3給出的是K1分潮計(jì)算所得同潮圖, 其中實(shí)線表示遲角, 虛線表示振幅。由圖3清晰可見, 在北部灣存在一個(gè)退化了的無潮點(diǎn), 位于 107°5′E,16°40′N。在潮波由南向北的傳播過程中, 振幅不斷增加, 至北部灣灣頂振幅增至80 cm左右。

圖4給出的是M2分潮計(jì)算所得同潮圖, 其中實(shí)線表示遲角, 虛線表示振幅。由圖4可見, 北部灣M2分潮的振幅要比全日分潮小得多, 但作為前進(jìn)波,M2分潮明顯地由灣口沿灣軸向?yàn)稠攤鞑? 遲角不斷增大。在北部灣西北部的沿岸海域存在一個(gè)退化了的無潮點(diǎn), 這與孫洪亮等[2]計(jì)算結(jié)果基本一致。

表3 潮汐調(diào)和常數(shù)計(jì)算值與實(shí)測值比較Tab. 3 Comparison between the simulated and observed constituent constants

圖3 K1 分潮同潮圖Fig. 3 Cotidal chart of K1 constituent

圖4 M2 分潮同潮圖Fig. 4 Cotidal chart of M2 constituent

3.2 潮流

同北部灣的潮汐現(xiàn)象一樣, K1分潮潮流占主導(dǎo)地位, K1分潮的最強(qiáng)流速區(qū)出現(xiàn)在瓊州海峽與海南島西側(cè), 最大達(dá)70 cm/s。M2分潮的最強(qiáng)流速區(qū)則出現(xiàn)在瓊州海峽與海南島西北側(cè), 最大達(dá)20 cm/s左右。無論 K1還是 M2分潮, 它們的潮流長軸方向都與灣軸方向基本一致, 在瓊州海峽則呈東西向分布。圖5、圖6分別列出了北部灣 K1與M2分潮流表層的橢圓長短軸分布, 其中箭頭表示流速大小。從圖5, 6上可見, 在北部灣大部分海區(qū)無論 K1分潮還是M2分潮, 潮流基本呈往復(fù)式。除在瓊州海峽西側(cè), K1和M2分潮流表現(xiàn)出很強(qiáng)的旋轉(zhuǎn)性外, K1分潮流在海南島的西側(cè)也表現(xiàn)出較強(qiáng)的旋轉(zhuǎn)性。

圖5 K1分潮表層潮流橢圓分布Fig. 5 Distribution of K1 tidel current ellipses in the surface layer

圖6 M2分潮表層潮流橢圓分布Fig. 6 Distribution of M2 tidel current ellipses in the surface layer

3.3 潮余流

由于在北部灣中, K1分潮流占主導(dǎo)地位, 所以作者考慮K1分潮余流。圖7給出了K1分潮余流的分布,由圖7可見, 瓊州海峽的潮余流由東向西, 進(jìn)入北部灣后向西北而后向西南, 瓊州海峽的潮余流最大,約為15 cm/s; 海南島西岸有一由南向北的沿岸潮余流, 在 19°N附近與折向西南的潮余流匯合, 流速約10 cm/s。這也與孫洪亮等[2]計(jì)算結(jié)果基本一致。

圖7 表層潮余流分布Fig. 7 Distribution of tidel residual current in the surface layer

3.4 水平速度的垂向分布

在北部灣中選取H, G兩個(gè)點(diǎn)(點(diǎn)位置見圖1), 水深分別為60 m和24.9 m。圖8給出了H, G點(diǎn)水平速度的垂向分布特征, 垂向采用σ坐標(biāo), 水深范圍規(guī)范化為 0～1, 其中最大流速表層至底層基本一致,反映了潮流是正壓流動(dòng)的特點(diǎn); 最大流速方向隨深度稍微左偏; 海底附近流速變小, 反映了底摩擦的影響。最大流速發(fā)生時(shí)刻從表層到底層H, G兩點(diǎn)分別提前0.1 h和0.2 h。

圖8 水平速度垂向分布Fig. 8 Vertical distribution of the horizontal velocity

3.5 湍流黏滯系數(shù)的垂直分布

由湍能封閉模型得到的垂直黏滯系數(shù)同時(shí)具有時(shí)間和空間的變化, 圖9給出了 H點(diǎn)湍流黏滯系數(shù)的垂直分布在1～6 h內(nèi)的變化, 可以看出垂直分布為海底、海面小、中間大的拋物型, 最大值約在深度的中間偏下取得。

4 結(jié)論

利用 ECOMSED三維水動(dòng)力模式, 采用正交曲線坐標(biāo), 較好地?cái)M合地形, 對(duì)北部灣的潮汐、潮流和潮余流等特征進(jìn)行了分析, 可以得到以下結(jié)論：(1)當(dāng)海底摩擦系數(shù)、水平湍流摩擦系數(shù)、海底粗糙度系數(shù)分別取 2×10–3, 1×102m2/s, 1×10–3m 時(shí), 計(jì)算結(jié)果與驗(yàn)潮站觀測擬合較好, 其中 K1和 M2兩個(gè)分潮的絕對(duì)平均誤差：振幅分別為5.4 cm, 3.5 cm; 遲角分別為 5.4°, 7.5°。(2)在 107°5′E, 16°40′N 有一個(gè)退化了的 K1分潮無潮點(diǎn), 在西北部沿岸海域存在一個(gè)退化了的M2分潮無潮點(diǎn)。(3)瓊州海峽的K1分潮余流由東向西, 進(jìn)入北部灣后向西北而后向西南; 海南島西岸有一由南向北的沿岸潮余流。(4)最大流速表層至底層基本一致, 海底附近流速變小, 反映了底摩擦的影響, 最大流速的方向自上而下稍向左偏,最大流速發(fā)生時(shí)刻從表層到底層H, G兩點(diǎn)分別提前0.1 h和0.2 h。(5)湍流黏滯系數(shù)的垂直分布為海底、海面小、中間大的拋物型, 最大值約在深度的中間偏下取得。

圖9 H點(diǎn)垂直黏滯系數(shù)在6 h內(nèi)的變化(每隔1 h)Fig. 9 Variation of the vertical turbulent viscosity in 6 hours at point H (every 1 h)

[1]夏華永, 殷忠斌, 郭芝蘭, 等. 北部灣三維潮流數(shù)值模擬[J]. 海洋學(xué)報(bào), 1997, 19(2)：21-31.

[2]孫洪亮, 黃衛(wèi)民. 北部灣潮汐潮流的三維數(shù)值模擬[J]. 海洋學(xué)報(bào), 2001, 23(2)：1-8.

[3]殷忠斌, 陳明劍, 李樹華, 等. 北部灣潮汐數(shù)值計(jì)算參數(shù)的試驗(yàn)[J]. 廣西科學(xué), 1996, 3(2)：71-74.

[4]Blumberg A F. A primer for Ecomsed[M]. America：Hydroqual, Inc, 2002. 1-188.

[5]Chen Changsheng, Beardsley R, Franks P J S. A 3-D prognostic numerical model study of the Georges Bank ecosystem. Part I：physical model[J]. Deep Sea Research Part II: Topical Studies in Oceanography,2001, 48(1-3)：419-456.

[6]Ezer T, Mellor G L. Sensitivity studies with the North Atlantic sigma coordinate Princeton Ocean Model[J].Dynamics of Atmospheres and Oceans, 2000, 32(3-4)：185-208.

Numerical simulation of the tide and tidal current in Beibu Bay

XU Zhen-hua1, LEI Fang-hui2, LOU An-gang3,4, CAO Sheng-shan1
(1. Department of Mathematics, Ocean University of China, Qingdao 266100, China; 2. China Offshore Oil Research Center, Beijing 100027, China; 3. College of Environmental Science and Engineering, Ocean University of China, Qingdao 266100, China; 4. Key Laboratory of Marine Environmental Science and Ecology, Ocean University of China, Qingdao 266003, China)

Mar. , 6, 2006

Beibu Bay; tide; tidal current; numerical simulation

The model ECOMSED with the orthogonal curvilinear grid system is applied to simulating the tide and tidal current in Beibu Bay. We use different coefficients such as bottom friction, bottom roughness and horizontal turbulent friction for the numerical test. The results show that, when coefficient of bottom friction is 0.002 to 0.003,coefficient of bottom roughness is 1×10–3to 2×10–3m and coefficient of horizontal turbulent friction is 1×102to 5×103m2/s, the simulations of the tide and tidal currents agree generally with the observation. The characteristics of the tide, tidal current, and residual current in Beibu Bay are investigated based on the simulation. The vertical turbulent coefficients were found to display a parabolic structure. This study offers more understanding about the hydrodynamic condition in Beibu Bay and provides a scientific basis for the program of environmental protection.

P731.2 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1000-3096(2010)02-0010-05

2006-03-06;

2008-09-09

國家自然科學(xué)基金資助項(xiàng)目(40376034)

徐振華(1980-), 男, 山東青島人, 碩士, 主要從事海洋環(huán)境數(shù)值模擬研究, 電話：0532-87452228, E-mail：xuzhh@2008.sina.com

劉珊珊)