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

（1.自然資源部第二海洋研究所，浙江杭州 310012；2.自然資源部海洋空間資源管理技術(shù)重點(diǎn)實(shí)驗(yàn)室，浙江杭州 310012）

0 引言

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

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

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

1 研究區(qū)域概述

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

2 研究方法與模型構(gòu)建

2.1 研究方法

基于構(gòu)建的水動力水質(zhì)模型，在研究區(qū)域上游設(shè)定一點(diǎn)源污染排放口，各種情況下污染物排放時間統(tǒng)一設(shè)置為24 h，通過調(diào)整點(diǎn)源的排放濃度，研究污染物在不同擴(kuò)散時間內(nèi)、不同排放濃度下的擴(kuò)散面積，分析污染物擴(kuò)散面積與排放濃度及擴(kuò)散時間之間的響應(yīng)關(guān)系。另外保持污染物排放總量和排放濃度不變，通過調(diào)整污染物排放持續(xù)時間，研究污染物擴(kuò)散面積與排放持續(xù)時間的關(guān)系。

2.2 二維耦合模型

2.2.1 水動力模塊

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

式中：t為時間，s；x、y為笛卡爾坐標(biāo)系；h為水深，m；η為水位，m；u、v為x、y方向的垂線平均流速，m/s；ρ為水的密度；ρ0水的相對密度；Pa為當(dāng)?shù)卮髿鈮海籪=2Ωsinφ為科式力參數(shù)（Ω為地球自轉(zhuǎn)角速率，φ為地理緯度）；g為重力加速度；Sxx、Sxy、Syx、Syy為輻射應(yīng)力分量；S為源匯項(xiàng)；Us、Vs為源匯項(xiàng)水流流速，m/s；τsx、τsy為表面風(fēng)應(yīng)力；τbx、τby為底摩擦應(yīng)力。

2.2.2 對流擴(kuò)散模塊

對流擴(kuò)散模塊基本方程：

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

2.3 模型構(gòu)建

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

圖1 模型計(jì)算域Fig.1 Model calculation domain

2.3.1 邊界條件

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

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

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

圖2 模型網(wǎng)格劃分Fig.2 Model grid division

2.3.2 模型計(jì)算參數(shù)

（1）糙率：計(jì)算域內(nèi)河網(wǎng)為城市河道，考慮到有疏浚整治工程對其進(jìn)行維護(hù)，并結(jié)合試算結(jié)果，河床糙率取0.02。

（2）干濕邊界：為避免模型中干濕交替區(qū)域出現(xiàn)非穩(wěn)流態(tài)，保證模型運(yùn)算穩(wěn)定［29］，干濕水深取模型推薦值，干水深0.005，淹沒水深0.05 m，濕水深0.1 m。

（3）擴(kuò)散系數(shù)：擴(kuò)散系數(shù)取1 m2/s。

（4）降解系數(shù)：通過模型試算結(jié)果以及計(jì)算域內(nèi)實(shí)測水質(zhì)數(shù)據(jù)，最終確定計(jì)算域內(nèi)NH3-N 的降解系數(shù)為0.03/d，BOD5的降解系數(shù)為0.05/d［30］。

3 結(jié)果與分析

3.1 水動力模型驗(yàn)證結(jié)果

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

表2 水位驗(yàn)證誤差計(jì)算表Tab.2 Calculation table of water level verification error

圖3 模型計(jì)算域內(nèi)水位驗(yàn)證Fig.3 Verification of water level in model calculation domain

3.2 水質(zhì)模型驗(yàn)證結(jié)果

選取主河道上僅有的一個水質(zhì)監(jiān)測站點(diǎn)（萬程橋）的數(shù)據(jù)進(jìn)行驗(yàn)證，實(shí)測水質(zhì)數(shù)據(jù)每月一個，監(jiān)測時間為2019 年，選取NH3-N、BOD5兩種污染物對水質(zhì)模型進(jìn)行驗(yàn)證，結(jié)果如圖4所示。

圖4 模型計(jì)算域內(nèi)水質(zhì)驗(yàn)證（萬程橋）Fig.4 Verification of water quality in model calculation domain（Wancheng Bridge）

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

表3 水質(zhì)驗(yàn)證誤差計(jì)算表%Tab.3 Calculation table of water quality verification error

3.3 污染物擴(kuò)散過程研究

3.3.1 污染物擴(kuò)散過程分析

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

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

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

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

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

表4 點(diǎn)源排放濃度對照表Tab.4 Comparison table of point source emission concentration

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

3.3.2 討論分析

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

圖8 不同擴(kuò)散時間內(nèi)污染物擴(kuò)散面積與排放濃度的關(guān)系Fig.8 Relationship between pollutant diffusion area and emission concentration in different diffusion time

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

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

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

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

圖11 相同排放總量下污染物擴(kuò)散面積與排放持續(xù)時間的關(guān)系（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

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

4 結(jié)論

（1）基于MIKE21建立了平原河網(wǎng)水動力水質(zhì)模型，使用實(shí)測水位數(shù)據(jù)和水質(zhì)數(shù)據(jù)對模型進(jìn)行了參數(shù)率定和檢驗(yàn)，水位驗(yàn)證最大相對誤差為3.23%，水質(zhì)驗(yàn)證除個別數(shù)據(jù)點(diǎn)相對誤差較大外，其余數(shù)據(jù)均驗(yàn)證良好，模型參數(shù)較為可靠。

（2）使用率定后的模型對污染物擴(kuò)散面積與事故排放濃度、污染物擴(kuò)散時間之間的響應(yīng)關(guān)系進(jìn)行分析，發(fā)現(xiàn)污染面積與污染物排放濃度呈四次相關(guān)關(guān)系；對污染物擴(kuò)散面積與擴(kuò)散時間之間的關(guān)系進(jìn)行分析，發(fā)現(xiàn)它們之間呈現(xiàn)較為明顯的線性相關(guān)，受到了排放總量的影響，部分?jǐn)?shù)據(jù)在超過一定時間后，污染物擴(kuò)散面積不再隨擴(kuò)散時間增加而增加，即在其他條件一定的情況下，在一定時間內(nèi)，污染物擴(kuò)散面積與擴(kuò)散時間呈較為明顯的線性相關(guān)。

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