膠州灣潮汐潮流高分辨率數(shù)值模擬研究*

陳金瑞,陳學恩,于華明,閆運偉,單士亮,趙 健

(中國海洋大學海洋環(huán)境學院,山東青島266100)

近年來,隨著數(shù)值模式的完善和超級計算能力的形成,使對區(qū)域海洋進行高分辨率數(shù)值模擬研究成為可能。對膠州灣及近海的潮汐、潮流變化開展高分辨率數(shù)值模擬研究,不僅可以認識膠州灣及近海自身環(huán)流結構的特征與變化,膠州灣與黃海的相互影響,灣內污染物的擴散輸運,泥沙的輸運和沉積等科學問題;更有著指導港口航運,港口建設,水產養(yǎng)殖,近岸海上運動等實用價值。

膠州灣位于黃海之濱,山東半島的南岸,以團島頭(36°02′36″N,120°16′49″E)與薛家島腳子石(36°00′53″N,120°17′30″E)連線為界,是一個與黃海相通的半封閉式海灣。海灣東西寬27.8 km,南北長33.3 km;灣口開向東南,口門最窄處為3.1 km,岸線長187 km,平均水深約7 m,最大水深64 m。迄今,前人對膠州灣潮汐潮流的數(shù)值模擬研究做出卓有成效的工作[1-9],例如,閆菊等[4]曾模擬了膠州灣M 2分潮的潮汐潮流變化,給出了潮致余流的水平結構;孫英蘭等[5]基于ECOM模式,引入變邊界處理技術,指出潮流的量值及最大流發(fā)生時刻隨深度增加而減小,在底層附近,減小比率較大,主要是因為底摩擦的作用;高大魯?shù)萚8]利用POM模式對膠州灣進行的多個分潮的潮位和潮流模擬,討論了膠州灣的潮位、水平流速和垂直流速的性質;呂新剛等[9]從能量的角度探討了潮波的傳播特性,并且討論了漫灘過程對潮流模擬的重要性;趙亮等[7]利用ECOM對膠州灣潮波系統(tǒng)及其驅動下的標識質點運移規(guī)律進行數(shù)值模擬,指出膠州灣平均水存留時間為80 d,不同的區(qū)域水交換能力相差較大。然而,由于膠州灣地形和岸線條件復雜,受當時計算機性能的限制,上述數(shù)值研究采用的模式網(wǎng)格較粗,對岸線和底形的分辨率不夠高,故其所提供的該區(qū)的潮汐潮流分布信息不夠細致。

本文應用Finite-Volume Coastal Ocean Model(FVCOM)模式[10],采用2008年的最新岸線地形,設計高分辨率的三角形網(wǎng)格,考慮漫灘效應,用以研究膠州灣潮汐潮流的精細化結構及其時空分布規(guī)律。

1 模型配置和驗證

FVCOM陸架海洋模式的特點是,在水平方向上采用非結構化三角形網(wǎng)格,垂向采用σ坐標變換,數(shù)值方法采用有限體積方法,通過積分三角形控制體通量來解控制方程,將有限元擬合曲折岸線的優(yōu)點與有限差分方法高效計算的優(yōu)點結合起來,并且對于近岸、河口等具有復雜地形和岸線的區(qū)域來說能更好地保證質量、動量的守恒性。另外,FVCOM還包含三維干/濕網(wǎng)格處理等模塊。本文采用FVCOM對膠州灣進行數(shù)值模擬研究。對于像膠州灣這樣1個具有海底地形復雜,灘涂面積大,岸線彎曲不規(guī)則,潮差較大的海灣,選擇FVCOM模式進行水動力模擬是非常適合的。

1.1 模型配置

為從數(shù)值上再現(xiàn)膠州灣近岸海域潮汐潮流,本文研究區(qū)域設定為120.00°E～121.20°E,35.50°N～36.50°N。模式的水深資料來自中華人民共和國海事局出版的海圖,灣內的水深主要取自2008年的海圖數(shù)據(jù)。結合FVCOM的特點,本文對膠州灣口岸線變化復雜及水深坡度較大的區(qū)域和一些較小的島嶼進行了網(wǎng)格加密(見圖1)。模型水平采用無結構三角形網(wǎng)格系統(tǒng),垂向分7個σ層,計算區(qū)域內共有25 956個網(wǎng)格節(jié)點(node),49 729個三角形單元(cell),最小網(wǎng)格步長為50 m,模型外模積分步長為0.8 s,內模步長為4 s。從圖1中可見,三角形網(wǎng)格較好地概括了計算域內復雜的島嶼岸線和地形特征,對復雜岸線的擬合較為精確。對于水平湍流混合系數(shù),FVCOM中采用了Smago rinsky渦動參數(shù)化方案。

式中:C為常數(shù),Ω為每個三角形網(wǎng)格的面積。

圖1 膠州灣模式的計算范圍和計算網(wǎng)格Fig.1 Calculation range and mesh in Jiaozhou Bay

1.2 邊界條件

本文主要進行數(shù)值診斷研究,不考慮溫鹽變化,溫鹽設為常數(shù),32 PSU和20℃;也不考慮大氣對海洋的影響。模式采用零初始條件(t=0時,U=V=η=0),潮位和流速均設為0。開邊界的73個節(jié)點采用8個主要天文分潮(M 2、S2、N2、K2、K1、O1、P1、Q1)的調和常數(shù),然后利用t_tide軟件后報出開邊界的水位,來驅動本文的數(shù)值模式。調和常數(shù)采用OTIS(OSU Tidal Inversion Software)數(shù)據(jù),并融合膠州灣近海的歷史同潮圖數(shù)據(jù)。數(shù)值積分時間跨度為35 d,2008年7月25日～2008年8月31日。

圖2 膠州灣海域的水深(m)和水位潮流驗證點分布Fig.2 Distribution of water dep th and observation stations in Jiaozhou Bay

1.3 模型驗證

本文選取青島五號碼頭和小麥島驗潮站2008年8月1日00··00～2008年8月30日23··00的逐時觀測資料(見圖3),來驗證數(shù)值結果,水位對比曲線表明二者吻合較好。五號碼頭和小麥島2個站位的模擬值與觀測值后報的潮位的平均絕對誤差均在10%以內。表1給出了2個潮位驗證點的振幅和遲角對比,除了S2和K1的遲角偏差稍大,M 2、N 2、O1分潮模擬與實測結果吻合較好,模擬精度可以接受。

圖3 五號碼頭和小麥島潮位模擬與實測值比較Fig.3 Comparison between the observed and computed tidal elevation at Xiaomaidao and Wuhaomatou

表1 觀測與計算的調和常數(shù)對比表Table 1 Comparison between the observed and computed tidal constituent

綜上所述,通過對流速流向對比表明,本文所建立的膠州灣潮流模式是可信的,模擬結果能較好地體現(xiàn)膠州灣海域的潮汐潮流特征。

2 結果分析

2.1 潮汐分析

圖4給出了膠州灣的潮汐性質判別數(shù),膠州灣及灣口近海的潮汐性質判別數(shù)介于0.35～0.40之間,屬規(guī)則半日潮,該海區(qū)半日潮占優(yōu)。對模擬結果的調和分析表明(見圖5),膠州灣M 2、S2、K1、O1分潮波都主要從外海傳入,均不構成完整獨立的潮汐系統(tǒng)。在膠州灣外,分潮波自東北向西南沿逆時針方向傳播,振幅由外海向灣內增大。灣口與灣頂M 2分潮的振幅差約10 cm,4個分潮的潮汐最大振幅分別為1.34、0.51、0.262、0.208 m。

圖4 膠州灣海域潮汐性質Fig.4 Distribution of calculated tidal characteristics in Jiaozhou Bay

圖5 膠州灣海區(qū)分潮同潮圖Fig.5 Sketch map of co-tidal lines in Jiaozhou Bay

2.2 潮流分析

對數(shù)值結果計算表明,膠州灣海域潮流類型判別數(shù)介于0.2～0.3之間,也具有規(guī)則的半日潮流的性質。圖6給出了小麥島水文站大潮時段的漲落急時刻潮流分布圖,外海的海水偏西向進入膠州灣外灣口后分叉:一支偏西南向進入黃島前灣和海西灣;主流繞過團島咀,偏西北向進入膠州灣內灣口再分叉:一支北偏東向進入滄口水道;一支北偏西向,由中沙礁西側進入內灣;余部由中央水道北進[14],漲、落潮最大流速均發(fā)生在膠州灣灣口。

圖7給出了研究區(qū)域M 2、K1分潮的表層潮流橢圓水平分布(紅線:順時針旋轉;黑線:逆時針旋轉)。從圖7a可以看出,在灣外,M 2和K1分潮的潮流橢圓旋轉率較大,以旋轉流為主,潮流橢圓主要呈逆時針方向旋轉;在灣內則主要呈順時針方向旋轉;由于灣口附近岸線復雜和水深坡度大,潮流以往復流為主,潮流橢圓的旋轉方向變化復雜。灣口及灣內,潮流橢圓旋轉率較小且為負值,以往復流為主。各個分潮的橢圓長軸最大值都出現(xiàn)在灣口附近的水道上,M 2、S2、K1和O1的最大橢圓長軸分別為1.13,0.41,0.126和0.12 m/s。

圖6 大潮漲落急時刻的潮流分布Fig.6 Distribution of tidal current

圖7 潮流橢圓分布Fig.7 Distributions of tidal current ellipse

2.3 最大可能流速

最大可能流速由對模擬結果進行調和分析所得的調和常數(shù)來計算,對于規(guī)則半日潮海區(qū),最大可能流速的算式為[11]:

其中WM2、WS2、WK1、WO1分別為M 2、S2、K1、O1分潮潮波的潮流橢圓長半軸。根據(jù)上述公式,計算膠州灣海域的最大可能流速,繪制最大可能流速分布于圖8,從圖8可以看出,最大可能流速分布基本與等深線一致,與岸線平行。在膠州灣內最大可能流速介于0.5～0.9 m/s之間,灣口處因應岸線和地形變化,最大可能流速的變化梯度較大,最大可能流速可達2.14 m/s。

圖8 最大可能流速分布圖Fig.8 Distributions of maximum probable current velocities

2.4 潮致余流

潮致歐拉余流速度是指海域內某一確定點在1個潮周期內潮流速度的時間平均值,它表示在確定位置上流體周期平均的遷移趨勢,因在淺海中潮流的非線性項受到側向岸線和底摩擦的作用,一部分周期性能量會轉變成非周期性能量。余流的大小與地形和岸線有密切的聯(lián)系。根據(jù)模式的結果計算潮致余流(見圖9),在灣口處,由于岬角地形和岸線的共同作用,水平結構顯著地表達為“團團轉”的多渦結構[4],在團島咀處的岬角鋒,形成4個主要渦旋環(huán)流系統(tǒng)(A、B、C、D),且潮汐余流的最大流速在團島附近可達0.43 m/s,與呂新剛等[9]根據(jù)實測資料推算所得團島嘴西側近岸水域最大潮余流0.5 m/s吻合。除了這4個強環(huán)流系統(tǒng)[4、9],本研究在膠州灣中部,由于模式精度的提高,還首次揭示了2個相對較弱的余環(huán)流系統(tǒng)(E和F),流速小于0.1 m/s。

圖9 膠州灣的余流場Fig.9 Distribution of residual currents in Jiaozhou bay

2.5 潮流能通量

為反映膠州灣及其近海潮汐能的分布和傳播,本文通過潮流能通量來估算潮流能[12-13]:

P潮流能通量,ρ密度,g重力加速度,U體積輸運矢量,v流速,H水深,ζ水位,<>代表時間平均。潮流能通量不僅可以反映潮流能量的分布,而且可以反映潮流能量的傳播。結果表明(見圖10),膠州灣外海海區(qū)的潮流能通量主體自東北向西南方向傳播,基本上沿岸線;在接近灣口時,一部分潮流能沿著海岸線繼續(xù)向西南傳播,另一部分潮流能向灣口匯集并繼續(xù)向灣內傳播,經(jīng)過灣口,潮流能由匯聚轉為發(fā)散,傳向內灣。灣口附近由于地形和岸線的共同作用,潮流能通量在內外灣口分別做順時針旋轉,呈“左進右出”的結構。潮流能通量在內外灣口“左進”的量值比較大,“右出”的量值比較小,灣口附近存在1個潮能通量的高值區(qū),這是因為,灣口水深梯度較大和岸線共同導致流速較大,進而潮流能通量也較大。內灣大部分區(qū)域潮能通量比較弱。

圖10 膠州灣的潮流能通量的水平分布Fig.10 Horizontal distribution of tidal energy flux in Jiaozhou Bay

2.6 海灣納水體積與納潮量的計算

納潮量是表征半封閉海灣生命力的重要指標,其大小直接影響到海灣與外海的海水交換強度和污染物的遷移擴散,從而制約著海灣的自凈能力和環(huán)境容量,對于維護海灣良好的生態(tài)環(huán)境至關重要,納潮量的變化直接影響到海灣的潮流特性,還可能打破水動力條件與海灣形態(tài)之間的動態(tài)平衡,使海灣槽灘形態(tài)隨之進行調整。海灣面積的減小會造成納潮量的減小。

根據(jù)模式的計算結果,分析納潮量列于表2。結果表明,團島-黃島油碼頭以內(內灣)的海域面積為3.09億m2,平均納潮量為7.36億m3;團島-薛家島以內(外灣加灣口)的海域面積為3.487億m3,平均納潮量為8.31億m3。

2.7 Dye濃度的擴散

為分析內灣與外灣的水交換,本文設計了數(shù)值染色實驗來直觀地反應水體交換情況。

在數(shù)值模式中采用下述濃度擴散方程:

其中,C染色實驗的濃度,D總水深,u、v和ω為流速的分量,Kh垂直熱擴散系數(shù),Fc水平擴散項,C0為初始濃度或者點源的濃度。

取團島、薛家島之間的連線作為內、外灣分界線,將內灣初始污染物濃度設置為1,外灣初始污染物濃度設為0,模式其它設置與水動力配置完全相同,進行污染物擴散數(shù)值模擬實驗。圖12b給出了數(shù)值積分30 d后內灣污染物濃度的分布:30 d后,內灣余留污染物的平均濃度為63.2%,即30 d的水交換率為36.8%,可估算灣內整體水交換約需要82 d(30÷36.8%≈82),這與趙亮等[7]膠州灣平均水存留時間80 d這一結論吻合。

表2 大小潮膠州灣水體體積和納潮量Table 2 The volume of water and tidal prism of Jiaozhou Bay at spring and neap tidal/億m3

圖11 污染物濃度分布Fig.11 Distributions of dye

3 結論

本文基于三維變邊界水動力模型FVCOM,考慮8個主要天文分潮,對膠州灣及其近海海域的潮汐、潮流進行了高精度的數(shù)值模擬研究,三角形網(wǎng)格很好地擬合了不規(guī)則岸線,高精度地刻畫了較小尺度的島嶼和復雜岸線,提高了模式的模擬精度,干濕網(wǎng)格技術處理漫灘問題,使潮汐潮流的模擬結果更加可靠。

在結合歷史觀測資料和前人工作,對模式進行對比驗證的基礎上,揭示了膠州的潮波主要特征:

(1)膠州灣的潮汐類型為屬于規(guī)則半日潮。灣外分潮的潮流橢圓旋轉率較大,以旋轉流為主,逆時針方向旋轉;在灣口及灣內,潮流橢圓旋轉率較小且為負值,往復流為主,各個分潮的橢圓長軸最大值都出現(xiàn)在灣口附近的水道上,M 2、S2、K1和O1的最大橢圓長軸分別為1.13、0.41、0.126和0.12 m/s。

(2)高精度的數(shù)值模擬揭示了6個較大的潮汐余環(huán)流系統(tǒng),余流最大值達0.43 m/s;本研究還在膠州灣中部首次揭示了2個相對較弱的余環(huán)流系統(tǒng)(E和F),流速小于0.1 m/s。膠州灣的最大可能流速的分布基本與等深線一致,最大值可達2.14 m/s。

(3)灣口附近由于地形和岸線的共同作用,潮流能通量在內外灣口分別做順時針旋轉,呈“左進右出”的結構,從外灣口進入的潮波能量有41.7%耗散在內外灣口之間的海區(qū)。

(4)納潮量研究結果表明,膠州灣的平均納潮量為8.31億m3,染色試驗表明,膠州灣30天的水交換率為36.8%。

致謝:在本文研究過程中,作者們還得到中國海洋大學海洋環(huán)境學院計算中心的支持,在此表示感謝。

[1] Chen Changsheng,Ji Rubao,Zheng Lianyuan.Influencesof physical p rocesseson the ecosystem in Jiaozhou Bay:A coupled physical and biological model experiment[J].J Geophys Res,1999,104:29925-29949.

[2] 王化桐,方欣華,匡國瑞,等.膠州灣環(huán)流和污染擴散的數(shù)值模擬Ⅰ膠州灣潮流數(shù)值計算[J].中國海洋大學學報:自然科學版,1980,10(01):26-63.

[3] 鮑獻文,閆菊,趙亮,等.ECOM模式在膠州灣潮流計算中的應用[J].海洋科學,1999,(05):56-60.

[4] 閆菊,王海,鮑獻文.膠州灣三維潮流及潮致余環(huán)流的數(shù)值模擬[J].地球科學進展,2001,16(02):172-177.

[5] 孫英蘭,張越美.膠州灣三維變動邊界的潮流數(shù)值模擬[J].海洋與湖沼,2001,32(04):355-362.

[6] 孫長青,王學昌,孫英蘭.填海造地對膠州灣污染物輸運影響的數(shù)值計算[J].海洋科學,2002,26(10):47-50.

[7] 趙亮,魏皓,趙建中.膠州灣水交換的數(shù)值研究[J].海洋與湖沼,2002,33(01):23-29.

[8] 高大魯,魏澤勛,華鋒.膠州灣多分潮漫灘數(shù)值模擬研究[J].海洋科學進展,2007,25(02):131-138.

[9] 呂新剛,喬方利,夏長水.膠州灣潮汐潮流動邊界數(shù)值模擬[J].海洋學報,2008,30(04):21-29.

[10] Chen Changshen,Huang Haosheng,Beardsley Robert C,et al.A finite volume numerical approach for coastal circulation studies:Comparisons with finite difference models[J].Journal of Geophysical Research,2007,112(C03018):1-34.

[11] 宋德海,鮑獻文,朱學明.基于FVCOM的欽州灣三維潮流數(shù)值模擬[J].熱帶海洋學報,2009,28(02):7-14.

[12] Garrett C.Tides in gulfs[J].Deep-Sea Research,1975,22:23-35.

[13] Zu Tingting,Gan Jianping,Erofeeva Svetlana Y.Numerical study of the tide and tidal dynamics in the South China Sea[J].Deep-Sea Research,2008,55:137-154.

[14] 中國海灣志編纂委員會.中國海灣志第四分冊(山東半島南部和江蘇省海灣)[M].北京:海洋出版社,1993:448.