陳金瑞,陳學(xué)恩,于華明,閆運(yùn)偉,單士亮,趙 健

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

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

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

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

1 模型配置和驗(yàn)證

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

1.1 模型配置

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

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

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

1.2 邊界條件

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

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

1.3 模型驗(yàn)證

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

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

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

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

2 結(jié)果分析

2.1 潮汐分析

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

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

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

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

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

2.3 最大可能流速

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

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

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

2.4 潮致余流

潮致歐拉余流速度是指海域內(nèi)某一確定點(diǎn)在1個(gè)潮周期內(nèi)潮流速度的時(shí)間平均值,它表示在確定位置上流體周期平均的遷移趨勢,因在淺海中潮流的非線性項(xiàng)受到側(cè)向岸線和底摩擦的作用,一部分周期性能量會轉(zhuǎn)變成非周期性能量。余流的大小與地形和岸線有密切的聯(lián)系。根據(jù)模式的結(jié)果計(jì)算潮致余流(見圖9),在灣口處,由于岬角地形和岸線的共同作用,水平結(jié)構(gòu)顯著地表達(dá)為“團(tuán)團(tuán)轉(zhuǎn)”的多渦結(jié)構(gòu)[4],在團(tuán)島咀處的岬角鋒,形成4個(gè)主要渦旋環(huán)流系統(tǒng)(A、B、C、D),且潮汐余流的最大流速在團(tuán)島附近可達(dá)0.43 m/s,與呂新剛等[9]根據(jù)實(shí)測資料推算所得團(tuán)島嘴西側(cè)近岸水域最大潮余流0.5 m/s吻合。除了這4個(gè)強(qiáng)環(huán)流系統(tǒng)[4、9],本研究在膠州灣中部,由于模式精度的提高,還首次揭示了2個(gè)相對較弱的余環(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體積輸運(yùn)矢量,v流速,H水深,ζ水位,<>代表時(shí)間平均。潮流能通量不僅可以反映潮流能量的分布,而且可以反映潮流能量的傳播。結(jié)果表明(見圖10),膠州灣外海海區(qū)的潮流能通量主體自東北向西南方向傳播,基本上沿岸線;在接近灣口時(shí),一部分潮流能沿著海岸線繼續(xù)向西南傳播,另一部分潮流能向?yàn)晨趨R集并繼續(xù)向?yàn)硟?nèi)傳播,經(jīng)過灣口,潮流能由匯聚轉(zhuǎn)為發(fā)散,傳向內(nèi)灣。灣口附近由于地形和岸線的共同作用,潮流能通量在內(nèi)外灣口分別做順時(shí)針旋轉(zhuǎn),呈“左進(jìn)右出”的結(jié)構(gòu)。潮流能通量在內(nèi)外灣口“左進(jìn)”的量值比較大,“右出”的量值比較小,灣口附近存在1個(gè)潮能通量的高值區(qū),這是因?yàn)?灣口水深梯度較大和岸線共同導(dǎo)致流速較大,進(jìn)而潮流能通量也較大。內(nèi)灣大部分區(qū)域潮能通量比較弱。

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

2.6 海灣納水體積與納潮量的計(jì)算

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

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

2.7 Dye濃度的擴(kuò)散

為分析內(nèi)灣與外灣的水交換,本文設(shè)計(jì)了數(shù)值染色實(shí)驗(yàn)來直觀地反應(yīng)水體交換情況。

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

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

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

表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 結(jié)論

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

在結(jié)合歷史觀測資料和前人工作,對模式進(jìn)行對比驗(yàn)證的基礎(chǔ)上,揭示了膠州的潮波主要特征:

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

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

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

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

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

[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)流和污染擴(kuò)散的數(shù)值模擬Ⅰ膠州灣潮流數(shù)值計(jì)算[J].中國海洋大學(xué)學(xué)報(bào):自然科學(xué)版,1980,10(01):26-63.

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

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

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

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

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

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

[9] 呂新剛,喬方利,夏長水.膠州灣潮汐潮流動邊界數(shù)值模擬[J].海洋學(xué)報(bào),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] 宋德海,鮑獻(xiàn)文,朱學(xué)明.基于FVCOM的欽州灣三維潮流數(shù)值模擬[J].熱帶海洋學(xué)報(bào),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.