閩江河口三維潮流和余流特征及污染物運動軌跡的數值模擬

林金城， 肖桂榮， 林建偉

(1. 福州大學 數字中國研究院(福建)， 福建 福州 350108；2. 福建省水產研究所， 福建 廈門 361012)

閩江是中國福建省最大獨流入海(東海)河流，全長559 km，是沿海典型的山溪性河流[1].閩江河口屬于強潮三角洲型河口，閩江東流入海，受瑯岐、粗蘆、川石、壺江諸島阻隔，形成了五口入海的復雜河網[2].隨著閩江流域經濟社會的不斷發展，在沿岸地區形成了密集的居住區，大量的污染物隨著閩江東流入海，給閩江河口海域的生態環境帶來巨大的壓力.潮流是河口海灣物質輸運的動力基礎，對污染物的遷移擴散起著決定性的作用.余流則體現了水體的輸運和交換過程，與河口內物質的長期輸運、擴散、沉積等有著密切聯系[3].利用Lagrangian粒子示蹤法模擬水體中污染物的運動過程，通過計算水體中示蹤粒子的運動軌跡能夠直觀地模擬污染物在水體中的三維運移軌跡、滯留時間及歸宿等[4].

目前，已有一些針對閩江河口潮流及污染物運移規律的數值模擬研究.季杜鑫等[5]通過實測閩江河口的水文、泥沙資料等，指出閩江河口的潮汐特征為漲潮歷時短、退潮歷時長，且以往復流為主；劉梅冰[6]建立閩江下游河道一維水動力水質模型，對閩江下游水流、污染物運動規律及影響因素進行研究；文獻[7-8]通過建立閩江河口二維水動力模型，模擬閩江河口主要污染物的時空分布特征及運動變化規律；湯軍健等[9]選用二維淺水方程組模擬閩江河口的潮流場，建立拉格朗日質點跟蹤方法，近似模擬閩江河口海域的泥沙輸運特征；夏澤宇等[10]基于非結構網格、有限體積的近岸、河口海洋模型(FVCOM)建立閩江河口三維潮流數值模型，討論海底摩阻系數的選取，并分析閩江河口的水動力特征.

學者們的工作雖然較好地揭示了閩江河口潮流運動及污染物運移特征，但仍有不足之處：早期研究僅限于一維與二維模型，不能較好地模擬河口內潮；絕大部分研究忽視了閩江河口廣闊的潮間帶，未考慮漫灘效應，而漫灘過程的引入對閩江河口潮流速度的模擬至關重要[11]；大部分研究側重于對閩江河口潮流場的模擬，鮮有對余流特征及污染物粒子運移規律的研究.因此，深入研究閩江河口三維潮、余流及污染物運動規律特征，充分認識閩江河口物質輸運的動力基礎，對于閩江河口環境污染的控制和治理具有重要的指導意義.

本文選用三維水動力模型，即河口、陸架和海洋沉積物(ECOMSED)模型[12]，對閩江河口進行三維斜壓潮流數值模擬，將干濕網格判別技術引入潮汐潮流的漫灘過程，在模型驗證良好的基礎上，分析研究閩江河口的潮流、余流等特征，深入了解閩江河口的水動力狀況及物質輸運的動力基礎；采用Lagrangian粒子示蹤技術，模擬閩江河口保守污染物運移過程，以揭示污染物在閩江河口的運移規律.

1 模型的構建與驗證

國內外學者基于不同模型對潮流運動與物質輸運進行研究[11，13-18].ECOMSED模型是以普林斯頓海洋模型(POM)[19]與三維河口、海岸與海洋模型(ECOM)[20]為基礎發展起來的較為成熟的淺海三維水動力學模型，適用于河口及近岸海域的海洋模擬.該模型主要包括水動力模塊、沉積物輸運模塊、風浪模塊、熱通量模塊、水質模塊和粒子追蹤模塊；模型采用模態分離技術，考慮垂向分層.文中模型是基于ECOMSED的水動力模塊和粒子追蹤模塊建立的.模型的主要控制方程詳見文獻[21].

1.1 動邊界處理

閩江河口具有大面積的淤泥質潮灘，潮間帶上的水陸邊界隨著潮漲、潮落不斷變化，使計算區域不斷發生改變.然而，ECOMSED模型采用固定邊界技術，無法模擬漫灘過程，使計算區域隨著潮漲、潮落而發生改變，因此，引入動邊界處理技術，即干濕網格判別法.在若干時間步長內，對水深做一次判斷，若發現網格點總水深H小于某一臨界水深Hdry，則認為該水位網格為干網絡；若網格總水深H大于臨界水深Hdry，則認為其為濕網格.在模型運行過程中，當判斷為濕網格時才參與計算，而干網格不參與計算，且其流速應為零.臨界水深的大小可根據潮間帶坡度和時間步長確定，取Hdry=0.4 m[13].

1.2 數值模擬與驗證

考慮閩江河口河網復雜、岸線曲折、灘涂面積大等特點，模型的經緯度范圍為119°23′～119°50′E，25°55′～26°18′N.閩江河口模型的計算網格、水深示意圖，如圖1所示.模擬采用正交矩形網格，網格分辨率為200 m×200 m，共有網格數量222×208個，垂向采用sigma坐標，等間距分為11層.底摩擦系數IBFR根據文獻[10]取為0.007，再微調至恰當值；底粗糙系數BZO為0.006；水平對流擴散參數NHORCO取0.1；垂向紊動參數LUMO取1×10-6m2·s-1.內模時間步長為5 s，內外模時間步長分裂比為5.模擬采用冷啟動條件，不考慮波浪及風應力的影響.模擬時間為2014年6月11日0:00至2014年6月25日0:00，考慮模型穩定所需時間，提取2 d后的數據作為驗證數據.

(a) 計算網格 (b) 計算域水深

計算水域地形數據采用海軍航保部1∶30 000的閩江口海圖(圖號為13991)與1∶20 000的金牌門至馬尾海圖(圖號為13992A，13992B)，經數字化獲得網格水深后，利用內插方法進行計算.

1.2.1 邊界條件 模型通過潮汐調和常數法，取S2，M2，N2，K1，P1，O1這6個分潮，采用美國俄勒岡州立大學建立的中國海區域潮汐潮流模型的分潮調和常數，在最小二乘法的規則下，擬合TOPE/Poseidon and Jason沿軌跡平均數據和拉普拉斯潮汐方程獲得的模擬結果.

在閩江河口，模型采用開邊界處理，在河流相應的網格上加入2014年6月閩江平均徑流量形成的“源”，約為3 770 m3·s-1，以反映閩江入海流量.

1.2.2 模型的驗證 選取2014年6月13-24日的實測資料驗證模型.潮位及潮流的測點位置，如圖2所示.潮流觀測設4個站點(P1～P4)；潮位觀察設2個站點(T1和琯頭).

圖2 潮位及潮流的測點位置

對6個觀測站連續觀測的潮位、潮流資料進行對比驗證，結果如圖3所示.圖3中：h為潮位；v為流速；α為流向.由圖3可知：各觀測站的潮位、流速、流向都與實測結果較為吻合.模型驗證結果，如表1，2所示.表1中：vf，ve分別為漲潮、落潮的平均流速；ε1，ε2分別為漲潮、落潮平均流速的觀測值和模擬值的偏差；αf，αe分別為漲潮、落潮的平均流向；E1，E2分別為漲潮、落潮平均流向的觀測值和模擬值的偏差.表2中：hh，hl分別為高、低潮位；δ1，δ2分別為高、低潮位觀測值和模擬值的偏差.由表1，2可知：P1～P4潮流觀測點的漲、落潮流速的偏差值較小，相對誤差均小于10%，其對應的流向誤差小于10°；琯頭、T1潮位站平均高、低潮位的偏差值在±0.1 m以內.表1和表2的結果均符合技術規范[22]對潮流場的驗證要求.綜上可知，模型采用的物理參數和計算參數基本合理，計算方法可靠，能夠模擬閩江河口內潮波運動特性.

(a) 琯頭站點 (b) T1站點

表1 P1～P4潮流站流速、流向的驗證結果

表2 T1和琯頭站潮位的驗證結果

2 試驗結果與討論

2.1 閩江河口潮流場數值模擬

閩江河口大潮漲急和落急的表層流速矢量圖，如圖4所示.圖4中：vf，max，ve，max分別為漲潮、落潮的最大流速.閩江河口潮流以半日潮為主，潮流運動形式以往復流為主.潮流受閩江徑流作用影響，落潮流速明顯大于漲潮流速.潮流特征受地形變化影響較大，在外海10 m水深等值線與潮流流速呈現出的分界線相一致.

(a) 漲急 (b) 落急

由圖4(a)可知：外海潮波向河口內陸傳遞，沿水道流速逐漸增大，在10 m水深等值線可見明顯1.0 m·s-1流速分界.在漲潮過程中，外海潮波受粗蘆島、川石島與瑯岐島阻隔分叉，分別從烏豬水道、熨斗水道、川石水道和梅花水道進入向福州近岸海域傳播.其中，漲潮流在川石水道(深槽)的流場較強，由于水道水深相較周圍深，因此，流速增大，最大漲潮流速約為1.8 m·s-1；在烏豬水道上的最大漲潮流速約為1.4 m·s-1；在熨斗水道上的最大漲潮流速約為1.3 m·s-1；在金牌門峽道口，匯攏了從烏豬水道、熨斗水道、川石水道的潮流，水道狹窄且深，產生高流速區，流速最大可達1.7 m·s-1，向長門水道推進；由于梅花水道區域淤積嚴重、水深較淺，潮流流場較弱，流速約為0.7 m·s-1；南、北兩支水道主流在亭江區域匯攏，沿著河岸線朝西南方向向閩江上游推進.

由圖4(b)可知：落潮時，落潮方向與閩江徑流方向一致，與漲潮流方向相反流出閩江河口，最大流速發生在長門水道與川石水道之間，可達1.8 m·s-1，落潮流速大于漲潮流速.落潮時，閩江河口與外海未見明顯整齊的流速分界線，流速分界主要表現在沿水道深槽處.烏豬水道、熨斗水道的深槽流速比漲潮時刻更快，且作用明顯.落潮時，潮流整體偏向東南處流動.東北部海域在漲落潮時流速相差不大，而東南處海域落潮流速稍大于漲潮時流速，表明東北部海域比東南部海域的水體交換能力弱，容易滯留污染物；污染物更容易向東南部海域輸運.

通過閩江河口4個斷面(D1～D4)大潮期間的流速垂向剖面圖，分析其三維潮流特性，如圖5所示.圖5中：H為水深.由圖5可知：漲急時，在長門水道至川石水道區域，由于航道水深遠深于近岸兩側，D1～D4斷面流速從表層至中層均保持較高流速，且斷面中間流速高于兩側，為1.2～1.8 m·s-1，隨著水深的增加，流速逐漸減小且減小幅度較小，表層流速比底層大約0.5 m·s-1；由于地形變化，梅花水道開口處流速較為平緩，水深較淺，總流速明顯低于北支水道，而在水深較深的水道，流速則明顯增大，表明流速受地形影響較大；北支水道流速均大于南支水道流速.

(a) D1斷面漲急 (b) D1斷面落急

落急時，D1～D4斷面流速自表層隨水深的增大而逐漸減小，呈現出明顯的分層現象，垂向分布較為均勻.受徑流作用影響，落急時表層流速大于漲急時流速，但中層流速小于漲急時流速.南支水道在落潮時流速為0.43～1.40 m·s-1，流速減小的幅度較大，表層流速比底層大約0.80 m·s-1.

由于南北支水道水深地形差異，D1～D4斷面北支水道的漲落流速均大于南支水道.落急時流速呈現明顯的分層現象，流速隨水深的增大而減小，流速減小的幅度較大.南支水道漲潮流速大于落潮流速.

2.2 閩江河口余流場數值模擬

余流是指海流剔除掉周期性潮流運動后剩下的部分，與河口污染物的長期輸移有著密切關系.余流分析是研究海岸帶物質輸移方向的重要手段[23].采用歐拉余流計算方法，取兩個半日潮周期(約25 h)，計算大潮表層、中層、底層的余流流速.歐拉余流場分布示意圖，如圖6所示.圖6中：vr為余流流速.

(a) 表面余流

由圖6(a)可知：閩江河口大潮表層余流流向與閩江河口水道走向一致，以徑流為主要動力，沿河道向外海傳遞，以落潮余流為主；由于閩江徑流輸入影響，歐拉余流流速較大的區域主要在長門-川石水道(航道深槽)；在熨斗水道島群周圍出現較強余流區，余流流速約為0.40 m·s-1；烏豬水道水深較淺是弱余流區，隨著向外海流動，水深增大，河道變寬，余流流速提高至0.27 m·s-1；在南、北兩支水道上，余流強度南北不對稱，梅花水道余流強度隨著水深變淺，余流流速由0.40 m·s-1減弱至0.15 m·s-1，而北支水道的平均余流流速則保持在0.50 m·s-1左右；川石水道與梅花水道流出的余流在瑯岐島東側近岸區域相互作用，產生一條狹長的滯留區域；粗蘆島東北側及熨斗島東側為水深較淺的灘涂區域，余流流速較低，出現大面積的滯留區域.

由圖6(b)可知：中層余流流速弱于表層余流流速，仍沿河道向外海傳遞，以落潮余流為主；不同于表層，南支梅花水道中層余流比北支長門水道余流強度大，主要發生在猴嶼鎮東側、粗蘆島、熨斗島南側水道，余流流速約為0.28 m·s-1；在水深較淺的灘涂區域，余流流速明顯降低，僅為0.04 m·s-1左右；由于閩江徑流輸入作用，馬尾至閩安附近水道的中層余流場與表層基本一致.

由圖6(c)可知：底層余流比中層余流弱，但余流基本趨勢差別明顯，表現出以漲潮余流為主導地位；在馬尾附近水道，余流流向為東北方向，但在亭江附近水道，余流流向為西南方向，南、北兩道余流在閩安附近相互作用；長門水道深槽沒有出現較明顯的強余流，底層余流方向為西南方向；烏豬水道與熨斗水道的余流由外海傳入，流向為西南方向；在瑯岐島東側以及熨斗島東側附近海域，底層余流形成余流渦旋的趨勢；金牌門峽道窄口附近為強余流區，余流流速最大為0.35 m·s-1，流向為西向；梅花水道底層余流整體上朝東向，但在梅花鎮西北側余流流向表現為漲潮向西.

閩江河口的余流呈明顯的層化現象，余流流速由表層至底層逐層遞減，表層和底層余流在長門水道、川石水道等區域表現出不同的流向趨勢.在亭江附近、北支水道，表層余流表現為落潮流，底層余流則表現為漲潮流.在梅花水道，表層、中層余流整體上表現為落潮流，底層余流在猴嶼鎮東側表現為落潮流，但在梅花鎮西北側余流表現為漲潮流.在東北部海域，余流場整體表現較弱，余流流速僅有0.02～0.03 m·s-1，表明污染物在此容易滯留，不易向外海繼續遷移；而在東北部海域，表層余流則較強，余流流速約為0.15 m·s-1，表明污染物在此容易向外海輸移擴散，這與節2.1潮流場分析結果一致.

2.3 Lagrangian粒子運移軌跡模擬

2.3.1 粒子源設置 從排污口空間分布上看，閩江下游排污口主要集中在福州市區北支臺江與南支烏龍江；在閩江河口，長樂區梅花鎮附近有2個排污口，如圖7所示.

圖7 閩江河口主要排污口位置示意圖

為研究污染物的運動規律，分別在潭頭、梅花排污口與研究區域臺江、烏龍江的開邊界處設置模擬粒子釋放點源，通過Lagrangian粒子示蹤法模擬計算粒子的運動軌跡.由于粒子的運動軌跡不僅與釋放點位置有關，還與釋放時刻有關，研究中選擇在漲憩和落憩時刻釋放.閩江河口地形復雜，五口入海，且單個粒子運動具有較強的隨機性，因此，通過對污染物粒子群進行數值模擬，從而更準確地研究污染物的整體運移規律.梅花鎮附近的排污口數目較少，僅有2個(潭頭排污口與梅花排污口)，分別將粒子釋放在表、中、底3個水層，每層釋放10個粒子；在研究區域臺江與烏龍江的開邊界處，于水體表層分別同時釋放200個粒子.

2.3.2 粒子運移軌跡 漲憩和落憩時刻，分別在潭頭排污口與梅花排污口2個釋放點的水體表層、中層、底層釋放粒子，每10 min間隔的粒子運動軌跡圖，如圖8所示.由圖8可知：粒子在2 d內的運動軌跡基本符合閩江河口內流場狀況，受往復潮流的影響，粒子在近岸海域做往復運動，最終在徑流與潮流的作用下向東南方向遷移至外海.

(a) 漲憩時刻釋放的粒子

2個排污口處的粒子運動軌跡粒子均向東南方向遷移.潭頭與梅花排污口釋放點處，表層的粒子運動軌跡比中層、底層的粒子簡單，在漲憩時刻釋放的粒子未做過多往復運動，往東南方向海域向下遷移，離開閩江河口，而落憩時刻釋放的粒子，受漲潮流影響，先在梅花水道附近做往復運動，移動距離較長，后慢慢前進向東南方向遷移.而中層與底層的粒子由于流速較慢，受潮流場往復作用明顯，粒子在梅花水道向東南方向海域往復前進；落憩時刻釋放的粒子受潮流往復流作用更為顯著.2個排污口釋放的粒子由表及底向外海遷移的運動軌跡逐漸復雜，往復運動越來越強烈.

統計在表、中、底3個水層釋放10個粒子的遷移時間(t)與軌跡長度(L)，結果如表3所示.

表3 粒子的遷移時間和軌跡長度

由表3可知：表層的粒子比中層、底層的粒子更快向閩江河口外遷移.粒子的軌跡長度表明了粒子的滯留能力，軌跡越長，表明粒子更易滯留；反之，表明粒子容易向外遷移.落憩時刻釋放粒子的運動軌跡長度均明顯大于漲憩時刻，表明漲憩時刻釋放的粒子更容易向閩江河口外遷移，不易滯留；漲憩時刻釋放粒子運動軌跡長度由表層及底層逐漸變長，表明表層污染物更易向閩江河口外遷移.

圖9為閩江河口粒子群位置示意圖.由圖9可知：在徑流與潮流作用下，粒子群在漲憩與落憩時刻釋放時，都呈現向外海輸移擴散的態勢.由于北支長門水道表層余流流速較快，在1.5 d后，受河口開口方向影響，北支長門水道流速比南支梅花水道快，漲憩時刻釋放的一部分粒子群已經運移到金牌門附近，而梅花水道流速較慢，較少粒子運移到猴嶼鄉右側水道，個別粒子經由梅花水道運移到瑯岐島東側；而落憩時刻釋放的粒子群遷移速度比漲憩時刻釋放的粒子群更快，長門水道裹挾大部分的粒子到達川石水道，再分別從烏豬水道、熨斗水道遷移，而漲憩時刻釋放的粒子群在2 d后才發生這一過程.3 d后，絕大部分粒子均離開長門水道，一部分粒子從烏豬水道、熨斗水道遷移進入東北部海域，而進入東北部海域的粒子位移距離變短，容易滯留較長時間，向東北方向遷移緩慢，這是由于該海域余流流速較弱，水體交換能力弱，污染物不易向外海遷移，需要較長時間通過海水交換凈化，這與節2.2分析余流場的結果一致；從川石水道與梅花水道遷移出的粒子均向東與東南方向移動，由于此處余流流速較大，使其運動范圍擴大，向外海遷移擴散，5 d后，此海域的大部分粒子均離開閩江河口.大概6 d后，閩江河口的全部粒子均遷移至外海.由此可見，烏豬水道與熨斗水道沿岸即川石水道以北沿岸不是較理想的排污口，排污口應設于梅花水道或瑯岐島東側，便于污染物向外海輸移擴散.

(a) 表層漲憩時刻釋放

在閩江河口自梅花鎮向寨洋鎮方向設置一處斷面(如圖9紅色線段所示)，分別統計在2個時刻釋放粒子群后，半數粒子群通過斷面的時間.經統計，在漲憩時刻釋放粒子群，大概3 d后，半數粒子群通過斷面；在落憩時刻釋放粒子群，大概2.5 d后，半數粒子群通過斷面.結果表明，落憩時刻釋放粒子群的遷移速度更快，更利于污染物向閩江河口外運移.

統計漲憩與落憩時刻，在5條水道釋放的粒子數量(N)，結果如表4所示.由表4可知：與南支梅花水道相比，北支長門水道輸移污染物能力更強，能夠裹挾較多粒子向外海遷移擴散；在不同時刻釋放粒子群，南北支水道表現的輸移污染物能力相差不大，但在漲憩時刻釋放的粒子群比落憩時刻更容易向梅花水道輸運遷移.表現差距較明顯的是熨斗水道，相比漲憩時刻釋放的粒子群，在落憩時刻釋放的粒子群更容易進入熨斗水道.相比漲憩時刻，雖然落憩時刻釋放的粒子群能更快地向閩江河口外遷移，但落憩時刻釋放的粒子群更易向東北部海域遷移，滯留時間更長.綜上所述，漲憩時刻是較理想的排污時間，該時刻污染物易于向外海進行遷移擴散，不易滯留于東北部海域.

3 結論

基于ECOMSED模型和干濕網格判斷法，考慮漫灘效應，建立閩江河口水動力數值模型，研究閩江河口潮流和余流結構的三維特征，并在此基礎上建立粒子追蹤模型，模擬閩江河口內粒子運動軌跡，得出以下4點主要結論.

1) 閩江河口潮流運動以往復流為主，河道內漲落潮流方向與岸線走向基本一致，航道深槽通道流速大于兩側.潮流受閩江徑流、地形水深的影響較為明顯，落潮流速明顯大于漲潮流速.北支水道的漲落潮流速均大于南支水道.落急時流速呈現明顯的垂向分層現象，流速隨水深的增大而逐漸減小，流速減小的幅度較大.

2) 根據余流場模擬結果可知，北支水道余流強度大于南支水道，表明北支水道的水體輸移和交換能力強于南支水道；垂向上，表層余流流速大于底層余流流速.閩江河口的表、中層余流場以落潮余流為主，而底層余流場則以漲潮余流為主.

3) 粒子運移軌跡模擬結果表明，位于梅花水道的2個排污口釋放的粒子在水體表層、中層、底層的運動軌跡有些不同，在漲憩時刻排放的污染物較不易滯留；表層排放的污染物比底層更容易向閩江河口外遷移.

4) 川石水道以北沿岸不是較理想的排污口位置，此處釋放的粒子容易滯留在東北部海域；排污口應設于梅花水道或川石水道以南沿岸，此處潮、余流較為強烈，水體交換能力較強，便于污染物向外海輸移擴散.通過統計粒子群在各個水道的數量，得出漲憩時刻是較理想的排污時間，污染物易于向外海進行遷移擴散，不易滯留于東北部海域.由此可為排污口選址、陸源污染治理方案提供科學依據.