閩江感潮河段水動(dòng)力及污染物NH3-N特性

1、閩江感潮河段水動(dòng)力及污染物NH3-N特性閩江是福建省最大的河流，徑流量次于長(zhǎng)江、珠江，位居全國(guó)第3；流域面積60 992 km2，約占福建全省面積的一半，流域范圍內(nèi)人口和工農(nóng)業(yè)區(qū)眾多，其河口地區(qū)又是福建省會(huì)城市福州所在地，交通和經(jīng)濟(jì)發(fā)達(dá)。強(qiáng)烈的人類活動(dòng)影響通過閩江干流向下游河口三角洲和外海傳遞，引起諸如富營(yíng)養(yǎng)化、赤潮、生物多樣性減少、海岸侵蝕、海水入侵、河床下切和水位不斷下降等一系列嚴(yán)重的環(huán)境和生態(tài)問題。因此，開展閩江下游感潮河段的水動(dòng)力學(xué)特性研究，對(duì)閩江河道、河口及近海生態(tài)環(huán)境的保護(hù)、合理開發(fā)閩江河道、河口及近海資源具有重要的應(yīng)用價(jià)值，并可為提高河道、河口及近海環(huán)境生態(tài)變化預(yù)測(cè)能力和海岸帶開

2、發(fā)與綜合管理水平提供理論基礎(chǔ)和科學(xué)依據(jù)，實(shí)現(xiàn)河口及其近海地區(qū)社會(huì)經(jīng)濟(jì)的可持續(xù)發(fā)展。河口是由河流向海洋過渡的區(qū)域，河流與海洋動(dòng)力條件不盡相同，水質(zhì)情況亦不一樣。影響河口地區(qū)鹽水入侵的動(dòng)力因素主要有上游徑流、外海潮汐、風(fēng)和波浪。早在2022年，鄭慧瓊等1調(diào)查統(tǒng)計(jì)，福州境內(nèi)通過各種途徑將廢水排入閩江的工業(yè)企業(yè)共450家，以閩江干流流域統(tǒng)計(jì)，2022年閩江福州段工業(yè)企業(yè)污染源廢水排放總量為1 868.05萬t，COD排放總量為3 326.61 t，NH3N排放總量為470.20 t。潘東曦2和武晶3基于20222022年實(shí)測(cè)河道地形資料，分析南北港河道橫向、縱向及平面的演變特征；研究成果均表明：202

3、22022年，南北港河道總體呈現(xiàn)刷深趨勢(shì)，其中20222022年下切速度較緩，20222022年、20222022年下切速度較快。倪冠韜4利用營(yíng)養(yǎng)指數(shù)法、潛在性營(yíng)養(yǎng)化評(píng)價(jià)法、有機(jī)質(zhì)指數(shù)清臺(tái)和模糊數(shù)學(xué)法的評(píng)價(jià)結(jié)果，認(rèn)為：閩江研究區(qū)域富營(yíng)養(yǎng)化現(xiàn)象較為嚴(yán)重?；陂}江下游感潮河段地形不斷調(diào)整，以及水質(zhì)不容樂觀的情勢(shì)下，對(duì)其水動(dòng)力及污染物NH3N特性進(jìn)行研究。國(guó)內(nèi)外研究較多涉及河流水動(dòng)力、感潮河段污染物輸移。1986年，Van 等5采用水深平均模型模擬了北海南部營(yíng)養(yǎng)物（N、P）和重金屬（Cd，Hg，Cu，Pb，Zn，Cr）等物質(zhì)的輸移。李有為等6以汊道眾多的揚(yáng)中場(chǎng)河段為例，建立了感潮河段二維淺水非恒定潮

4、流數(shù)學(xué)模型，提出：潮流界附近的揚(yáng)中河段漲潮動(dòng)力不強(qiáng)，落潮期則持續(xù)時(shí)間較長(zhǎng)；汊道分流比漲潮期趨向于無序狀態(tài)且穩(wěn)定在落急時(shí)刻附近。潘紅波7針對(duì)長(zhǎng)江南京段建立平面二維水動(dòng)力模型，分析其水流運(yùn)動(dòng)特征，基于拉格朗日法的隨機(jī)游走模型，模擬污染物的對(duì)流輸移及擴(kuò)散輸運(yùn)，分析了典型工況條件下污染物的回落規(guī)律與停留時(shí)間。朱紅偉等8根據(jù)實(shí)驗(yàn)結(jié)果，分別針對(duì)靜態(tài)水流和動(dòng)態(tài)水流條件，利用菲克第1定律和第2定律，推導(dǎo)了感潮河段反映底泥污染物釋放的公式。以淺水緩坡海灘為研究對(duì)象，Sun等9模擬了波浪對(duì)污染物輸移的影響，其結(jié)果與實(shí)測(cè)數(shù)據(jù)吻合良好。基于隨機(jī)流動(dòng)粒子法，Tompson 等10建立了多孔介質(zhì)中溶質(zhì)輸移數(shù)學(xué)模型，認(rèn)為：

5、該隨機(jī)流動(dòng)粒子法可較好地模擬飽和或非飽和介質(zhì)中的對(duì)流擴(kuò)散過程。很多學(xué)者通過修正方程11、模擬流線12、實(shí)驗(yàn)觀測(cè)13、化學(xué)反應(yīng)14、考慮吸附動(dòng)力學(xué)15和粗糙邊界16等深入研究了溶質(zhì)輸移。針對(duì)感潮河段水動(dòng)力與污染物擴(kuò)散特性兩方面的研究較少，該研究分析了閩江感潮河段水動(dòng)力與污染物NH3N擴(kuò)散時(shí)空分布的特性，以期為閩江感潮河段污染物控制提供理論基礎(chǔ)和科學(xué)依據(jù)。1 研究方法以水口電站壩下至入海口為研究對(duì)象，建立水動(dòng)力-水質(zhì)耦合數(shù)學(xué)模型。因河道不同高程物質(zhì)組成不同，結(jié)合分層率定糙率驗(yàn)證水動(dòng)力要素特性的研究成果17，以2022年特征斷面處污染物濃度，率定污染物NH3N擴(kuò)散衰減系數(shù)。在2022年污染物NH3N

6、排放規(guī)劃的邊界條件下，研究閩江下游感潮河段污染物濃度分布情況。1.1模型基本情況閩江下游感潮河段水動(dòng)力水質(zhì)模型基本方程包括水動(dòng)力和水質(zhì)控制方程，其中水動(dòng)力模塊主要采用動(dòng)量方程和連續(xù)性方程，水質(zhì)模擬主要采用質(zhì)量守恒控制方程，建立二維數(shù)學(xué)模型。1.1.1數(shù)學(xué)模型水流連續(xù)方程：ht+hux+hvy=hS（1）x、y方向水流運(yùn)動(dòng)方程：hut+huux+huvy=fvhghxh0paxgh220 x+sx0bx010(Sxxx+Sxyy)+x(hTxx)+y(hTxy)+husS（2）hvt+huvx+hvvy=fuhghyh0paygh220y+sy0by010(Syxx+Syyy)+x(hTxy)+

7、y(hTyy)+hvsS（3）式中：t為時(shí)間；x、y為笛卡爾坐標(biāo)系坐標(biāo)；為水位；h為靜止水深；pa為當(dāng)?shù)卮髿鈮簭?qiáng)；u、v分別為x、y方向上的速度分量；f為科氏力系數(shù)，f=2sin，為地球自轉(zhuǎn)角速度，為當(dāng)?shù)鼐暥?；g為重力加速度；、0為水的密度和水的相對(duì)密度；Sxx、Sxy、Syx、Syy分別為輻射應(yīng)力分量；S為源項(xiàng)；Tij為剪應(yīng)力；bx、by為河床切應(yīng)力；sx、sy為風(fēng)剪應(yīng)力；（us，vs）為源項(xiàng)水流流速。字母上帶橫杠的是平均值，例如，u，v為沿水深平均的流速。紊動(dòng)擴(kuò)散系數(shù)反映不同水動(dòng)力條件下污染物的擴(kuò)散現(xiàn)象，對(duì)于有降解過程的污染物來說，河道污染物對(duì)流擴(kuò)散降解的濃度變化規(guī)律的關(guān)系式：(hC)t+

8、(huC)x+(hvC)y=x(hExCx)+y(hEyCy)hKC+q0C0（4）式中：C為水中污染物的濃度；Ex、Ey為橫向和縱向紊動(dòng)擴(kuò)散系數(shù)；K為衰減系數(shù)；q0為源匯處單位面積上的流量；C0為初始斷面的污染物濃度。降解過程滿足一級(jí)反應(yīng)方程式：dC/dt=KC。1.1.2數(shù)值計(jì)算方法直角坐標(biāo)系下，水流運(yùn)動(dòng)的控制方程可用如下通用形式表示：(H)t+(uH)x+(vH)y=x(Hx)+y(Hy)+S（5）式中：為通用變量；為廣義擴(kuò)散系數(shù)；S為源項(xiàng)。以三角形與四邊形單元所構(gòu)成的混合網(wǎng)格為控制體，待求變量存儲(chǔ)于控制體中心。采用有限體積法對(duì)控制方程進(jìn)行離散，用基于同位網(wǎng)格的壓力耦合方程組的半隱式方法

9、（semi-implicit method for pressure linked equations，SIMPLE）處理水流運(yùn)動(dòng)方程中水深和速度的耦合關(guān)系。離散后的代數(shù)方程組可以寫成如下形式：APP=j=1NeAEjEj+b0（6）式中：AP、AEj與b0為系數(shù)；Ne表示控制體的邊數(shù)，對(duì)不同的計(jì)算網(wǎng)格Ne的取值不同；P、Ej為通量。用Gauss迭代法求解離散后的代數(shù)方程組，具體步驟如下：1）給全場(chǎng)賦以初始的猜測(cè)水位。2）計(jì)算動(dòng)量方程系數(shù)，求解動(dòng)量方程。3）計(jì)算水位修正方程的系數(shù)，求解水位修正值，更新水位和流速。4）根據(jù)單元?dú)堄噘|(zhì)量流量和全場(chǎng)殘余質(zhì)量流量判斷是否收斂。在工程計(jì)算中，一般來說當(dāng)單

10、元質(zhì)量流量達(dá)到全局質(zhì)量流量的0.01%，全場(chǎng)殘余質(zhì)量流量達(dá)到進(jìn)口流量的0.5%時(shí)即可認(rèn)為迭代收斂。1.1.3研究區(qū)域考慮局部區(qū)域加密，生成三角形網(wǎng)格。網(wǎng)格大小為5100 m，網(wǎng)格總數(shù)為11 883個(gè)。采用2022年10月水下地形（1500）數(shù)據(jù)，采用最臨近插值點(diǎn)法進(jìn)行插值，得到研究區(qū)域每個(gè)網(wǎng)格單元高程，如圖1所示。圖1研究區(qū)域（單位：m）Fig.1Study area（unit: m）1.2邊界條件因重點(diǎn)研究污染物的擴(kuò)散范圍，故主要關(guān)注多年平均流量和最小下泄流量等不利條件下的水動(dòng)力和污染物濃度時(shí)空分布情況。上邊界為水口電站的流量邊界，下邊界為外海潮位邊界，具體工況條件如表1所示。曼寧系數(shù)M采用

11、分層糙率取值。簡(jiǎn)言之，因不同高程處河道植物特性不同，其阻水不同，產(chǎn)生的床面糙率亦不同，故將河道分為若干層，糙率分別取值。表1工況條件Table 1Working conditions工況流量/（m3s-1）外海實(shí)時(shí)潮位時(shí)間11 741.22022/1/202022/1/2421 741.22022/9/272022/10/13308.02022/1/202022/1/244308.02022/9/272022/10/11.3模型驗(yàn)證引用福州市環(huán)境科學(xué)研究院編制水口水庫(kù)下泄流量、北港分流比河污染源對(duì)閩江下游溶解氧影響的物理模型研究提供的2022年閩江氨氮入河污染綜合普查數(shù)據(jù)作為2022年驗(yàn)證數(shù)據(jù)

12、；引用福州市環(huán)境總體規(guī)劃規(guī)劃年2022年氨氮污染源作為現(xiàn)況污染源數(shù)據(jù)。國(guó)控及省控相關(guān)監(jiān)測(cè)點(diǎn)位位置如圖2所示。圖2閩江下游省控?cái)嗝嫖恢梅植紙DFig.2Location distribution of provincial control section in the downstream of Minjiang River根據(jù)各省控?cái)嗝妫ㄈ鐖D2所示）的實(shí)測(cè)污染物濃度數(shù)據(jù)，以實(shí)時(shí)下泄流量（上邊界）和潮位（下邊界）為邊界條件，將水動(dòng)力與水質(zhì)耦合數(shù)值模擬。進(jìn)行污染源的排污口概化時(shí)應(yīng)遵循如下原則：1）當(dāng)工業(yè)企業(yè)排污口污染物排放流量較大（超過單元總量的10 %），必須作為獨(dú)立的概化排污口處理。2）其他排污

13、口若距離較近，可把多個(gè)排污口簡(jiǎn)化成集中的排污口。3）距離較遠(yuǎn)并且排污量均比較小的分散排污口，可概化為非點(diǎn)源入河。4）大型的污水處理廠需作為概化排污口考慮。5）城市人口聚集地需概化排污口。從而，得出概化排污口分布，如圖3所示；概化排污口濃度如表2所示。圖3概化排污口位置分布圖Fig.3Distribution of generalized sewage outlet表2現(xiàn)狀2022年與規(guī)劃年2022年排污概化成果Table 2Result of generalized sewage outlet in 2022 (current situation) and 2022 (planning yea

14、r)排污口編號(hào)污水廠編號(hào)面源范圍氨氮年排放量/m32022年2022年1W10閩清縣、梅城鎮(zhèn)、梅溪鎮(zhèn)、白樟鎮(zhèn)、云龍鄉(xiāng)469.89413.382-鴻尾鄉(xiāng)88.1277.523-白沙鎮(zhèn)88.5177.874W7竹崎鄉(xiāng)、甘蔗街道220.69194.155-上街鎮(zhèn)、荊溪鎮(zhèn)340.88299.896W1建新鎮(zhèn)、鼓樓區(qū)、臺(tái)江區(qū)部分245.47215.957W2鼓樓區(qū)、晉安區(qū)、臺(tái)江區(qū)、蓋山鎮(zhèn)部分2 118.81 864.018W3、W6鼓樓區(qū)、晉安區(qū)、臺(tái)江區(qū)、蓋山鎮(zhèn)部分1 788.131 573.109W4城門鎮(zhèn)、羅洲鎮(zhèn)部分226.58199.3310W9-263.96232.2211-南通鎮(zhèn)、南嶼鎮(zhèn)278

15、.06244.6212W8尚干鎮(zhèn)、祥謙鎮(zhèn)、青口鎮(zhèn)、羅洲鎮(zhèn)部分518.03455.7413W5、W11羅興街道、首占鎮(zhèn)、吳航街道、航城街道、營(yíng)前街道1 140.231 003.11總計(jì)-7 787.366 850.91污染物NH3N平均濃度對(duì)比成果如圖4所示。由于魁岐站點(diǎn)靠近福州主城區(qū)，內(nèi)河排污量較大，尤其是光明港內(nèi)河污染較大且較復(fù)雜，且隨著漲退潮的波動(dòng)，水質(zhì)波動(dòng)范圍較大，計(jì)算最大值為0.6 mg/L，最小值為0.26 mg/L。此外，由于人工監(jiān)測(cè)取樣時(shí)刻受附近污染源的影響，存在較大的波動(dòng)性，魁岐站點(diǎn)的月平均值的計(jì)算結(jié)果跟單次取樣存有一定的差距?？傮w看來，驗(yàn)證成果與實(shí)際情況基本相符，可進(jìn)行深入研

16、究。圖4驗(yàn)證成果對(duì)比圖Fig.4Comparison of verification results2 計(jì)算成果分析2.1污染物NH3-N濃度最大值分布以閩江下游河潮河段福州段為研究對(duì)象，根據(jù)非恒定流條件下污染物濃度分布成果，從中繪制了現(xiàn)狀排污條件下研究區(qū)域內(nèi)污染物NH3N濃度最大值分布圖，如圖58所示。經(jīng)數(shù)據(jù)分析，得出：圖5污染物NH3N濃度最大值分布圖（工況1）（單位：km）Fig.5Maximum concentration distribution of pollutant NH3N (in the first working condition) (unit: km)圖6污染物NH3

17、N濃度最大值分布圖（工況2）（單位：km）Fig.6Maximum concentration distribution of pollutant NH3N (in the second working condition) (unit: km)圖7污染物NH3N濃度最大值分布圖（工況3）（單位：km）Fig.7Maximum concentration distribution of pollutant NH3N (in the third working condition)(unit: km)圖8污染物NH3N濃度最大分布圖（工況4）（單位：km）Fig.8Maximum concent

18、ration distribution of pollutant NH3N (in the fourth working condition)(unit: km)1）排污時(shí)空分布特點(diǎn)。閩江北港污染物排污口集中，排放量較大，并且北港排污量是南港排污量的3.9倍。2）河道束窄，限制分流。在南北港分汊口處，北港平均河寬僅南港的1/4.99。3）南港沖淡作用明顯。據(jù)福建省水文局分析，10年一遇洪水頻率條件下，北港過流僅是南港的1/2.4。2.2NH3-N濃度變化的面積對(duì)時(shí)間的響應(yīng)本小節(jié)以福州段（如圖9所示）的污染物各濃度所占的面積為研究對(duì)象，相應(yīng)濃度面積百分比隨時(shí)間的變化如圖1013所示。圖9研究區(qū)域

19、示意圖Fig.9Schematic diagram of study area圖10相應(yīng)濃度面積百分比的變化（工況1）Fig.10The area percentage change of corresponding concentration（in the first working condition）圖11相應(yīng)濃度面積百分比的變化（工況2）Fig.11The area percentage change of corresponding concentration（in the second working condition）圖12相應(yīng)濃度面積百分比的變化（工況3）Fig.12The

20、area percentage change of corresponding concentration（in the third working condition）圖13相應(yīng)濃度面積百分比的變化（工況4）Fig.13The area percentage change of corresponding concentration（in the fourth working condition）從圖1013和數(shù)據(jù)分析，可知：1）總體上，污染物濃度越小，相應(yīng)的濃度面積百分比越大。2）4個(gè)工況變化趨勢(shì)一致；隨著時(shí)間的推移，感潮河段污染物濃度是一個(gè)動(dòng)態(tài)變化的過程。3）對(duì)于4個(gè)工況，各自的2個(gè)峰值

21、之間的時(shí)差為12 h，與潮位峰值時(shí)差幾乎相等。4）在4個(gè)工況條件下，研究區(qū)域NH3N濃度變化平均值分別為：-0.015 4、-0.014 6、-0.028 9和-0.027 6 mg/L。由此可以看出：1）4個(gè)工況的污染沖淡程度，分別是：工況1工況3，工況2工況4。2）在不同流量和同潮周期條件下，流量對(duì)污染物濃度的影響非常大。3）從工況1與工況2計(jì)算成果對(duì)比可以看出：潮水有一定的沖淡效果，但不是十分顯著。4）在4個(gè)工況條件下，研究區(qū)域NH3N濃度減少0.03 mg/L（含）以上的面積百分比的平均值分別為：1.18%、1.15%、10.08%和9.91%；對(duì)于小流量（308 m3/s）而言，20

22、22年后污染物NH3N濃度變化更大。因此，再次說明了潮水有一定的沖淡作用，但并不是非常顯著，比如：工況1與工況2僅差0.03 %。2.3污染物面積變化對(duì)潮位的響應(yīng)圖1417給出了污染物面積變化對(duì)潮位的響應(yīng)。從圖1417可以看出：圖14污染物面積變化對(duì)潮位的響應(yīng)（工況1）Fig.14Response of the change of pollutant area to tide level（in the first working condition）圖15污染物面積變化對(duì)潮位的響應(yīng)（工況2）Fig.15Response of the change of pollutant area to tide level（in the second working condition）圖16污染物面積變化對(duì)潮位的響應(yīng)（工況3）Fig.16Response of the change of pollutant area to tide level（in the third working condition）圖17污染物面積變化對(duì)潮位的響應(yīng)（工況4）Fig.17Response of the change of pollutant area to tide level（in the fourth working condition）1）隨著潮位變化，污染物濃度減少0.03 mg/L的面