基于水動力水質模型的平原河網排污模擬分析

袁行知，許雪峰，2，俞亮亮，楊萬康，壽鹿

（1.自然資源部第二海洋研究所，浙江杭州 310012；2.自然資源部海洋空間資源管理技術重點實驗室，浙江杭州 310012）

0 引言

平原河網地區水系縱橫交錯，在防洪排澇、生態健康和文化景觀等方面發揮著重要作用［1］。河網地區經濟文化相對發達，但同時也是水質污染高發區［2］，受潮汐影響，河流流動呈現往復性，污染物在河流內停留時間較長，遷移擴散過程相對較慢，不利于污染物的降解，水質變化情況復雜［3-5］，河網中修建的水閘工程又影響了水體交換更新的速度，進一步加劇了該過程［6］，汛期水位波動對水質也有一定影響［7］。隨著我國城鎮化步伐的加快，人口大量向城市涌入，這一方面促進當地經濟的發展，另一方面也使得城鎮工業污染日益嚴重，造成污染聚集［8］，這其中包括大量未經處理就直接排入河道的工業副產污染物，其中氨氮會影響魚類攝食與生長，含量高時甚至會導致魚類死亡，人類長期飲用氨氮超標水體易引發癌癥；BOD5能夠衡量水體受有機物污染的程度，是評價水質的重要指標。為了在事故排放后盡快確定污染物的污染范圍，并進一步采取相關措施降低損失，對排放入河道污染物遷移擴散規律的研究以及污染物影響區域的分析研究對于水質改善至關重要［9，10］。

對于平原河網地區點源污染的過程和影響一般使用數值模擬的方法進行研究。曹霞等［11］使用二維EFDC模型建立了賈魯河鄭州段的污染源與水質關系響應模型，為其他河網水質模型的建立提供了指導；戴君等［12］以松花江哈爾濱段作為研究對象，建立了水動力-水質模型，研究了多情景污染負荷下河流水質的響應關系；李丹等［13］建立了珠江流域河網水動力水質模型，研究了污染排放與水質響應關系；包子云等［14］建立了太湖流域一維河網水動力模型，研究了入河排污口設置對環境的影響；趙士文等［15］以梁溪區河網為研究對象，對河網污染負荷進行了估算分析。對水質污染特征的分析以及對水質模擬結果量化表示是常用的研究方法［16，17］，常用的平面二維數值模型有MIKE21、Deflt3D、EFDC等［18-20］，其中MIKE21 由于其對各種地表水環境的適用性較強而被廣泛應用于工程問題和學術研究中［21］。在水庫總氮預測研究［22］、河道水流特性分析［23］以及污染物擴散數值分析［24］等問題中均有應用。但總體來看，對污染源濃度與污染物擴散面積以及污染源排放時間與污染物擴散面積之間的響應關系研究較少。

嘉興地區河網總體隸屬于太湖流域水系［25］，區域內地勢較為平坦，湖泊、河流遍布，在水源供應、航運和文化旅游等方面發揮著重要作用。但由于區域內工業園區眾多，易發生污染物事故排放事件，給環境保護造成巨大壓力［26，27］。本文通過結合實測數據，建立嘉興地區河網水動力-水質模型，對點源污染的排放過程進行模擬分析，研究了污染物擴散面積與污染源排放模式及時間之間的響應關系。

1 研究區域概述

研究區域位于嘉興市主城區東南側，地處北緯30°34'18″～30°46'9″、東經120°52'17″～121°14'45″之間，區域內河道眾多，主要支流有上海塘、海鹽塘兩條河流，總體上屬于太湖流域。研究區域主要包括平湖塘流域和黃姑塘流域兩部分，平湖塘位于南湖區和平湖市之間，黃姑塘位于平湖市下游。由于流域內水系復雜，主河道兩側工業園區密布，事故情況下可能出現污染物泄露入河的情況，且研究區域橫跨嘉興市南湖區、平湖市兩個主城區，流域內人口密度大，一旦出現較為嚴重污染物泄漏事件，將對河網水環境、生態系統健康造成較大影響。

2 研究方法與模型構建

2.1 研究方法

基于構建的水動力水質模型，在研究區域上游設定一點源污染排放口，各種情況下污染物排放時間統一設置為24 h，通過調整點源的排放濃度，研究污染物在不同擴散時間內、不同排放濃度下的擴散面積，分析污染物擴散面積與排放濃度及擴散時間之間的響應關系。另外保持污染物排放總量和排放濃度不變，通過調整污染物排放持續時間，研究污染物擴散面積與排放持續時間的關系。

2.2 二維耦合模型

2.2.1 水動力模塊

水動力模塊（HD）是MIKE21 的基本模塊，包括平面二維連續方程和動量方程。

式中：t為時間，s；x、y為笛卡爾坐標系；h為水深，m；η為水位，m；u、v為x、y方向的垂線平均流速，m/s；ρ為水的密度；ρ0水的相對密度；Pa為當地大氣壓；f=2Ωsinφ為科式力參數（Ω為地球自轉角速率，φ為地理緯度）；g為重力加速度；Sxx、Sxy、Syx、Syy為輻射應力分量；S為源匯項；Us、Vs為源匯項水流流速，m/s；τsx、τsy為表面風應力；τbx、τby為底摩擦應力。

2.2.2 對流擴散模塊

對流擴散模塊基本方程：

式中：c為污染物濃度；Dx、Dy為x、y方向上的擴散系數，m2/s；F為線性彌散系數，s-1；S為源匯項。

2.3 模型構建

本模型以平湖塘、黃姑塘為主河道，向下游延伸至入海閘門附近，將主河道兩側3 km內的主要河流納入計算域。河網區域河流密度大，河寬較窄，因此使用四邊形網格劃分更適宜；網格采用二維模型進行構建，網格劃分過程中，在河道橫向上盡量保證有兩個網格，但對于部分過窄的支流河段，在橫向上僅劃分一個網格，否則會嚴重影響計算效率。模型共劃分為26 222 個網格節點，16 956 個網格單元；汊道相交處采用三角形網格連接，模型計算域及網格如圖1、2所示。

圖1 模型計算域Fig.1 Model calculation domain

2.3.1 邊界條件

（1）水動力邊界（HD）：模型內共設置六條開邊界，平湖塘、鳳橋、陳良及步云開邊界采用2019 年水位實測值，每小時記錄一個水位；上海塘及海鹽塘開邊界水位采用預報數據。

（2）水質邊界（AD）：通過對2019 年計算域內測站實測數據求平均值，并結合實際情況對數據進行修正，確定除平湖塘上游邊界外，其他邊界及初始濃度為：NH3-N 為0.1 mg/L、BOD5為2 mg/L，初始污染物濃度也依據此數據進行賦值。平湖塘邊界污染物濃度依據實測數據進行賦值如表1所示。

表1 平湖塘邊界2019年各月份入河污染物濃度Tab.1 Concentration of pollutants entering the river at Pinghutang boundary in each month in 2019

圖2 模型網格劃分Fig.2 Model grid division

2.3.2 模型計算參數

（1）糙率：計算域內河網為城市河道，考慮到有疏浚整治工程對其進行維護，并結合試算結果，河床糙率取0.02。

（2）干濕邊界：為避免模型中干濕交替區域出現非穩流態，保證模型運算穩定［29］，干濕水深取模型推薦值，干水深0.005，淹沒水深0.05 m，濕水深0.1 m。

（3）擴散系數：擴散系數取1 m2/s。

（4）降解系數：通過模型試算結果以及計算域內實測水質數據，最終確定計算域內NH3-N 的降解系數為0.03/d，BOD5的降解系數為0.05/d［30］。

3 結果與分析

3.1 水動力模型驗證結果

模型計算域內有新豐、平湖兩個水位測站，位置如圖1 所示。選取2019 年5 月25 日至2019 年6 月6 日之間的實測水位數據進行驗證，結果如圖3 所示，誤差如表2 所示。兩個站點水位驗證相對誤差均小于4%，表明本模型模擬精度高，模擬結果可靠。

表2 水位驗證誤差計算表Tab.2 Calculation table of water level verification error

圖3 模型計算域內水位驗證Fig.3 Verification of water level in model calculation domain

3.2 水質模型驗證結果

選取主河道上僅有的一個水質監測站點（萬程橋）的數據進行驗證，實測水質數據每月一個，監測時間為2019 年，選取NH3-N、BOD5兩種污染物對水質模型進行驗證，結果如圖4所示。

圖4 模型計算域內水質驗證（萬程橋）Fig.4 Verification of water quality in model calculation domain（Wancheng Bridge）

模型水質驗證結果見圖4，各月份相對誤差計算如表3 所示，對于氨氮，除去5 月份和11 月份的數據外，其他十個月的平均相對誤差為18.4%，推測五月份計算值偏小的原因為計算域中的河段存在臨時的污染物排放，從而導致計算值較實測值偏??；推測11 月份相對誤差較大的原因為實測污染物濃度較小，導致在絕對誤差不大的情況下，其對應的相對誤差較大；BOD5的驗證結果要優于氨氮的驗證結果，平均相對誤差為13.4%；整體而言，表明模型模擬精度較高，可以較為精確地重現污染物輸移過程。

表3 水質驗證誤差計算表%Tab.3 Calculation table of water quality verification error

3.3 污染物擴散過程研究

3.3.1 污染物擴散過程分析

為了研究：①污染物擴散面積（指污染物最大濃度掃過的包絡面積，下同）與其排放濃度之間的關系；②污染物擴散面積與其擴散時間之間的關系；③排放總量一定時，污染物擴散面積與其排放持續時間之間的關系；使用率定好的模型對其進行研究，各開邊界污染物濃度設定為NH3-N 為0.1 mg/L、BOD5為2 mg/L。在研究區域上游設置一個點源排放口，污染物事故排放濃度設置如表4 所示，分析污染物的擴散面積與事故排放濃度、擴散時間之間的響應關系。另外，設置點源排放工況如表5所示，研究在排放總量一定的情況下，保持污染物排放濃度不變，污染物擴散面積與排放持續時間的關系，以污染物開始排放為時間起點，計算72 h 內的污染物擴散面積。以氨氮為例，計算72 h 內污染物的擴散面積，在其他條件相同的情況下，污染物擴散面積隨排放濃度變化的結果如圖5 所示，展示結果為N1、N3、N6、N9 四種具有代表性的工況；以氨氮為例，排放濃度為500 mg/L（N9）時，在其他條件相同的情況下，污染物擴散面積隨擴散時間變化的結果如圖6 所示，選取擴散時間為12、36、48、72 h 四種工況進行展示。以氨氮為例，污染面積統計時間為排放開始后的72 h 內，污染物排放總量一定時，污染物擴散面積隨污染物排放持續時間變化的結果如圖7所示。

圖6 污染物排放后不同時間內的擴散范圍（N9）Fig.6 Diffusion range of pollutants in different time after discharge（N9）

圖7 不同排放持續時間下污染物擴散范圍（排放總量一定）Fig.7 Diffusion range of pollutants under different discharge durations（the total discharge amount is certain）

表5 污染源排放工況設置Tab.5 Setting of pollution source discharge conditions

圖5 不同排放濃度下72 h內污染物擴散范圍Fig.5 Diffusion range of pollutants within 72 hours under different emission concentrations

表4 點源排放濃度對照表Tab.4 Comparison table of point source emission concentration

以氨氮為例，污染物在不同排放濃度、排放后不同時間內含量達到1.5 mg/L（Ⅳ類水標準限值）的污染面積如圖8 所示；在排放濃度為500 mg/L（N9）的情況下，污染物擴散面積隨時間的變化情況如圖9 所示，其中“S0.5”表示計算域內氨氮含量達到0.5 mg/L的水域面積，其他以此類推；在排放總量不變的情況下，排放濃度設置為500 mg/L（N9），污染物擴散面積與排放持續時間的關系數據如圖11所示。

3.3.2 討論分析

（1）模型中氨氮濃度達到1.5 mg/L的污染面積，在不同時間內不同排放濃度下的污染物的擴散面積如圖8所示。用四階曲線（圖中虛線）對數據點進行擬合，數據點與曲線擬合良好，除其中一條曲線對應的相關系數為0.987 7，其余各曲線對應相關系數均在0.99 以上，表明在其他條件相同的情況下，污染物擴散面積與污染物排放濃度的四次方呈正相關。

圖8 不同擴散時間內污染物擴散面積與排放濃度的關系Fig.8 Relationship between pollutant diffusion area and emission concentration in different diffusion time

（2）以氨氮排放濃度為500 mg/L（N9）為例，對污染物排放后的擴散面積與擴散時間的關系進行分析，關系曲線如圖9 所示，其中S0.5 表示計算域內氨氮濃度達到0.5 mg/L 及以上的區域面積，其他以此類推。從圖9中可以看出，不同污染物排放濃度下，污染物擴散面積與時間呈現明顯的線性關系，對圖9中各條曲線用線性關系進行擬合，各組數據的相關系數均在0.99以上；以氨氮排放濃度為250 mg/L（N4）為例，繼續對污染物排放后的擴散面積與時間的關系進行分析，如圖10 所示，可以看出S0.5、S1.0所對應的曲線仍表明，污染物擴散面積與時間有明顯的線性關系，對其進行直線擬合，兩組數據線性相關系數均在0.99 以上。而對于S1.5、S2 所對應的曲線，其后半部分對應污染面積不再隨時間增加而增加，這與常識相符合，即一次事故排放的污染物只能污染一定面積的水域，其污染面積不可能隨時間增加而無限增加。對于N4和N9兩種排放濃度，由于N9在相同時間內排放的污染物是N4 的兩倍，所以S1.5 和S2.0 兩條曲線在圖9 可以隨時間增加一直向上延伸，而在圖10 中則到達一定面積后不再增加。綜上可以得出，在其他條件相同的情況下，在一定的時間內，污染物擴散面積與擴散時間呈現較為明顯的線性相關。

圖9 相同排放濃度下污染物擴散面積與擴散時間的關系（N9）Fig.9 Relationship between diffusion area and diffusion time of pollutants at the same emission concentration（N9）

圖10 相同排放濃度下污染物擴散面積與擴散時間的關系（N4）Fig.10 Relationship between diffusion area and diffusion time of pollutants at the same emission concentration（N4）

（3）以氨氮排放濃度為500 mg/L（N9）為例，研究在污染物排放總量一定的情況下，排放持續時間與污染面積之間的關系，污染面積統計時間為開始排放后的72 h 內。根據圖11 所示，在污染物排放總量不變、排放濃度不變的情況下，污染物的擴散面積與污染物排放持續時間呈現明顯的負相關，即在排放總量不變的情況下，污染物排放速度越慢，其造成的污染面積越小。但是在對氨氮含量大于等于1.5 mg/L（Ⅳ類水標準限值）、2.0 mg/L（Ⅴ類水標準限值）的污染面積進行統計時發現，這兩組數據對應的曲線的后面一段存在一個“斷崖式”下降的過程。以排放持續時間等于32 h 的污染面積為基準值1，計算排放持續時間等于48 h的污染面積，數據如表6所示。

圖11 相同排放總量下污染物擴散面積與排放持續時間的關系（N9）Fig.11 Relationship between pollutant diffusion area and discharge duration under the same total discharge amount（N9）

表6 不同排放時間污染面積對比Tab.6 Comparison of polluted areas at different discharge times

數據表明，污染物排放持續時間從32 h 增加到48 h 的過程中，排放持續時間雖然只在32 h 的時間上增加了一半，但氨氮濃度達到1.5 mg/L 的污染面積卻減少了超過50%，而對于氨氮濃度達到2.0 mg/L 的污染面積，排放48 h 對應的數據僅僅只有排放32 h 對應數據的16%，由此可見，污染物排放后的污染面積除了與污染物排放濃度、污染物擴散時間有關外，還與排放模式有著關聯，尤其是對于污染物濃度含量較高的水域面積的統計，排放模式的影響更是不可忽略。這點應該引起有關監管部門的注意，對于一些有污染物排放需求的工業園區，在其附近設置的水質監測點應適當密集，避免其通過延長排放時間的手段向環境中排入超標污染物，從而躲避相關部門的監察。

4 結論

（1）基于MIKE21建立了平原河網水動力水質模型，使用實測水位數據和水質數據對模型進行了參數率定和檢驗，水位驗證最大相對誤差為3.23%，水質驗證除個別數據點相對誤差較大外，其余數據均驗證良好，模型參數較為可靠。

（2）使用率定后的模型對污染物擴散面積與事故排放濃度、污染物擴散時間之間的響應關系進行分析，發現污染面積與污染物排放濃度呈四次相關關系；對污染物擴散面積與擴散時間之間的關系進行分析，發現它們之間呈現較為明顯的線性相關，受到了排放總量的影響，部分數據在超過一定時間后，污染物擴散面積不再隨擴散時間增加而增加，即在其他條件一定的情況下，在一定時間內，污染物擴散面積與擴散時間呈較為明顯的線性相關。

（3）在污染物排放總量、排放濃度不變的情況下，通過改變污染源的排放模式發現，在保持污染物排放濃度不變的情況下，污染面積隨排放持續時間增加而減少，即單位時間內排放的污染物總量越少，其污染的面積越小。而對于濃度較高的污染面積進行統計時發現，當排放持續時間超過一定值時（本文中為32 h），污染面積會隨排放時間增加而出現“斷崖式”下降。