有限體積海岸海洋模型模擬泉州灣人工島建設(shè)對水交換的影響

謝毅暉， 盧毅敏

(福州大學(xué) 數(shù)字中國研究院(福建)， 福建 福州 350003)

位于福建省泉州市的泉州灣是晉江、洛陽江匯合入海的半封閉海灣，北與莆田市湄洲灣相接，南與廈門市圍頭灣相鄰，東瀕臺灣海峽.泉州灣具有豐富的漁業(yè)資源及獨(dú)特的港口區(qū)位優(yōu)勢，對泉州及其周邊地區(qū)的發(fā)展起到不可或缺的重要作用.為了進(jìn)一步提高泉州灣的對外交通能力及大宗貨物集散的樞紐效應(yīng)，承接泉州灣后渚港區(qū)的功能轉(zhuǎn)換，替代泉州中心城區(qū)及其周邊的貨物集散運(yùn)輸功能，促進(jìn)海上絲綢之路出海通道的形成[1]，經(jīng)多方論證，泉州市政府決定建設(shè)秀涂人工島[2]，秀涂人工島工程位于泉州灣中部，石湖作業(yè)港區(qū)以北，通過圍?；靥钚纬桑斯u長4 480 m，寬800 m，工程用海3.4 km2，可形成陸域0.227 km2.

目前,泉州灣海域的數(shù)值模擬研究多停留于水動(dòng)力層次，關(guān)于泉州灣水交換能力的研究較少.林作梁等[3]基于有限體積海岸海洋模型(FVCOM)分析泉州灣三維潮汐及潮流特性.楊晨等[4]采用二維數(shù)學(xué)模型，對泉州灣潮流運(yùn)動(dòng)進(jìn)行數(shù)值模擬.秀涂人工島的建立不可避免地影響海域的水動(dòng)力特征，進(jìn)而影響水交換能力.水交換能力是描述海洋動(dòng)力過程的參數(shù)，它直接影響灣內(nèi)水體與外海的交換強(qiáng)度，進(jìn)而影響海灣水質(zhì).不同海域的對流輸運(yùn)及稀釋擴(kuò)散能力導(dǎo)致其自我凈化能力各不相同.對于自凈能力較弱的海域，海灣內(nèi)的污染物長期不能與外海水體進(jìn)行交換，導(dǎo)致污染物在海灣內(nèi)持續(xù)積累，造成富營養(yǎng)化等問題.

由于單一概念無法準(zhǔn)確地描述海域的水交換能力.因此，學(xué)者定義了不同概念定量描述水交換能力.Bolin等[5]提出“壽命”的概念，用于描述水質(zhì)點(diǎn)在研究區(qū)域經(jīng)歷的時(shí)間.Zimmerman[6]在此基礎(chǔ)上引入“滯留時(shí)間”，即水質(zhì)點(diǎn)離開研究區(qū)域所需的時(shí)間.Luff等[7]定義了“半交換時(shí)間”，即污染物濃度變?yōu)槌跏紳舛鹊囊话胨璧臅r(shí)間.歐拉彌散方法可定量分析海域的水交換能力，相較于其他方法，歐拉彌散方法在考慮對流的同時(shí)，還考慮了濃度的擴(kuò)散過程，能夠更加合理、準(zhǔn)確地反映海灣水體污染物的擴(kuò)散情況[8-14].本文采用有限體積海岸海洋模型對泉州灣的潮汐、潮流進(jìn)行模擬，定量分析泉州灣秀涂人工島的建設(shè)對該海域水交換能力的影響.

圖1 泉州灣區(qū)域圖Fig.1 Quanzhou Bay area map

1 有限體積海岸海洋模型及其應(yīng)用

有限體積海岸海洋模型是基于不規(guī)則三角形網(wǎng)格、有限體積的三維數(shù)值模型[15].該模型兼具有限元法和有限差分法的優(yōu)點(diǎn)，可充分保證模型的動(dòng)量、能量和質(zhì)量，具有更好的守恒性.

泉州灣區(qū)域圖，如圖1所示.模型研究區(qū)域?yàn)槿轂持苓吅Ｓ颍髌饡x江大橋，東至道仕嶼與尖峰嶼之間的水道，南起深滬灣南側(cè)，北至龍嶼水域，包含大、小墜門島、大山嶼等眾多島嶼和暗礁.模型研究區(qū)域在水平方向劃分為不重疊的非結(jié)構(gòu)三角形網(wǎng)格，并在晉江、洛陽江流域及大、小墜門島附近對網(wǎng)格進(jìn)行加密，網(wǎng)格的最小分辨率為110 m，開邊界節(jié)點(diǎn)共95個(gè)，平均節(jié)點(diǎn)間距約為1.6 km，區(qū)域網(wǎng)格點(diǎn)數(shù)為8 797個(gè)，單元格數(shù)量為16 452個(gè)；模型研究區(qū)域在垂直方向采用σ坐標(biāo)系，平均分為5層，外模時(shí)間步長為1 s，內(nèi)模時(shí)間步長為10 s.

泉州灣不規(guī)則三角形網(wǎng)格，如圖2所示.建島前、后不規(guī)則三角形網(wǎng)格只在人工島附近有所不同，其余區(qū)域完全一致.模型采用的海岸線地形數(shù)據(jù)來自美國國家環(huán)境信息中心(NCEI)，通過GEODAS軟件進(jìn)行提取.灣內(nèi)水深數(shù)據(jù)采用中國人民解放軍海軍海道測量局(NGD)發(fā)布的1∶35 000海圖，灣外近海開邊界附近的水深數(shù)據(jù)采用美國國家海洋和大氣管理局(NOAA)公布的全球地形水深數(shù)據(jù)，通過SMS軟件將岸線數(shù)據(jù)閉合生成不規(guī)則三角形網(wǎng)格，并將水深數(shù)據(jù)插值到三角形網(wǎng)格節(jié)點(diǎn)中.開邊界處潮位驅(qū)動(dòng)采用潮汐預(yù)測軟件OTPS預(yù)報(bào)的實(shí)時(shí)水位，模型主要考慮K1，O1，M2，S2等4個(gè)主要分潮的影響，開邊界95個(gè)節(jié)點(diǎn)處的時(shí)間序列潮位數(shù)據(jù)由Matlab軟件的TMD程序包生成.外部強(qiáng)迫場數(shù)據(jù)主要考慮風(fēng)、蒸發(fā)、降水和熱通量的影響，數(shù)據(jù)來自再分析資料CFSR，其時(shí)間分辨率為6 h，空間分辨率為0.2°.模型采用正壓模式，整個(gè)海域內(nèi)初始溫鹽場為常數(shù)，不考慮溫鹽的時(shí)空變化.模式運(yùn)行采用冷啟動(dòng)，即假設(shè)初始場流速從0開始.

(a) 建島前 (b) 建島后圖2 泉州灣不規(guī)則三角形網(wǎng)格Fig.2 Irregular triangular mesh of Quanzhou Bay

2 有限體積海岸海洋模型的驗(yàn)證

水動(dòng)力模擬的時(shí)間段為2016年9月1日至2016年10月30日，共60 d.模型從9月4日0時(shí)0分0秒開始輸出，10月30日0時(shí)0分0秒結(jié)束，每隔1 h輸出1個(gè)結(jié)果，時(shí)長為1 345 h.通過石湖(站點(diǎn)xt1)和浮山(站點(diǎn)xt4)兩個(gè)潮流站點(diǎn)測定潮位、流速和流向?qū)崪y值[16]，并與模型模擬值進(jìn)行對比驗(yàn)證.

2.1 潮位的驗(yàn)證

泉州灣潮位(h)驗(yàn)證圖，如圖3所示.由圖3可知：在潮位漲落趨勢上，潮位模擬值與實(shí)測值基本吻合；24 h內(nèi)出現(xiàn)兩次高潮和兩次低潮，最高潮位接近4 m，充分體現(xiàn)了泉州灣海域半日潮的特性；除個(gè)別點(diǎn)和時(shí)刻外，潮位過程、高潮位和低潮位的模擬值與實(shí)測值較為一致，模擬結(jié)果能夠較好地反映潮位的變化特征.

(a) 站點(diǎn)xt4 (b) 站點(diǎn)xt1圖3 泉州灣潮位驗(yàn)證圖Fig.3 Tidal level verification map of Quanzhou Bay

2.2 流速、流向的驗(yàn)證

泉州灣潮流流速(v)驗(yàn)證圖，如圖4所示.泉州灣潮流流向(θ)驗(yàn)證圖，如圖5所示.

(a) 站點(diǎn)xt1 (b) 站點(diǎn)xt4圖4 泉州灣潮流流速驗(yàn)證圖Fig.4 Tidal current velocity verification map of Quanzhou Bay

(a) 站點(diǎn)xt1 (b) 站點(diǎn)xt4 圖5 泉州灣潮流流向驗(yàn)證圖Fig.5 Tidal current direction verification map of Quanzhou Bay

由圖4，5可知：流速和流向的模擬值與實(shí)測值基本吻合，轉(zhuǎn)流時(shí)刻基本一致，總體變化趨勢一致，潮流場模擬的重現(xiàn)度較好；最大流速和流向的模擬值與實(shí)測值略有偏差，這是由于流速、流向等數(shù)據(jù)在三角形網(wǎng)格的中心原點(diǎn)位置，劃分三角形網(wǎng)格時(shí)未能足夠縮小網(wǎng)格步長，造成模擬值與實(shí)際值略有偏差.由此可知，泉州灣海域潮流場的計(jì)算結(jié)果具有可靠性，模型能夠較好地模擬泉州灣海域的潮流特征.

3 水環(huán)境的變化

3.1 水動(dòng)力特征的變化

泉州灣潮流流速的分布情況，如圖6所示.由圖6可知：建島前、后流速的差異主要體現(xiàn)在灣口處，灣內(nèi)流速差異不明顯；漲急時(shí)，潮流由東南部外海涌入泉州灣，灣內(nèi)整體流向?yàn)槲飨蚝捅毕騕3]；漲憩時(shí)，潮流大致以大、小墜門島及大山嶼形成分界線，分界線左側(cè)的潮流流向與漲急時(shí)相似，以西向和北向?yàn)橹?，分界線右側(cè)的潮流沿東朝灣外流出；落急時(shí)，灣內(nèi)潮流呈現(xiàn)朝灣外東北部流動(dòng)的態(tài)勢，灣內(nèi)整體流向?yàn)闁|向和南向，與漲急時(shí)相反；落憩時(shí)，潮流大致以大、小墜門島及大山嶼形成分界線，分界線左側(cè)的潮流流向與落急時(shí)相似，以東向和南向?yàn)橹?，分界線右側(cè)的灣外潮流沿東南朝北涌入.

(a) 漲急(建島前) (b) 漲憩(建島前) (c) 落急(建島前) (d) 落憩(建島前)

(e) 漲急(建島后) (f) 漲憩(建島后) (g) 落急(建島后) (h) 落憩(建島后)圖6 泉州灣潮流流速的分布情況Fig.6 Tidal current velocity distribution of Quanzhou Bay

建島前、后的潮流流速差(Δv),如圖7所示.由圖7可知：秀涂人工島建成后，泉州灣內(nèi)秀涂通港引橋以東區(qū)域漲急、漲憩、落急、落憩4個(gè)時(shí)刻的流速差約為0.1 m·s-1，秀涂人工島的北側(cè)北烏礁水道一帶及南側(cè)石湖港區(qū)的流速差達(dá)到0.2 m·s-1，這是因?yàn)槿斯u建立后，原本直接通過泉州灣中部的潮流在秀涂島附近分為南、北兩支分流，導(dǎo)致流速差明顯；由于秀涂人工島對潮流的阻攔，崇武附近海域及泉州灣口的流速差明顯.

(a) 漲急 (b) 漲憩

(c) 落急 (d) 落憩圖7 建島前、后潮流流速差的對比Fig.7 Tidal current velocity difference before and after island construction

3.2 潮致余流的變化

潮致余流是指在潮汐運(yùn)動(dòng)周期中，海域空間中某點(diǎn)的潮流速度的平均值，它表明該點(diǎn)在潮周期過程中的平均遷移趨勢，其強(qiáng)弱往往受到復(fù)雜多變的地形和蜿蜒盤旋的岸線的影響[17].因此，雖然相較于潮流而言，潮致余流極小，但它對物質(zhì)輸運(yùn)至關(guān)重要.

潮致余流(vE)取模型穩(wěn)定后25 h的結(jié)果進(jìn)行計(jì)算,計(jì)算公式[11]為

(1)

式(1)中：T為時(shí)間；i為網(wǎng)格數(shù).

建島前、后潮致余流的對比,如圖8所示.

(a) 建島前 (b) 建島后圖8 建島前、后潮致余流的對比Fig.8 Comparison of tidal residual currents before and after island construction

由圖8可得以下3點(diǎn)結(jié)論.

1) 建島前，潮致余流一般為0.2 m·s-1，建島后，泉州灣潮致余流減小較為明顯，尤其是灣口海域；建島前，最大潮致余流普遍超過0.5 m·s-1，而建島后最大潮致余流不到0.5 m·s-1.

2) 受地形、水深、岸線等因素的影響，建島前、后潮致余流的強(qiáng)流區(qū)主要分布在石湖港區(qū)北側(cè)、大、小墜門島等附近水道，在崇武沿岸、浮山周邊形成較強(qiáng)的岬角渦旋，在小墜門島附近，形成封閉的反氣旋式余流渦旋，這與文獻(xiàn)[3]的實(shí)驗(yàn)結(jié)果相符.

3) 泉州灣的潮致余流場呈一定程度的分叉趨勢，即一股余流沿著浮山南側(cè)、崇武沿岸及赤星礁一側(cè)經(jīng)過北烏礁水道向西流進(jìn)泉州灣；另一股余流主要沿著石湖北側(cè)、北烏礁、大、小墜門島水道內(nèi)朝東流出泉州灣.兩股余流呈現(xiàn)不同的運(yùn)動(dòng)特征，在小墜門島的西南側(cè)及圍頭東北側(cè)海域各形成封閉的反氣旋式環(huán)流.余流場與潮流場強(qiáng)流區(qū)的分布基本一致，證明了潮致余流與潮流速度具有密切的關(guān)系[18].

3.3 納潮量的變化

納潮量是指一個(gè)海灣容納潮水的最大體積，用于表征一個(gè)半封閉海灣的生存能力[19-20].納潮量的變化對海灣水動(dòng)力特性的影響不可忽視，對維持海灣的生態(tài)環(huán)境平衡具有重要作用.一般而言，納潮量越大，表明灣內(nèi)水體與灣外水體的交換越快，自凈能力越強(qiáng)，水交換能力越強(qiáng)；反之則相反.通常納潮量W的計(jì)算公式[20]為

(2)

式(2)中：S1，S2分別為高水位和低水位時(shí)的水域面積；Δh為高、低潮位之間的差值.

然而，式(2)一般將海灣當(dāng)成一個(gè)理想的立方體，而泉州灣島嶼眾多、岸線蜿蜒曲折，不能簡單地視作理想立方體.因此，采用納潮量的實(shí)際計(jì)算公式[21]，即

(3)

式(3)中：n為網(wǎng)格總數(shù)；Si為第i個(gè)網(wǎng)格的面積；h1，i，h2，i分別為第i個(gè)網(wǎng)格的高潮位和低潮位.

為了研究建島前、后潮流強(qiáng)度對納潮量的影響，采用大、小潮期間的潮流動(dòng)力驅(qū)動(dòng)建島前、后的模型，并計(jì)算相應(yīng)的納潮量.建島前、后大潮期間納潮量分別為2.789 7×109，2.656 0×109m3，建島后大潮期間納潮量比建島前減少了0.133 7×109m3，變化率為4.79%;建島前、后小潮期間納潮量分別為1.526 4×109，1.372 3×109m3，建島后小潮期間納潮量比建島前減少了0.154 1×109m3，變化率為10.09%.由此可得以下3點(diǎn)結(jié)論.

1) 建島前、后大潮期間的納潮量都約為小潮期間納潮量的兩倍，這可能是由于大、小潮期間的潮位差異造成的，大潮期間的潮位差異約為小潮期間的兩倍，在總面積保持不變的情況下，納潮量的差異就等同于潮位差異，秀涂人工島的建立后并未改變兩者之間的倍數(shù)關(guān)系.

2) 建島后，隨著海域面積減少了3.22 km2，納潮量也呈現(xiàn)減少趨勢.納潮量的減少會削弱泉州灣內(nèi)水體與灣外水體的交換能力，繼而弱化灣內(nèi)水體的自凈能力，使污染物在灣內(nèi)堆積，無法及時(shí)、有效地與外海水體混合，這將增大灣內(nèi)水體遭受污染的可能性，影響灣內(nèi)的生態(tài)平衡.

3) 相較于大潮期間納潮量，小潮期間納潮量的減小程度更大，這說明秀涂人工島的建立對小潮的影響大于大潮，小潮期間灣內(nèi)水體與灣外水體的交換能力較弱，更容易遭受污染.

3.4 水交換能力的變化

采用Luff等[7]定義的“半交換時(shí)間”進(jìn)行泉州灣水交換能力的定量研究，相關(guān)計(jì)算公式[22]為

(4)

式(4)中：t表示平均滯留時(shí)間；C(t)表示當(dāng)前時(shí)刻海域內(nèi)的污染物濃度；C0表示t0初始時(shí)刻的污染物濃度；當(dāng)C(t)為初始污染物濃度的50%時(shí)，半交換時(shí)間th=t-t0.

海域污染物濃度Ca的計(jì)算公式[23]為

(5)

式(5)中：Ci,j表示第i個(gè)網(wǎng)格，第j層的污染物濃度；hj為第j層的高度.

在進(jìn)行污染物擴(kuò)散計(jì)算時(shí)，采用歐拉彌散方法計(jì)算灣內(nèi)污染物濃度，將灣內(nèi)污染物濃度設(shè)置為 1(以此為參照值，無單位)，灣外水體污染物初始濃度設(shè)置為 0.采用文中模型提供的水動(dòng)力特征條件進(jìn)行污染物擴(kuò)散的數(shù)值模擬，當(dāng)模型的水動(dòng)力場數(shù)值穩(wěn)定時(shí)，進(jìn)行56 d的污染物釋放和擴(kuò)散模擬.

水交換率δ=(1-Ca)×100%,由此計(jì)算灣內(nèi)水體的水交換率(以第15,35,55天數(shù)據(jù)為例)，結(jié)果如表1所示.

表1 建島前、后灣內(nèi)水體水交換率的變化Tab.1 Changes of water exchange rate in bay before and after island construction

建島前、后污染物濃度的變化，如圖9所示.

(a) 第15天(建島前) (b) 第35天(建島前) (c) 第55天(建島前)

(d) 第15天(建島后) (e) 第35天(建島后) (f) 第55天(建島后)圖9 建島前、后污染物濃度的變化Fig.9 Changes of pollutant concentration before and after island construction

由表1和圖9可得以下4點(diǎn)結(jié)論.

1) 污染物濃度從灣口到灣頂呈遞減趨勢，灣口處水動(dòng)力較強(qiáng),污染物濃度下降較快;灣頂處水深較淺，區(qū)域的半封閉使水動(dòng)力較弱，污染物濃度下降較慢.

2) 當(dāng)污染物釋放第15天時(shí)，建島前泉州灣內(nèi)以石湖港區(qū)為界，右側(cè)大、小墜門島附近的污染物濃度基本降至0.2以下，分界線左側(cè)污染物濃度分別朝晉江和洛陽江兩側(cè)遞減，水交換率達(dá)63.60%；隨著秀涂人工島建立，石湖港區(qū)分界線右側(cè)的污染物濃度普遍約為0.3，浮山一側(cè)的污染物濃度最大達(dá)到0.4，在石湖與人工島分界線上，由于人工島的存在，兩側(cè)航道為灣內(nèi)與灣外水交換的必經(jīng)之道，污染物濃度約為0.2，并以此朝晉江和洛陽江兩個(gè)灣頂水域逐漸遞減，灣內(nèi)水交換率達(dá)55.56%.

3) 當(dāng)污染物釋放第35天時(shí)，建島前、后泉州灣口部分海域的污染物濃度基本都在0.2以下，并大致沿著沙帽礁-南烏礁-北烏礁連接成的分界線向左側(cè)逐漸遞減，灣內(nèi)大部分海域的污染物濃度降至初始濃度的50%以下, 僅在一些灘涂、港區(qū)內(nèi)部因水動(dòng)力較弱,污染物濃度仍舊較高，建島前灣內(nèi)水交換率為77.44%，建島后水交換率為71.01%.

4) 當(dāng)污染物釋放第55天時(shí)，灣內(nèi)污染物濃度變化放緩，建島前灣內(nèi)水交換率為79.85%，建島后灣內(nèi)水交換率為72.43%，之后的濃度基本可視為定常.

圖10 建島前、后污染物的濃度差(第56天)Fig.10 Concentration difference of pollutant before and after island construction (day 56)

將建島后第56天的污染物擴(kuò)散結(jié)果減去建島前第56天的結(jié)果，得到污染物的濃度差(ΔC)，如圖10所示.由圖10可知：相較于建島前，建島后泉州灣海域的污染物濃度呈現(xiàn)增大趨勢，即擴(kuò)散減慢，水交換率降低；泉州灣北部及西南海域，特別是洛陽江流域和金嶼的污染物濃度差較大，前后變化量差值超過0.15，在臨近泉州灣灣口處，污染物濃度差略微出現(xiàn)遞減態(tài)勢，北烏礁附近海域達(dá)到最小值.

泉州灣建島前、后的半交換時(shí)間及其差值(Δth)，如圖11所示.

(a) 半交換時(shí)間(建島前) (b) 半交換時(shí)間(建島后) (c) 半交換時(shí)間差值 圖11 建島前、后的半交換時(shí)間及其差值Fig.11 Half exchange time and its difference value before and after island construction

由圖11(a),(b)可知：泉州灣水體半交換時(shí)間與污染物濃度差分布大體類似，泉州灣海域水體的半交換時(shí)間從灣口向?yàn)稠斨饾u增加；北部和西南部的半交換時(shí)間較長，大部分海域的半交換時(shí)間長達(dá)30 d以上；中部和東部的半交換周期較短，多為10 d以下.由圖11(c)可知：建島后，泉州灣北部洛陽江流域和西部金嶼靠近陸地部分區(qū)域的水體半交換時(shí)間出現(xiàn)了不同程度的增加，普遍增加了10～15 d，部分海域的半交換時(shí)間增加較多，增加了25 d以上.

泉州灣污染物濃度變化曲線，如圖12所示.由圖12可知：建島前，泉州灣內(nèi)污染物濃度約7 d后即可降至初始濃度的50%；建島后，污染物濃度降至初始濃度的50%大約需10 d，較建島前約增加了3 d.

圖12 泉州灣污染物濃度變化曲線Fig.12 Variation curves of pollutant concentration of Quanzhou Bay

因此，秀涂人工島的建立在一定程度上影響了泉州灣污染物的對流和擴(kuò)散，灣內(nèi)水體與灣外水體的水交換時(shí)間變長.

4 結(jié)論

1) 秀涂人工島的建設(shè)主要對潮流流速產(chǎn)生影響.泉州灣大部分海域的流速約減小0.1 m·s-1；在秀涂人工島的北側(cè)北烏礁水道一帶及南側(cè)石湖港區(qū)之間，流速差達(dá)到0.2 m·s-1；由于人工島對潮流的阻攔，使灣口及崇武附近海域的流速差較為明顯.

2) 秀涂人工島對泉州灣石湖港區(qū)人工島連線以西大部分區(qū)域的余流結(jié)構(gòu)及大小影響并不顯著,強(qiáng)流區(qū)主要分布在石湖港區(qū)北側(cè)、大、小墜門島等附近水道內(nèi)，在泉州灣灣口的余流減少最為明顯.此外，在小墜門島附近，還形成了封閉的反氣旋式余流渦旋.

3) 建島前、后大潮期間納潮量均約為小潮期間納潮量的兩倍，建島后的海域面積減少3.22 km2，由此導(dǎo)致納潮量呈減少趨勢，大潮期間納潮量減少0.133 7×109m3，小潮期間納潮量減少0.154 1×109m3.小潮期間納潮量減少程度較大，變化率為10.09%，表明小潮期間更容易遭受污染威脅.

4) 建島后，泉州灣整體水交換能力減弱.建島前，灣內(nèi)水體半交換時(shí)間約為7 d，建島后，灣內(nèi)水體的半交換時(shí)間增加約3 d.秀涂人工島在一定程度上降低了泉州灣與外海之間的水交換能力，使污染物在灣頂堆積，濃度稀釋速度減緩.

5) 在實(shí)施海洋工程前，應(yīng)將其對水交換能力的影響納入考慮，以避免因圍填海建島工程造成惡劣的環(huán)境影響.