蘇 博朱金龍王佳瑩趙玉庭王立明馬元慶*

(1.山東省海洋資源與環(huán)境研究院,山東 煙臺264006;2.山東省海洋生態(tài)修復(fù)重點實驗室,山東 煙臺264006)

萊州灣位于山東半島西北、渤海南部,與遼東灣和渤海灣并稱為渤海三大海灣。萊州灣的海灣面積占據(jù)山東各海灣面積之首,其灣口西起黃河新入?？?東至屺姆島高角,口門寬約96 km[1]。受渤海半封閉型內(nèi)海地理區(qū)位特點的影響,萊州灣內(nèi)海水交換能力差,沿岸排入污染物總量遠(yuǎn)超海域海水的自凈能力,導(dǎo)致水體質(zhì)量不斷惡化[2-3]。近年來,有學(xué)者和社會輿論提出修建膠萊運河[4](膠萊人工海河),運河可使膠州灣海水輸入萊州灣,增強萊州灣水交換,同時緩解萊州灣乃至渤海的海洋生態(tài)環(huán)境污染現(xiàn)狀[5]。然而,膠萊運河是在膠萊河基礎(chǔ)上從膠州灣到萊州灣開挖航道貫穿山東半島的,該工程體量和修建耗資都很龐大,須對其修建可行性和效果做充分、科學(xué)的研究。此前相關(guān)學(xué)者做過初步論證,如左軍成等[6]分析了兩灣的潮流和環(huán)流結(jié)構(gòu)以及膠萊運河對膠州灣和萊州灣之間的水交換的貢獻,認(rèn)為膠萊運河的開通將極大改善膠萊兩灣的水環(huán)境條件,不足的是其未采用數(shù)學(xué)模型進一步論證說明;王佳等[7]通過一、二維耦合模式模擬分析認(rèn)為工程的實施對渤海水質(zhì)影響效果不顯著,但其模型模擬網(wǎng)格分辨率較低。

以往研究成果未能明確修建膠萊運河對萊州灣海洋環(huán)境的具體影響,所以,本文利用渤海、黃海、萊州灣地形地形及海洋環(huán)境資料,選擇萊州灣作為重點研究區(qū)域,建立二維平面潮流和對流擴散模型,結(jié)合2018年萊州灣海水水質(zhì)監(jiān)測結(jié)果,研究膠萊運河修建對萊州灣水交換及海洋環(huán)境的影響,旨在為此重大項目的可行性及必要性進行初步的研究。

1 材料與方法

1.1 數(shù)據(jù)與資料來源

模型岸線提取自中國“高分一號”衛(wèi)星,影像地面分辨率為2 m,成像時間為2015年。應(yīng)用ArcGIS軟件對遙感影像進行校正和配準(zhǔn),經(jīng)校正配準(zhǔn)后的遙感影像依據(jù)多年平均大潮高潮時的海陸分界線所具備的光譜特征,采用人機交互解譯獲取岸線分布信息。模型水深選取中國人民解放軍海軍航海保證部10011號海圖的地形水深測量資料,該數(shù)據(jù)基于中國人民解放軍海軍海道測量局長期測量結(jié)果。以萊州灣水深為例,萊州灣水深自灣口向灣內(nèi)遞減,最大水深位于屺姆島高角附近,超過20 m,灣內(nèi)水深一般在10 m 以內(nèi)(圖1)。

圖1 模型區(qū)域的岸線與水深Fig.1 Coastline and water depth in the modeling area

1.2 數(shù)值模型

1.2.1 水動力控制方程

水動力模型的數(shù)值解基于二維淺水方程所得,原始方程的空間離散采用的是單元中心有限體積法[8],正交笛卡爾坐標(biāo)系。平面二維連續(xù)性方程和動量方程:

式中,h為總水深,η為水位,t為時間,分別為x和y向的水深平均流速,x和y為笛卡爾坐標(biāo),S為源強,f為科氏系數(shù),g為重力加速度,τsx和τsy分別為x和y向的風(fēng)剪切力,τbx和τby分別為x和y向的底床剪切力,ρ為水的密度,A為水平渦黏系數(shù)。

1.2.2 定解條件

數(shù)學(xué)模型計算域岸線、水深及范圍見圖1。膠萊運河橫跨膠東半島,連通黃、渤兩海,研究尺度較大,故本近海數(shù)值模型采用水位控制開邊界(AB連線),使得潮波傳播到研究區(qū)域時與實際情況相符。考慮到膠萊運河修建的經(jīng)濟性和施工難度,河道參數(shù)設(shè)置為河寬1 km,水深5 m。模型中海域采用非結(jié)構(gòu)三角形網(wǎng)格,運河河道采用四邊形網(wǎng)格。模型最大網(wǎng)格面積為100 000 m2,最小網(wǎng)格面積為500 m2。為了能夠精確地分析運河工程對周邊海域水動力環(huán)境的影響,模擬過程中,將工程附近海域進行網(wǎng)格加密,加密區(qū)域確定為整個萊州灣和膠州灣,工程附近計算網(wǎng)格見圖2。

初始條件:當(dāng)開始計算時,采用計算初始時刻相對應(yīng)的邊界條件。

邊界條件:將研究區(qū)域的陸地邊界作為固定邊界,假定垂直于固定邊界的法向速度為零。隨漲、落潮交替出現(xiàn)的邊灘,采用“干、濕”動邊界技術(shù)[9-10];開邊界水位根據(jù)MIKE Tidal工具箱提供的全球潮汐分潮數(shù)據(jù)(the global tide model data)推算,全球潮汐分潮數(shù)據(jù)來源于潮汐模型同化Topex/Poseidon衛(wèi)星的測高數(shù)據(jù)(altimetry data),是由全球10個主要天文分潮(日分潮5個,為S1,K1,O1,P1和Q1;半日分潮4個,為M2,S2,K2和N2;淺水分潮1個,為M4)[11-12]組成,分潮數(shù)據(jù)空間分辨率為0.125°×0.125°。

圖2 工程附近計算網(wǎng)格Fig.2 Computing grids in the areas near the project

1.3 模型驗證

底床阻力系數(shù)由曼寧阻力系數(shù)來確定。隨水深的增加,曼寧阻力系數(shù)取值為0.015～0.022[13]。采用Smagorinsky公式[8]計算水平渦流黏度,Smagorinsky系數(shù)Cs取值為0.24。風(fēng)場的時序資料來自美國國家海洋和大氣管理局(National Oceanic and Atmospheric Administration,NOAA)。風(fēng)力拖曳系數(shù)ca=1.255×10-3。

利用黃、渤海海域內(nèi)濰坊港和青島港等24個潮位站(圖3)的潮汐調(diào)和常數(shù)計算得到的潮位值與模型計算得到的潮位預(yù)測值進行對比驗證,限于篇幅,潮位驗證結(jié)果僅展示接近膠萊運河河口位置的濰坊港和青島港站(圖4)。驗證結(jié)果表明,所有潮位站潮汐相位誤差不超過0.25 h,潮位峰谷值誤差不超過0.10～0.15 m,這說明模型潮位能夠反映實際的潮汐特征。

利用大潮期萊州灣內(nèi)濰坊港、刁龍咀和屺姆島三個潮流觀測站以及膠州灣內(nèi)象頭東和青島港兩個潮流觀測站(圖3)與模型計算得到的潮流預(yù)測值進行對比驗證,限于篇幅,潮流驗證僅展示運河河口位置的濰坊港和象頭東站(圖5)。驗證結(jié)果分析,所有潮流觀測站位流速峰谷值的絕對誤差不大于漲落潮最大流速值絕對值之和的20%,流速相位差一般也不超過0.25 h,說明模型潮流能夠反映實際的潮流特征。

圖3 黃渤海海域內(nèi)24個潮位站的分布Fig.3 Distribution of 24 tide stations in the Yellow Sea and the Bohai Sea

圖4 濰坊港和青島港兩潮位站驗證結(jié)果Fig.4 Verifications of tide level at two stations:Weifang Port and Qingdao Port

圖5 濰坊港和象頭東兩潮位站流速、流向驗證結(jié)果Fig.5 Verifications of flow velocity and direction at two tide Stations:Weifang Port and Xiangtoudong

1.4 水交換計算

建立水交換模式采用的是保守物質(zhì)的對流-擴散方程:

式中,為保守物質(zhì)的垂向平均濃度,Dh為水平擴散系數(shù)。初始時刻,萊州灣內(nèi)水體保守物質(zhì)的相對濃度設(shè)為1,灣外水體的相對濃度為0。模擬海灣某一時刻平均濃度的計算公式為:

式中,n為海灣所占網(wǎng)格數(shù),C i為第i個網(wǎng)格的污染物濃度,H i為第i個網(wǎng)格的即時深度,A i為第i個網(wǎng)格的面積。灣內(nèi)某一時刻水交換率的計算公式為。根據(jù)半水交換[14]的概念,半交換周期為灣內(nèi)的平均濃度降為初始濃度一半所需的時間。

2 結(jié)果與討論

圖6 萊州灣和膠州灣河口處潮位過程線Fig.6 Tidal level process in the estuary areas of the Laizhou Bay and the Jiaozhou Bay

2.1 水動力計算結(jié)果與分析

由于缺乏觀測數(shù)據(jù),在此給出萊州灣和膠州灣膠萊運河河口處(圖2中1#點、5#點)大潮同時期的潮位過程模擬結(jié)果如圖6所示。由圖6可知,膠州灣河口為正規(guī)半日潮,萊州灣河口為不正規(guī)半日潮,膠州灣河口的潮差遠(yuǎn)大于萊州灣河口的潮差。兩灣河口處的潮位位相幾乎相反,即當(dāng)萊州灣河口漲潮時膠州灣河口落潮,當(dāng)萊州灣河口漲急時膠州灣河口落急,反之亦然。

由大潮期萊州灣河口漲急、膠州灣河口落急時刻潮流場(圖7a)可知:漲急時刻潮流由萊州灣外向灣內(nèi)運動,流速一般為0.10～0.78 m/s;最大流速出現(xiàn)在黃河口、萊州淺灘處,最小流速出現(xiàn)在灣頂附近;灣內(nèi)流速整體呈現(xiàn)西部流速略大于東部流速的空間特征。流向整體上呈現(xiàn)偏SW 向,在黃河口南部海域,潮流流向偏NW 向;在刁龍咀南部海域,潮流流向偏SE向。膠州灣河口落急時刻海水流向與漲急時刻相反,海水由灣內(nèi)經(jīng)灣口流向灣外,灣口流速大于灣內(nèi),最大流速出現(xiàn)在灣口附近,流速最大值約1.4 m/s,最小流速均出現(xiàn)在灣頂,流速約0.2 m/s。灣內(nèi)漲潮流向由灣頂指向灣口。

由大潮期萊州灣河口落急、膠州灣河口漲急時刻潮流場(圖7b)可知:大潮期間落急時刻萊州灣潮流分布情況與漲急時刻潮流場相似,只是流向與漲急時刻潮流場剛好相反,落急時刻潮流流速一般介于0.10～1.12 m/s,最大流速出現(xiàn)在黃河口、萊州淺灘處,最小流速出現(xiàn)在灣頂附近,灣內(nèi)流速整體呈現(xiàn)西部流速略大于東部流速的空間特征。流向整體上呈現(xiàn)偏NE向,在黃河口南部海域,潮流流向偏SE向;在刁龍咀南部海域,潮流流向偏NW 向。膠州灣河口漲急時刻海水經(jīng)灣口由灣外流入灣內(nèi),灣口流速大于灣內(nèi),最大流速出現(xiàn)在灣口附近,流速最大值約1.5 m/s,最小流速均出現(xiàn)在灣頂,流速約0.2 m/s。灣內(nèi)漲潮流向由灣口指向灣頂。

圖7 萊州灣和膠州灣的落急和漲急流場Fig.7 The flood and ebb flow fields in the Laizhou Bay and the Jiaozhou Bay

圖8 河道特征點流速與流向過程線Fig.8 Process of flow velocity and direction at character points in the river channel

為分析河道水動力環(huán)境特征,鑒于膠萊河沿線較長,選擇膠萊河兩端口門內(nèi)河段5個典型位置(圖2),繪制其48 h的流速與流向過程曲線圖,見圖8,其中1#點位于萊州灣河口、5#點位于膠州灣河口。膠州灣口門的5#點最大流速為1.1 m/s左右,其他各特征點最大流速不超過0.5 m/s。1#～5#點的流向每隔一天存在1 h左右由膠州灣至萊州灣的單向流動,并且每一轉(zhuǎn)流周期,1#點萊州灣入流歷時大于出流歷時約7 h。模型統(tǒng)計了膠州灣河口與萊州灣河口斷面海水進、出的累計流量。海水經(jīng)膠州灣河口斷面進入河道的累計流量與海水經(jīng)萊州灣河口斷面進入萊州灣的累計流量隨時間變化的量值一致,膠萊運河年徑流總量約為7.2×109m3,凈流方向由膠州灣流向萊州灣。

2.2 水交換能力影響分析

為研究膠萊運河工程實施對萊州灣納潮量的影響,在灣口斷面計算流量獲得整個斷面的流量時間序列,然后按漲潮期或落潮期進行流量累計可得海灣潮量。運河工程實施前后萊州灣海域大潮期和小潮期納潮量計算結(jié)果顯示,工程實施后萊州灣大小潮期納潮量均有所減少,其中大潮期納潮量減少3.2×107m3,小潮期納潮量減少1.7×107m3,大潮期納潮量減少值大于小潮期。

以位于119°24′E經(jīng)線附近、垂直于萊州灣灣口連線的直線為界將萊州灣分為西部和東部兩個部分(圖2),根據(jù)水交換率計算公式及半交換周期的概念表述,模擬得到工程前后萊州灣整體、西部、東部的海水水交換率隨時間變化與工程前后萊州灣半交換周期(圖9、表1),工程后整個萊州灣水交換能力有略微增強趨勢,全年平均與現(xiàn)狀相比增強約2.0%,工程后半交換周期(113 d)比現(xiàn)狀(127 d)縮短14 d;西部萊州灣水交換能力基本沒有變化,略微增強約0.5%;東部萊州灣受膠萊河入水口的影響,水交換比西部強,增強約2.7%。由圖10可見,萊州灣水體整體上是自西向東交換,灣內(nèi)運河處出水口交換率最高。運河工程實施后,濰坊港防波堤至刁龍咀海域水交換能力略有增強,濰坊港防波堤以西海域水交換能力基本未發(fā)生變化,即運河工程的實施對萊州灣水交換能力影響范圍主要位于濰坊港防波堤至刁龍咀之間海域。

圖9 水交換率隨時間變化Fig.9 Changes of water exchange rate with time

表1 萊州灣水交換率變化Table 1 Variations of the water exchange rate in the Laizhou Bay %

續(xù)表 %

圖10 工程前后在160 d萊州灣水交換物質(zhì)濃度分布(初始值為1)Fig.10 Concentration distribution of water exchange substances in 160 days before and after the project(initial value:1)

2.3 萊州灣海水質(zhì)量狀況的影響分析

2018年山東省海洋生態(tài)環(huán)境公報①山東省生態(tài)環(huán)境廳.山東省海洋生態(tài)環(huán)境狀況公報.2018.中指出,山東省主要的水質(zhì)污染(低于二類水質(zhì)標(biāo)準(zhǔn))區(qū)域在渤海沿岸,其中萊州灣的海水水質(zhì)尤其差,相對應(yīng)的,膠州灣內(nèi)雖有部分區(qū)域水質(zhì)較差,但總體相對于萊州灣水質(zhì)要好(圖11)。從海水水質(zhì)分布趨勢來看,萊州灣東部海水質(zhì)量普遍高于西部,萊州灣西部陸源污染強度要比東部大,主要是因為黃河及小清河入海徑流大,攜帶陸源污染較多,另外萊州灣西部水交換效率也明顯低于東部。從膠萊河工程前后160 d水交換物質(zhì)濃度差異(圖11a)來看,工程前后萊州灣入流口水交換物質(zhì)濃度差最大為0.5(初始為1),但是影響范圍較小,僅僅影響河口向海5 km 的扇面海域。通過數(shù)模計算出年輸送海水量,約為7.2×109m3(南膠萊河年入海徑流量為2.49×108m3,北膠萊河年入海徑流量為4.96×108m3),黃海海水經(jīng)膠州灣由膠萊運河輸送至萊州灣,沖淡萊州灣內(nèi)非優(yōu)良水體,可減少水質(zhì)污染面積約1.5 km2,在一定程度上可緩解萊州灣東部河流入?？诟浇》秶Ｓ蛑械暮Ｋ|(zhì)污染現(xiàn)狀。

圖11 2018年萊州灣及膠州灣海域水質(zhì)等級及工程前、后160 d水交換物質(zhì)濃度差異Fig.11 Seawater quality grades in the Laizhou Bay and the Jiaozhou Bay in 2018 and difference in concentration of water exchange substances in 160 days before and after the project

3 結(jié) 論

建立了二維潮流和對流擴散模型,結(jié)合海水水質(zhì)現(xiàn)狀,初步模擬分析了膠萊運河的貫通對萊州灣水動力條件及海洋環(huán)境現(xiàn)狀的影響,得到結(jié)論:膠萊運河開通后年徑流總量約為7.2×109m3,凈流方向為膠州灣流向萊州灣。萊州灣水體交換整體上是自西向東,灣內(nèi)運河河口附近水交換率最高。工程后萊州灣整體年平均水交換能力與現(xiàn)狀相比增強約2.0%,半交換周期縮短11.0%;其中萊州灣西部水交換能力變化較小,增強約0.5%;萊州灣東部水交換能力變化較大,增強約2.7%。運河工程的實施對萊州灣水交換能力提升主要位于濰坊港防波堤至刁龍咀之間海域,作用范圍僅為河口向海5 km 的扇面海域。

需要指出的是,黃海海水經(jīng)由膠州灣、運河進入到萊州灣的海水水質(zhì)優(yōu)良,同時保證運河兩側(cè)無陸源污染物匯入的前提下,以目前設(shè)計的工況貫通運河,萊州灣東部運河河口附近水質(zhì)污染會略有緩解,但無法從根本上解決海洋生態(tài)環(huán)境問題,萊州灣污染治理的關(guān)鍵應(yīng)該是污染源頭的控制和管理。