基于水動力數值模擬的羅源灣水體交換周期及納潮量計算

陳志琦，逄 勇，張 倩，孫 凡，吳昌淦

(1.河海大學淺水湖泊綜合治理與資源開發教育部重點實驗室，江蘇 南京 210098；2.河海大學環境學院，江蘇 南京 210098)

隨著經濟的飛速發展，環境污染問題日益突出[1-2]。沿海城市在依靠優越的地理位置發展沿海經濟的同時，廢水和廢熱不斷排入海灣，影響了海灣的生態環境。羅源灣位于福州市東北部濱岸處，隨著羅源灣灘涂圍墾面積的增加和臨港工業及養殖業發展，羅源灣水質逐步惡化，生態系統的敏感性越來越強。近年來，福建省政府對羅源灣開展了退養工作，羅源灣水質雖有所改善[3-4]，但仍存在超標現象。水體交換周期及納潮量是衡量海灣自凈能力及海灣生存能力的重要指標[5-6]，研究海灣的納潮量及水體交換周期對海域的合理開發利用具有重要意義。

目前，對于水體交換周期的研究主要采用箱式模型、粒子追蹤模型和對流-擴散水質模型。蔣磊明等[7]根據FVCOM(finite-volume coastal ocean model)模擬結果計算了欽州灣的納潮量及水體交換周期。秦曉等[8]基于Mike21軟件對東山灣進行了水動力數值模擬，并計算得到東山灣納潮量及水體交換周期。周爭橋等[9]采用Lagrange質點追蹤方法對防城灣豐、枯水期水體交換能力進行了評價。1988年開始，國內學者對羅源灣開展了研究，但多集中在灣內魚貝類等水生生物的研究，對羅源灣水動力模擬較少。胡建宇[10]根據實測資料及單箱模型計算公式計算得出羅源灣納潮量較大，大潮期間納潮量可達9.6億m3，并假設羅源灣灣外海水與灣內海水均勻混合，得出羅源灣水體半交換周期約為17個潮周期。杜伊等[11-12]應用ECOMSED(estuarine, coastal and ocean modeling system with sediments)模擬了羅源灣水動力情況，并用示蹤粒子三維追蹤模塊分析羅源灣水體交換情況，結果表明，羅源灣內不同海區海水存留時間存在差異，主要由歐亞余流場結構差異所致；西南風有利于水體交換，東北風阻礙水體交換，高潮及落潮期間水體交換能力強，表層水體交換快，底層較慢。王勇智等[13]采用Mike3軟件，分析了羅源灣3個典型圍填海時期岸線與地形條件下的納潮量和水體交換能力變化，結果表明，1996年和2012年的海灣全潮平均納潮量與20世紀60年代相比分別減少了約20.59%和28.38%，灣內30 d的平均水體交換率則分別減小了19.17%和21.42%，水體半交換時間延長了約1.74 d和2.42 d。

納潮量是海灣環境評價的重要指標，近幾年，國內學者應用Mike等數值模擬軟件對國內外各個海域的納潮量進行了模擬分析。婁安剛等[14]運用二維潮波運動方程建立了丁字灣附近海域的水動力模型，模擬規劃實施前后水動力特性及納潮量的變化；陳靜等[6]應用Mike21軟件計算了大連灣及3個內灣的動態納潮量；王誠超等[5]基于樂清灣近50年水域面積的演變特征，對樂清灣納潮量的演變規律進行了研究，得出樂清灣納潮量存在遞減趨勢，納潮量對水域面積改變的響應關系有階段性特征。

目前，應用數值模擬軟件模擬羅源灣水動力情況，并計算其水體交換周期以及納潮量的研究較少，本文應用Mike21軟件模擬羅源灣水動力特性，計算和分析羅源灣納潮量及水體交換周期，可為羅源灣污染物總量控制提供基礎性資料。

1 研究區概況

羅源灣是一個半封閉港灣，口窄腹大。海灣面積約154 km2，位于東經119.60°～119.86°、北緯26.33°～26.49°之間，北部與羅源縣相連，南部與連江縣毗鄰，羅源縣縣城與羅源灣西北角相接。匯入羅源灣的地表徑流主要有羅源縣的起步溪、護國溪、南門溪，連江縣的鯉溪、牛柘溪等，同時，羅源灣沿海陸域因受地形自然切割的影響，形成眾多短促的獨自入海溪流。灣內有1個潮位站——跡頭站，羅源灣地理位置及地形如圖1所示。

圖1 羅源灣地理位置及地形

2 研究方法

主要采用Mike21軟件的水動力模塊和對流擴散模塊進行數值模擬。Mike21軟件的計算基于三向不可壓縮和Reynolds值均勻分布的Navier-Stokes方程，并采用靜水壓力的假定和Boussinesq假定[15]。其水動力模塊基本計算方程包括連續性方程和動量方程，對流擴散模塊采用二維對流擴散方程進行計算。本文根據模型計算的潮位值確定羅源灣納潮量，依據不同時刻保守物質質量濃度確定羅源灣水體交換周期。

2.1 納潮量計算方法

在一個潮周期內海灣所能接納的海水體積即為海灣的納潮量。為獲得精確的計算結果，在有限體積方法中，將納潮量概念應用到每個網格范圍內：

(1)

式中:P為納潮量；Si為第i個網格的面積；Hhi和Hli分別為第i個網格的高、低潮水位；n為所選取海域的網格總數。

2.2 水體交換周期計算方法

設羅源灣保守物質的初始質量濃度為1 mg/L，外灣質量濃度為0 mg/L，無其他源項匯入，保守物質衰減系數取0。灣內污染物通過水體交換向灣外輸移，當灣內保守物質質量濃度降至0.5 mg/L時，認為水體完成了半交換，此時所用的時間即為該海灣水體的半交換周期。水體交換率指某點保守物質質量濃度的減小值與初始質量濃度的比值：

(2)

式中:γ為交換率；Δρ為保守物質質量濃度的減小值，mg/L；ρ1為保守物質初始質量濃度，mg/L。

3 模型構建與參數率定

模型計算范圍是羅源灣灣口至羅源灣整個海域，采用三角形網格對羅源灣海域進行剖分，網格總數6 821，最小的三角形單元面積為2 321 m2。模型計算的時間步長為Δt= 300 s，計算時間為2019年9月1日至2019年12月31日，總時長為122 d。

初始水位取模型計算起始時刻邊界潮位均值，設為-1.7 m，起始時刻流速設為0 m/s，降雨、蒸發、風速和風向資料來源于國家氣象信息中心氣象監測站——寧德站(區站號58846)2019年逐日實測數據。模型開邊界采用潮位驅動，灣口處潮位由Mike21軟件Toolbox中的全球潮汐預測程序給出，該潮汐預測程序采用預測區域8個分潮的調和常數進行疊加獲得潮位，再根據計算域內潮位驗證情況進行調整。岸灘邊界的變化由計算格點水深來判別，當某個網格點的水深小于某一極限值(取0.001 m)，且其周圍網格點的水位均低于該網格點的水位時，該網格點被認為是“干單元”，不參與計算；當“干單元”周圍任意一點的水位高于該點水位時，“干單元”即轉化為“濕單元”，重新參與計算，模型計算區域網格劃分如圖2所示。

選取羅源灣跡頭站大小潮潮位對模型進行率定。小潮潮位率定采用2019年9月7—9日潮位資料，大潮潮位率定采用2019年10月27—29日潮位資料。模型對大潮潮位模擬的平均絕對誤差為26 cm，平均相對誤差為14%；模型對小潮潮位模擬的平均絕對誤差為14 cm,平均相對誤差為19%。模型率定結果如圖3所示。

圖2 羅源灣模型計算區域與網格劃分

(a)大潮

(b)小潮

余東等[16]研究表明，羅源灣可門水道處流速約為1.0 m/s，灣內流速約為0.5 m/s；漲潮時，灣外水流經灣口可門水道流向灣內后分為兩支，一支流向西北灣跡頭，另一支流向西南；落潮時，羅源灣海域流場是典型的往復流，與本文流場模擬結果基本相符。

根據率定結果得到模型灣內海域Smagorinsky系數為0.28，糙率為0.024～0.035，風拖曳系數為0.001 0～0.001 5。

4 模擬結果與分析

4.1 流場模擬

羅源灣海域流場變化模擬結果如圖4所示。

a.漲潮時，在灣外海域潮波作用下，大量海水涌入灣內，潮流在可門水道轉向西南進入羅源灣，灣口處流速較大，最大處可達1.24 m/s。水流進入灣口后，水域面積增大，流速降低，為0.60～1.24 m/s，羅源灣內部流速為0.20～0.60 m/s。漲潮時羅源縣城區附近海域以及將軍帽工業園區附近海域流速受潮波影響較小，低于0.20 m/s。總體來說，羅源灣灣口處水動力條件較好，北岸水動力條件較差，北岸流場受潮波影響程度小于南岸流場。

(a)大潮(落潮) (b)大潮(漲潮)

(c)小潮(落潮) (d)小潮(漲潮)

b.落潮時，潮流從羅源灣退至外海，流向與漲潮時基本相反，海水由北岸流經南岸，再從南岸流向灣口。灣內流速為0.26～0.90 m/s，可門水道內流速較大，最大可達1.56 m/s，水動力條件較好，北岸水動力條件劣于南岸。落潮時南北岸流場密集，大量海水從可門水道流出羅源灣。

4.2 水體交換周期

按照水體交換周期研究方法進行參數設置，計算可得15 d、30 d、45 d、60 d后羅源灣水體平均交換率分別為48.5%、69.3%、78.8%和85.5%，半交換周期約為15～16 d (表1)。由圖5可知，在相同時間內灣口及可門水道內水體交換率最大；灣中部水體交換率中等；灣頂處(羅源灣西北角)水體交換率最小。15 d后灣中部水體交換率在30%～80%之間,30 d后可達50%～90%;15 d后灣頂處水體交換率低于20%,30 d后交換率可達44%,60 d 后可達72%。這表明隨著時間的推移，受灣外海域潮波影響，灣內的保守物質濃度不斷降低，其中距離外海最近的灣口及可門水道水體交換能力最強，污染物濃度迅速降低。灣中部水體交換能力良好，而西北角灣頂處水體交換能力最弱，這主要是由于羅源灣西北角三面環陸,而整個海灣僅通過東南方可門水道與外界進行水體交換,導致此處受外海潮波影響較小，且流速較低，不利于污染物擴散。根據1998年胡建宇[10]對羅源灣水體半交換周期的研究，羅源灣水體半交換周期為17個潮周期，約為9 d；本文得到羅源灣水體半交換周期約為15 d，說明羅源灣海域更新能力減弱，半交換周期時間延長。

表1 羅源灣水體交換率

(a)15 d(落潮) (b)15 d(漲潮)

(c)30 d(落潮) (d)30 d(漲潮)

(e)45 d(落潮) (f)45 d(漲潮)

(g)60 d(落潮) (h)60 d(漲潮)

為定量分析灣內污染物交換情況，在灣口、灣中部、灣西部各取一點位(圖1)，分析其保守物質質量濃度的變化規律。由圖6可知，受漲落潮影響，A、B、C點保守物質質量濃度呈波動變化，灣口保守物質質量濃度受漲落潮影響最為明顯。位于灣口A點的保守物質質量濃度急劇降低后趨于穩定，3 d后A點的保守物質質量濃度低于0.5 mg/L。位于灣中部B點的保守物質質量濃度降低速率無明顯的突變過程，26 d后B點的保守物質質量濃度基本低于0.5 mg/L，50 d后B點的保守物質質量濃度低于0.2 mg/L。位于西北角C點的保守物質質量濃度的變化受漲落潮影響最小，C點的水體半交換周期約為43 d，60 d后C點的保守物質質量濃度約為0.3 mg/L，水體交換率僅為70%。

圖6 A、B、C點保守物質質量濃度隨時間變化曲線

4.3 納潮量

利用羅源灣模型模擬結果計算9—10月(秋季)的納潮量，羅源灣平均納潮量為5.73 億m3，大潮納潮量為7.07億m3，小潮納潮量為4.40億m3；大小潮時期納潮量差距較大，兩者比值為1.61∶1。根據胡建宇[10]的研究成果，1998年羅源灣平均納潮量為7.2億m3，大潮納潮量為9.6億m3，小潮納潮量為5.1億m3。與1998年相比，羅源灣納潮量有減少的趨勢，2019年平均納潮量減少了20.4%，大潮納潮量減少了26.4%，小潮納潮量減少了13.7%。根據杜伊等[11-12]對羅源灣納潮量的研究，2007年羅源灣平均納潮量為7.2億m3，與該結果相比，羅源灣納潮量仍呈減少的趨勢。

綜合考慮羅源灣海域形狀及水動力特點，將羅源灣分為Ⅰ區、Ⅱ區和Ⅲ區，分別計算其納潮量。Ⅰ區、Ⅱ區和Ⅲ區面積占羅源灣總面積的比例分別為16%、66%和18%。Ⅲ區位于羅源灣西部，受外海水動力條件影響較小，水體交換率較小，15 d后水體交換率低于20%，大潮潮差為5.55 m，小潮潮差為3.42 m，平均納潮量為0.95 億m3；Ⅱ區范圍較大，水體交換能力中等，大潮潮差為5.15 m，小潮潮差為3.28 m，平均納潮量為3.65 億m3；Ⅰ區包括灣口及可門水道，受外海區域潮波影響大，潮差較大，大潮潮差為6.34 m,小潮潮差為3.7 m，平均納潮量為1.13億m3。羅源灣秋季納潮量計算結果見表2。

表2 羅源灣秋季納潮量計算結果

在模型計算基礎上，分析秋季(9—10月)、冬季(12月)海灣納潮量變化。結果表明，冬季Ⅲ區、Ⅱ區和Ⅰ區的平均納潮量分別為秋季的87%、89%和86%，整個海灣冬季納潮量減少，為秋季的88%。

5 結 論

a.構建了羅源灣水動力數學模型，率定得到模型灣內海域Smagorinsky系數為0.28，糙率為0.024～0.035，風拖曳系數為0.001 0～0.001 5。

b.漲潮時，灣口處流速達到1.24 m/s，海水流過灣口后，流速降低為0.60～1.24 m/s，羅源灣內部流速為0.20～0.60 m/s；落潮時，海水流向與漲潮時基本相反，灣內流速為0.26～0.90 m/s，灣口處最大可達1.56 m/s。羅源灣灣口處水動力條件較好，北岸水動力條件較差，北岸流場受外海潮波影響程度小于南岸流場。

c.羅源灣水體平均半交換周期約15 d，15 d、30 d、45 d和60 d后羅源灣水體平均交換率分別為48.5%、69.3%、78.8%和85.5%，灣口及可門水道處水體交換率最大，交換能力最強，羅源灣西北角和西南處的灣頂處水體交換率最小。

d.秋季時，羅源灣平均納潮量為5.73億m3，大潮納潮量為7.07億m3，小潮時納潮量為4.40億m3。Ⅰ區平均納潮量為1.13億m3，占總納潮量的19.7%；Ⅱ區平均納潮量為3.65億m3，占總納潮量的63.7%；Ⅲ區平均納潮量為0.95億m3，占總納潮量的16.6%。

e.羅源灣西北角流速小，水體交換率低，納潮量有限，水動力條件較差。建議在羅源灣入海污染物通量管控過程中，綜合考慮海域自身水動力條件，合理控制羅源灣西北角入海污染物排放量。