河流交匯區(qū)污染物擴(kuò)散數(shù)值模擬研究

楊志山，米潭，姚建，張弢

(1.四川大學(xué) 建筑與環(huán)境學(xué)院，四川 成都 610065；2.成都市生態(tài)環(huán)境科研監(jiān)測所，四川 成都 610072)

河流交匯現(xiàn)象普遍存在于自然的河流水系中，交匯河道的邊界岸線、河床地形等具有不斷變化的特點(diǎn)[1]，導(dǎo)致交匯區(qū)污染物的擴(kuò)散規(guī)律難以明晰。不同的交匯流量是影響污染物擴(kuò)散的主要因素[2]，研究不同交匯流量下污染物的擴(kuò)散特征，對工程實(shí)踐中涉及的污水排放方式、排污口位置的選擇等具有重要指導(dǎo)意義。目前，水質(zhì)數(shù)值模型逐漸成為研究水體中污染物擴(kuò)散規(guī)律的主要工具[3]，其中，平面二維自由表面流MIKE 21 FM模型在河流、湖泊、水庫的污染物擴(kuò)散研究中已得到廣泛應(yīng)用。如：張守平等[4]通過MIKE 21 FM模型研究了入河排污口對水體水質(zhì)的影響，結(jié)果表明該水質(zhì)模型能準(zhǔn)確反映污染物在水體中的擴(kuò)散規(guī)律；賈瑞鵬等[5]建立了萬寶湖水動力水質(zhì)二維耦合模型，分析了萬寶湖各類污染物的遷移轉(zhuǎn)化機(jī)制；楊晨等[6]考慮風(fēng)場對污染物擴(kuò)散的影響，模擬了汾河水庫11種污染物的遷移擴(kuò)散規(guī)律。

關(guān)于河流交匯區(qū)不同匯流比情況下污染物擴(kuò)散特征方面的研究，多為建立水槽物理模型進(jìn)行探討。如：魏娟等[7]選用水氣兩相流數(shù)學(xué)模型，通過物理實(shí)驗(yàn)，對河流交匯區(qū)不同流量情況進(jìn)行實(shí)驗(yàn)和數(shù)值模擬研究，結(jié)果表明污染帶的狹長程度受流量比的影響較大；袁航[8]通過建立45°支流斜接干流型的交匯水槽物理模型，探究了交匯河道不同匯流比對污染物混合特性的影響，結(jié)果表明水流交匯區(qū)會產(chǎn)生二次環(huán)流，對污染物混合有很大的影響；HAN S等[9]對室內(nèi)交匯明渠的污染物混合過程進(jìn)行了數(shù)值模擬，發(fā)現(xiàn)交匯口附近污染物的混合速度主要由上游來流流量決定，來流流量較小時(shí)混合速度較快。但是，上述研究中物理實(shí)驗(yàn)?zāi)Ｐ屯鶎⒑拥篮喕癁轫樦焙拥?，難以反映自然河流交匯區(qū)邊界、地形等復(fù)雜多變的特點(diǎn)。

為此，本研究基于平面二維自由表面流MIKE 21 FM模型，充分考慮自然河流交匯區(qū)的特點(diǎn)，選取中國西南中小河流石亭江與其支流鴨子江的交匯區(qū)，研究交匯區(qū)COD、氨氮質(zhì)量濃度在不同匯流比情況下的擴(kuò)散規(guī)律，為交匯區(qū)污染物擴(kuò)散機(jī)理研究提供借鑒。

1 模型的建立與驗(yàn)證

1.1 模型建立

平面二維自由表面流MIKE 21 FM模型的控制方程為基于Boussinesq假定和流體靜壓假定的二維不可壓縮雷諾平均N-S方程。建立的模型包括二維水動力學(xué)模型MIKE 21 HD和水質(zhì)模型MIKE 21 AD。

網(wǎng)格劃分：本研究計(jì)算區(qū)域范圍為東西向2.00 km、南北向1.95 km。采用無結(jié)構(gòu)三角形網(wǎng)格對研究區(qū)域進(jìn)行網(wǎng)格劃分，最后共生成網(wǎng)格1 962個(gè)，計(jì)算節(jié)點(diǎn)1 300個(gè)。在搭建計(jì)算網(wǎng)格的基礎(chǔ)上，利用實(shí)測天然地形，對網(wǎng)格地形進(jìn)行插值，得到計(jì)算區(qū)域的地形模擬圖(圖1)。

圖1 地形模擬圖

邊界設(shè)置：本研究計(jì)算區(qū)域共設(shè)置3個(gè)開邊界：交匯區(qū)上游約1.0 km處石亭江省控?cái)嗝鍸設(shè)為上邊界1；交匯區(qū)上游約800 m處鴨子江國控?cái)嗝鍿設(shè)為上邊界2；交匯區(qū)下游約1.2 km處石亭江出境斷面W設(shè)為下邊界。上邊界采用實(shí)測的流量數(shù)據(jù)作為控制條件，下邊界采用實(shí)測的水位數(shù)據(jù)作為控制條件。

斷面設(shè)置：為方便模型驗(yàn)證及污染物擴(kuò)散特征分析，交匯區(qū)域共設(shè)置6個(gè)斷面：Z1、Z2、Z3共3個(gè)橫向斷面，反映污染物橫向擴(kuò)散特征；C1、C2、C3共3個(gè)縱向斷面，反映污染物縱向擴(kuò)散特征。邊界及斷面設(shè)置如圖2所示。

圖2 邊界及斷面設(shè)置圖

計(jì)算條件設(shè)定：①初始條件：二維水動力模型選取2015年1月1日的水動力學(xué)指標(biāo)實(shí)測資料；水質(zhì)模型選取該日各監(jiān)測斷面COD、氨氮的質(zhì)量濃度實(shí)測資料。②邊界條件：二維水動力模型選取2015年1月1—30日L、S斷面的實(shí)測流量過程及W斷面的實(shí)測水位過程；水質(zhì)模型選取2015年1月1—30日L、S及W斷面COD、氨氮的質(zhì)量濃度實(shí)測資料。

1.2 模型驗(yàn)證

采用2015年1月1—30日Z1、Z2、Z3斷面的流量及COD、氨氮質(zhì)量濃度的模擬值，分別計(jì)算其平均絕對誤差(Mean Absolute Error，MAE)、決定系數(shù)(R2)、均方根誤差(Root Mean Square Error，RMSE)、平均相對誤差(Mean Relative Error，MRE)4種擬合優(yōu)度指標(biāo)來評價(jià)模型的性能。通過查閱文獻(xiàn)與模型調(diào)試，確定石亭江、鴨子江匯流前河床糙率為0.03，匯流后河床糙率為0.05；渦黏系數(shù)采用Smagorinsky公式估算，取值為0.28 m2/s；干水深為0.005 m，淹沒水深為0.05 m，濕水深為0.1 m；風(fēng)阻系數(shù)為0.001；污染物水平擴(kuò)散系數(shù)為0.05 m2/s；COD、氨氮質(zhì)量濃度的衰減系數(shù)分別為0.5、0.1 d-1。

Z1至Z3斷面上的流量、COD和氨氮質(zhì)量濃度的MAE、R2、RMSE、MRE值見表1。由表1可知，各監(jiān)測指標(biāo)在不同斷面處的MAE、R2、RMSE值均比較接近，說明各監(jiān)測指標(biāo)測量值與模擬值的誤差在同一數(shù)量級；各監(jiān)測指標(biāo)的R2值均在Z3斷面處最大，MRE值在Z3斷面處均為最小，表明Z3斷面處測量值與模擬值擬合最好。MRE值越小，模型的精度越高，本文的MRE值均在30%以內(nèi)，說明用實(shí)測值標(biāo)定的MIKE 21 FM模型可以保證精度，滿足模擬需求。

表1 各變量測量值與模擬值比較的指標(biāo)計(jì)算

2 結(jié)果與分析

采用本文建立的MIKE 21 FM模型對研究區(qū)域不同匯流比情況下的COD、氨氮擴(kuò)散特征進(jìn)行數(shù)值模擬，模擬工況設(shè)置見表2，其中Q在3種工況下分別選取石亭江、鴨子江在枯水期、平水期及豐水期時(shí)的平均流量，匯流比R=Q鴨子江/Q石亭江，模擬初始污染物質(zhì)量濃度C選取2015年1月1日L、S斷面實(shí)測的COD、氨氮質(zhì)量濃度。

表2 數(shù)值模擬工況參數(shù)

2.1 污染物濃度場分析

采用水質(zhì)模型MIKE 21 AD對不同工況下的COD、氨氮質(zhì)量濃度進(jìn)行模擬計(jì)算，獲得3種工況下COD、氨氮的質(zhì)量濃度場如圖3所示。根據(jù)《地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)》(GB 3838—2002)[10]中Ⅲ、Ⅳ類水質(zhì)的COD、氨氮質(zhì)量濃度范圍，得到COD、氨氮質(zhì)量濃度的分布面積見表3。

圖3 不同工況下COD、氨氮質(zhì)量濃度場

表3 COD、氨氮質(zhì)量濃度分布面積統(tǒng)計(jì)

由圖3可知，在石亭江和鴨子江交匯前，兩條河流的污染物質(zhì)量濃度均逐漸降低。石亭江COD的質(zhì)量濃度在不同匯流比下均大于22 mg/L，匯流比越大，COD質(zhì)量濃度分布在22～24 mg/L范圍內(nèi)的面積越小。鴨子江COD質(zhì)量濃度變化受匯流比變化的影響較大，當(dāng)R=0.33時(shí)，COD質(zhì)量濃度自上游至下游依次分布在18～20 mg/L、16～18 mg/L的范圍內(nèi)，在交匯口附近COD質(zhì)量濃度低于16 mg/L；當(dāng)R=0.50時(shí)，COD質(zhì)量濃度分布在18～20 mg/L的范圍內(nèi)，在交匯口附近分布在16～18 mg/L的范圍內(nèi)；當(dāng)R=0.67時(shí)，COD質(zhì)量濃度僅分布在18～20 mg/L范圍內(nèi)。由此可見，在兩條河流交匯前，匯流比越大，污染物越不易降解。石亭江、鴨子江氨氮質(zhì)量濃度在交匯之前隨匯流比變化不明顯，均分布在1.4～1.5 mg/L、低于1.0 mg/L的范圍內(nèi)。在兩條河流交匯后，污染物混合情況各有不同，當(dāng)R=0.33時(shí)，COD質(zhì)量濃度自上游至下游依次分布在20～22 mg/L、18～20 mg/L、16～18 mg/L的范圍內(nèi)；當(dāng)R=0.50時(shí)，COD質(zhì)量濃度分布依次為20～22 mg/L、18～20 mg/L的范圍內(nèi)；當(dāng)R=0.67時(shí)，COD質(zhì)量濃度分布依次為22～24 mg/L、20～22 mg/L、18～20 mg/L的范圍內(nèi)?？梢姡谙嗤幕旌虾佣尉嚯x內(nèi)，匯流比越小，污染物混合質(zhì)量濃度也越低。不同匯流比下，在交匯后氨氮質(zhì)量濃度均分布在1.4～1.5 mg/L、1.2～1.4 mg/L、1.0～1.2 mg/L的范圍內(nèi)。由于初始濃度CM>CN，交匯后，受石亭江影響，右岸(鴨子江側(cè))污染物的質(zhì)量濃度沿程呈上升趨勢，河流左岸至右岸污染物的質(zhì)量濃度逐漸降低。

研究區(qū)域所在水功能區(qū)水質(zhì)管理目標(biāo)為Ⅲ類。由表3可知，隨著匯流比的增加，COD、氨氮質(zhì)量濃度超過《地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)》Ⅲ類水(COD質(zhì)量濃度>20 mg/L，氨氮質(zhì)量濃度>1 mg/L)標(biāo)準(zhǔn)部分在研究區(qū)域內(nèi)形成的污染帶面積明顯增加，當(dāng)R=0.33時(shí)，COD形成的污染帶面積為159 844 m2，氨氮形成的污染帶面積為351 155 m2；R=0.50時(shí)，COD形成的污染帶面積為207 650 m2，增加30%，氨氮形成的污染帶面積為366 708 m2，增加4.4%；當(dāng)R=0.67時(shí)，COD形成的污染帶面積為309 728 m2，增加93.8%，氨氮形成的污染帶面積為388 465 m2，增加10.6%?？梢?，匯流比越大，污染帶分布越廣，適當(dāng)降低匯流比，對交匯區(qū)水質(zhì)可起到一定的改善作用。

2.2 匯流區(qū)污染物擴(kuò)散特征分析

圖4為河流交匯后COD、氨氮質(zhì)量濃度的分布情況，取不同橫向、縱向斷面分析其污染物擴(kuò)散特征。

圖4 匯流區(qū)COD、氨氮質(zhì)量濃度場

2.2.1 交匯區(qū)污染物縱向擴(kuò)散特征分析

圖5為不同匯流比下COD、氨氮質(zhì)量濃度在同一縱向斷面上隨不同橫向斷面的變化情況，選取C1、C3斷面進(jìn)行縱向擴(kuò)散分析。

由圖5(a)可知：由于石亭江初始污染物的質(zhì)量濃度較高，COD、氨氮沿C1縱向斷面的質(zhì)量濃度逐漸下降；當(dāng)R=0.33時(shí)，Z1至Z2斷面、Z2至Z3斷面COD的質(zhì)量濃度分別降低14.2%、10.4%，氨氮質(zhì)量濃度分別降低12.3%、6.7%；當(dāng)R=0.50時(shí)，COD質(zhì)量濃度分別降低11.7%、7.4%，氨氮質(zhì)量濃度分別降低11.4%、5.0%；當(dāng)R=0.67時(shí)，COD質(zhì)量濃度分別降低11.3%、3.0%，氨氮質(zhì)量濃度分別降低11.0%、3.4%。由圖5(b)可知：由于鴨子江初始污染物的質(zhì)量濃度較低，受石亭江污染物擴(kuò)散的影響，COD、氨氮沿C3縱向斷面的質(zhì)量濃度逐漸升高；當(dāng)R=0.33時(shí)，Z1至Z2斷面、Z2至Z3斷面COD質(zhì)量濃度分別增加6.5%、2.4%，氨氮質(zhì)量濃度分別增加20.2%、8.0%；當(dāng)R=0.50時(shí)，COD質(zhì)量濃度分別增加5.5%、1.7%，氨氮質(zhì)量濃度分別增加18.0%、5.0%；當(dāng)R=0.67時(shí)，COD質(zhì)量濃度分別增加4.9%、1.5%，氨氮質(zhì)量濃度分別增加15.4%、4.2%。根據(jù)C1、C3斷面污染物縱向質(zhì)量濃度變化情況，當(dāng)匯流比分別為0.33、0.50、0.67時(shí)，COD縱向質(zhì)量濃度的梯度變化分別為8.4%、6.6%、5.2%，氨氮縱向質(zhì)量濃度的梯度變化分別為11.8%、9.9%、8.5%。模擬發(fā)現(xiàn)，在交匯后未出現(xiàn)COD明顯的混合均勻位置；當(dāng)R=0.33時(shí)，氨氮在1 820 m處均勻混合；當(dāng)R=0.50時(shí)，氨氮在1 960 m處均勻混合；當(dāng)R=0.67時(shí)，氨氮在2 050 m處均勻混合。可見，匯流比越大，污染物縱向質(zhì)量濃度梯度越小，在交匯區(qū)出現(xiàn)的污染帶越長，出現(xiàn)混合均勻的位置越遠(yuǎn)離交匯口。

圖5 不同匯流比下交匯區(qū)污染物質(zhì)量濃度縱向變化

2.2.2 交匯區(qū)污染物橫向擴(kuò)散特征分析

圖6為不同匯流比下COD、氨氮質(zhì)量濃度在同一橫向斷面上隨不同縱向斷面的變化情況，選取Z1、Z3斷面進(jìn)行橫向擴(kuò)散分析。

由圖6(a)可知，C1至C3斷面污染物質(zhì)量濃度逐漸降低：當(dāng)R=0.33時(shí)，C1至C2斷面、C2至C3斷面COD質(zhì)量濃度分別降低13.6%、5.3%，氨氮質(zhì)量濃度分別降低14.8%、7.8%；當(dāng)R=0.50時(shí)，COD質(zhì)量濃度分別降低14.1%、6.6%，氨氮質(zhì)量濃度分別降低19.3%、10.7%；當(dāng)R=0.67時(shí)，COD質(zhì)量濃度分別降低15.2%、8.0%，氨氮質(zhì)量濃度分別降低21.7%、11.1%。由圖6(b)可知，C1至C3斷面污染物質(zhì)量濃度也逐漸降低：當(dāng)R=0.33時(shí)，C1至C2斷面、C2至C3斷面COD質(zhì)量濃度分別降低3.4%、1.8%，氨氮質(zhì)量濃度分別降低2.3%、1.6%；當(dāng)R=0.50時(shí)，COD質(zhì)量濃度分別降低4.6%、3.2%，氨氮質(zhì)量濃度分別降低3.1%、2.4%；當(dāng)R=0.67時(shí)，COD質(zhì)量濃度分別降低5.7%、4.1%，氨氮質(zhì)量濃度分別降低3.7%、3.1%。根據(jù)Z1、Z3斷面污染物橫向質(zhì)量濃度變化情況，當(dāng)匯流比分別為0.33、0.50、0.67時(shí)，COD橫向質(zhì)量濃度的梯度變化分別為6.0%、7.1%、8.3%，氨氮橫向質(zhì)量濃度的梯度變化分別為6.6%、8.9%、9.9%。模擬發(fā)現(xiàn)，橫向斷面越靠近交匯口，污染物橫向質(zhì)量濃度梯度越大，匯流比的增加會增大污染物的橫向質(zhì)量濃度梯度，尤其對于交匯口附近斷面的影響更為明顯，故而會加劇污染物在交匯口的橫向擴(kuò)散趨勢。

圖6 不同匯流比下交匯區(qū)污染物質(zhì)量濃度橫向變化

3 結(jié)論

較大的匯流比會導(dǎo)致污染帶分布變廣、污染物縱向質(zhì)量濃度梯度變小、縱向擴(kuò)散距離變長，同時(shí)，匯流比的增加會增大污染物橫向質(zhì)量濃度的梯度變化，加劇污染物在交匯口的擴(kuò)散趨勢；MIKE 21 FM模型可較好地反映自然河流交匯區(qū)的特點(diǎn)，可應(yīng)用于自然河流交匯區(qū)污染物擴(kuò)散機(jī)理的研究；其他影響自然河流交匯區(qū)污染物擴(kuò)散規(guī)律的因素如流速、污染物質(zhì)量濃度等有待開展深入研究。