長江下游典型河段入江污染因子擴散模擬研究

摘要：基于長江下游一沿江企業(yè)擴建優(yōu)化項目，采用CJK3D-WEM建立了水動力-水質(zhì)模型，模擬和評估不同水文條件下（枯水期和豐水期）主要污染因子（COD、氨氮、總磷）的擴散情況。結果表明：企業(yè)擴排后，污染因子擴散影響范圍較小，主要局限于排口上下游近岸水域；枯水期的最大擴散距離為2 300 m，豐水期為350 m，均符合所在水功能區(qū)的水質(zhì)管理要求；污染因子的影響距離與濃度增量之間存在較好的函數(shù)關系；豐水期稀釋作用明顯，污染物濃度降低所需距離較短，而枯水期則相反。為兼顧環(huán)境保護和經(jīng)濟發(fā)展，建議優(yōu)化水功能區(qū)的排污總量控制，實時監(jiān)測水質(zhì)和排污情況，實施動態(tài)監(jiān)管和調(diào)配。

關 鍵 詞：污染物因子；擴散模擬；水動力-水質(zhì)模型；CJK3D-WEM；數(shù)學模型；長江

中圖法分類號：X52 文獻標志碼：ADOI：10.16232/j.cnki.1001-4179.2024.S2.004

0 引言

為保證工業(yè)發(fā)展的先進性，增設新產(chǎn)業(yè)線完成內(nèi)部改革是較多企業(yè)優(yōu)化生產(chǎn)結構的重要手段。然而增設新產(chǎn)業(yè)線也增加了一定污染物的排放，從而給周邊水域水環(huán)境治理帶來一定壓力^［1］。因此論證產(chǎn)業(yè)結構優(yōu)化升級方案的可行性，需要模擬污染物排放對周邊水域造成的影響^［2-4］，以確保方案實施后，周邊水域的水質(zhì)等級仍滿足要求。污染物因子主要包括COD、氨氮與總磷。COD含量過高將導致水體富營養(yǎng)化^［5］、水生生物死亡、水質(zhì)惡化^［6］以及生態(tài)系統(tǒng)失衡等一系列惡劣影響。氨氮會影響魚類的攝食與成長，含量高時會導致魚類死亡，人類長期飲用氨氮超標水體易引發(fā)癌癥^［7］。總磷對魚類、底棲動物群落健康具有顯著影響，是水體富營養(yǎng)化的重要誘因^［8］。因此，針對污染物因子的擴散模擬研究對于周邊水域的水環(huán)境保護至關重要。

數(shù)學模型在理解和預測污染物行為方面發(fā)揮著重要作用，有助于環(huán)境管理和決策。這些模型從簡單的解析表達式到復雜的數(shù)值模擬，涵蓋了針對特定污染物和環(huán)境條件的各種工況。對流-擴散模型描述了污染物通過對流（由于平均流動的傳輸）和擴散（由于隨機運動的擴散）的運輸^［9］，適用于水污染場景。模型中的污染物濃度由偏微分方程控制，考慮了對流速度、擴散系數(shù)和反應速率。計算流體力學模型（CFD模型）使用數(shù)值方法求解描述流體運動的Navier-Stokes方程，并結合污染物傳輸?shù)膶α?擴散方程，模擬復雜幾何形狀和多變環(huán)境條件下的污染物擴散^［10］。CFD模型提供了詳細的污染物濃度時空分布，適用于城市空氣質(zhì)量研究和工業(yè)排放場景。拉格朗日模型追蹤單個污染物粒子或包裹在環(huán)境中的運動，適用于研究湍流中污染物的擴散，其中粒子的軌跡可能非常不規(guī)則，拉格朗日模型可以結合隨機過程以考慮粒子運動中的隨機波動^［11］。針對水污染，綜合水文模型結合了地表水和地下水系統(tǒng)，以模擬流域和含水層中的污染物運輸^［12］。該類模型通常需要求解對流-擴散-反應方程，并結合入滲、徑流和地下水流等水文過程，它們對于評估污染物對水質(zhì)和生態(tài)系統(tǒng)健康的影響至關重要。綜上所述，污染物擴散的數(shù)學模型為預測環(huán)境中污染物的傳播提供了寶貴工具，模型的選擇取決于特定污染物、環(huán)境介質(zhì)及其時空尺度。目前，針對COD、氨氮與總磷的擴散模擬，使用較多的為CFD模型，通過水動力-水質(zhì)數(shù)學模型進行污染物擴散模擬，該方法在萊州灣^［13］、渤海灣^［14］以及欽州灣^［15］等區(qū)域的污染物擴散研究中有較為成熟的應用。

本文選取長江下游一沿江企業(yè)擴建優(yōu)化項目開展入江排污影響研究，采用CJK3D-WEM模型模擬計算COD、氨氮和總磷等污染物因子的擴散。對各污染因子增量敏感度進行統(tǒng)計分析，研究結果可為后續(xù)水污染治理與水資源優(yōu)化配置提供參考。

1 研究范圍及基本概況

該沿江企業(yè)入江排口位于長江下游江西省九江市湖口縣境內(nèi)，排口段涉及長江下游九江河段中下段（張家洲河段、上下三號河段）。張家洲汊道自鎖江樓至八里江口，為彎曲分汊型河道，鄱陽湖在張家洲右汊末端匯入長江，右汊為主汊分流比約為59.9%。上下三號河段自八里江口至小孤山，為微彎多分汊型河道，自上而下分布有新洲、上三號洲、下三號洲，上三號洲已并入北岸。

位于排口上游的九江水文站為距其最近的長江干流水文站點，九江站地處長江中游干流江西省九江市，其下游32.1 km有鄱陽湖入?yún)R長江，對該站水位、流量有明顯頂托影響。據(jù)統(tǒng)計，1957～2023年，九江站歷年最高水位23.03 m（凍結基面，1998年8月2日），最低水位6.83 m（凍結基面，1963年2月12日）；三峽水庫蓄水前（1957～2002年），多年平均水位13.60 m；三峽水庫蓄水后（2003～2023年），多年平均水位12.94 m。1988～2023年，九江站歷年最大流量75 000 m³/s（1996年7月23日），最小流量5 850m³/s（1999年3月27日）；三峽水庫蓄水前（1988～2002年），多年平均流量23 800 m³/s；三峽水庫蓄水后（2003～2023年），多年平均流量22 300 m³/s。九江站全年徑流量主要集中在汛期，三峽水庫蓄水前和蓄水后5～10月的徑流量約占全年徑流量的71.3%和67.9%，全年7～9月為高水位期，1～3月為低水期。

該排口位于長江干流右岸長江湖口、彭澤保留區(qū)內(nèi)，水質(zhì)現(xiàn)狀與管理目標均是Ⅲ類。排口下游是長江彭澤飲用水源區(qū)，其中彭澤二水廠保護范圍內(nèi)的水質(zhì)現(xiàn)狀與管理目標均是Ⅱ類。本次模型研究范圍為長江湖口下游至長江彭澤飲用水源區(qū)終點處，研究區(qū)域與監(jiān)測站位分布見圖1。據(jù)長江湖口、彭澤保留區(qū)人工水質(zhì)監(jiān)測斷面2021年1月至2023年4月水質(zhì)監(jiān)測數(shù)據(jù)，該區(qū)域28個測次中，Ⅱ類水水質(zhì)27次，占總監(jiān)測頻次的96.4％，Ⅲ類水水質(zhì)1次，占總監(jiān)測頻次的3.6％，全部達Ⅲ類水評價標準；據(jù)集中式生活飲用水源地2021年1月至2023年4月水質(zhì)監(jiān)測數(shù)據(jù)，長江彭澤飲用水源區(qū)10個測次，全部達Ⅱ類水評價標準。監(jiān)測結果表明，研究區(qū)域水質(zhì)目前全部達標。然而區(qū)域內(nèi)企業(yè)產(chǎn)業(yè)結構亟需優(yōu)化升級，需增設生產(chǎn)線，為了促進區(qū)域水資源可持續(xù)利用和經(jīng)濟社會高質(zhì)量發(fā)展，控制擴排項目對水環(huán)境造成的不利影響，需在現(xiàn)狀基礎上對新增污染物進行擴散模擬，以保障受納水體所在水功能區(qū)的生活、生產(chǎn)和生態(tài)用水安全。

2 水域污染物排放數(shù)學模型

2.1 二維污染物排放模型

基于二維有限體積法的水環(huán)境數(shù)學模型CJK3D-WEM，建立研究區(qū)域二維COD、氨氮及總磷水質(zhì)模型，水流、水質(zhì)運動方程為

式中：z為水位；t為時間；H為總水深；u，v為流速矢量沿x，y方向的速度分量；f為科氏系數(shù)；g為重力加速度；N_x，N_y為x，y向水流紊動黏性系數(shù)；C為污染物濃度；A_x，A_y為x，y向水質(zhì)擴散系數(shù)；S為水質(zhì)源匯項；k為水質(zhì)衰減系數(shù)。

水動力模型主要為水質(zhì)模型提供水動力場，驅動COD、氨氮及總磷在時間和空間上產(chǎn)生變化，采用三角形網(wǎng)格對計算區(qū)域進行離散，將單一的網(wǎng)格單元作為控制單元，采用有限體積法對式（1）進行離散求解^［16］。

數(shù)學模型模擬范圍上溯八里江水位站西側約3 km，下至彭澤水位站下游約2 km，范圍總長約33 km，見圖1。采用三角形網(wǎng)格作為計算單元，模型網(wǎng)格總數(shù)為155 741個，工程附近水域網(wǎng)格加密，整個模型最大網(wǎng)格邊長54 m，最小網(wǎng)格邊長約8 m。模型水深采用大范圍實測地形，工程附近水域采用2023年最新實測地形，高程采用1985國家高程基準。開邊界給定潮位、溫度、鹽度、COD、氨氮及總磷污染物質(zhì)量濃度值，上游邊界采用流量控制，下游邊界采用水位控制。其中，計算時間步長取0.3 s，糙率采用附加糙率方式進行處理，基本糙率為0.021，附加糙率取值0.015，水流紊動黏性系數(shù)取0.1HU_*（H為水深，U_*為摩阻流速），動邊界水深取0.01 m。水質(zhì)計算的關鍵參數(shù)為不同污染物的降解系數(shù)選取，已有學者對中國21條河流的資料進行回歸分析后，得到有機污染物自然降解速率K的計算公式為

K_COD=0.65×Q^-0.15

K_氨氮=1.8×Q^-0.49

K_總磷=0.5586×Q^-0.15（2）

式中：K為河流中污染物降解系數(shù)，1/d；Q為河流流量，m³/s。

經(jīng)計算枯水期至豐水期，長江段K_COD=0.138～0.167，K_氨氮=0.011～0.021，K_總磷=0.119～0.143。水質(zhì)擴散系數(shù)取0.1HU_*。水動力模型驗證資料選取2023年最新水文測驗資料，水質(zhì)模型驗證資料選取2022年第四季度水質(zhì)監(jiān)測資料。

2.2 模型模擬驗證

2.2.1 水動力驗證

采用2023年水文測驗資料作為數(shù)學模型驗證資料，水動力驗證點位置如圖1所示。本次水文測驗期間上游來水流量約24 000 m³/s。表1為水位和流量計算值與實測值對比，圖2為斷面流速對比。由圖2和表1可知：數(shù)學模型水位計算值與實測值差值基本在0.04 m以內(nèi)，斷面流速分布趨勢數(shù)模計算值與實測值基本一致，河道主汊流量偏差在3%以內(nèi)，滿足規(guī)范要求。本文建立的水動力模型可以較好地反映研究水域的水動力情況，可為水質(zhì)模型的建立提供基礎。

2.2.2 水質(zhì)驗證

（1）排污口輸入。根據(jù)現(xiàn)狀排污口的排污情況，將現(xiàn)有排污口的排污量作為輸入條件，排污口位置見圖1，具體排污量如表2所列。排水污染物濃度按照一級A標準，COD取值50 mg/L，氨氮取值8 mg/L，總磷取值0.5 mg/L。

（2）上游來水條件。根據(jù)前期監(jiān)測結果輸入模型上邊界流域來水及水質(zhì)情況。由于枯水期多是水污染較為嚴重的時期，因此選用2022年第四季度的監(jiān)測結果。監(jiān)測結果顯示，參照現(xiàn)狀排污口上游500 m處的水質(zhì)監(jiān)測情況，上游來水水質(zhì)COD濃度為15.5 mg/L，氨氮濃度為0.1 1 mg/L，總磷濃度為0.10 mg/L。上游流量在9 000～10 000 m³/s之間。

（3）水質(zhì)驗證。在現(xiàn)狀排口上游500 m（1號）、下游500 m（2號）、下游1 500 m（3號）、下游3 000 m（4號）進行了水質(zhì)監(jiān)測，驗證點位置見圖1。污染物濃度數(shù)學模型計算值與實測值對比見圖3，可見除個別點偏差較大外，數(shù)學模型的計算值與實測值偏差較小，說明數(shù)學模型計算參數(shù)取值合理。

3 結果與討論

3.1 計算條件與工況

3.1.1 徑流條件

根據(jù)HJ 2.3—2018《環(huán)境影響評價技術導則地表水環(huán)境》，邊界輸入條件“河流不利枯水條件宜采用90%保證率最枯月流量或近10 a最枯月平均流量”。采用三峽水庫蓄水后九江水文站2002～2023年最枯月均流量進行P-Ⅲ曲線適線，90％保證率最枯月均流量為8 630 m³/s，根據(jù)2023年水位流量關系對應的水位z=6.05 m，九江站2023年水位流量關系擬合見圖4。據(jù)統(tǒng)計，九江站近10 a最枯月平均流量為8 890 m³/s，根據(jù)2023年水位流量關系對應的水位z=6.16 m。本次數(shù)模計算采用不利情況進行預測：90%保證率最枯月流量8 630 m³/s，對應的水位z=6.05 m。同時，根據(jù)河道的評價等級，按照HJ 2.3—2018要求，還應計算豐水期污染物擴散范圍，豐水期流量取2002～2023年九江站洪季5～10月平均流量，約為30 700 m³/s，對應水位約為13.70 m。

3.1.2 水質(zhì)因子與背景濃度

預測內(nèi)容為正常排放工況下的COD、NH₃-N、TP濃度。根據(jù)排放標準限值，污染物排放流量為40 000 m³/d（0.463 m³/s），COD、NH₃-N、TP排放源強依次為50，8，0.5 mg/L。污染物濃度計算至輸移穩(wěn)定。

將數(shù)學模型驗證的上游來流水質(zhì)條件計算結果作為水質(zhì)的背景濃度。排口附近水域COD濃度平均值約為14.5 mg/L，氨氮約為0.11 mg/L，總磷約為0.09 mg/L。

3.1.3 計算工況

根據(jù)徑流條件與項目廢水排水源強，確定本次污染物數(shù)值模擬的計算工況包括：①枯水期（8 630 m³/s）；②豐水期（30 700 m³/s）。

3.2 排污擴散結果

在數(shù)學模型中設置多個采樣點進行結果統(tǒng)計，采樣點分布如圖5所示。其中SY1、SY2距離排口250，500 m，SY2～SY10之間各采樣點距離500 m。

3.2.1 枯水期

污染物形態(tài)呈現(xiàn)出長條形，基本沿長江南岸水域分布。COD因子由于排放濃度高，因此污染物增量范圍最大，TP范圍最小。據(jù)統(tǒng)計，COD增量超過0.1 mg/L的水域長度約為446 m，寬度約為100 m，COD增量超過0.03 mg/L的水域長度約為2 300 m；NH₃-N增量超過0.0 1mg/L的水域長度約為783 m，寬度約為110 m，NH₃-N增量超過0.005 mg/L的水域長度約為2 250 m；TP增量超過0.001 mg/L的水域長度約為475 m，寬度約為80 m。根據(jù)數(shù)學模型計算結果，各水質(zhì)因子增量影響范圍分布在現(xiàn)有排口沿岸方向最大上溯距離為43 m，沿岸方向最大下延距離為2 300 m范圍內(nèi)（圖6～8）。

表3為正常排放條件下的敏感目標濃度增量值，可以看出：各水質(zhì)因子基本聚集在排口下游附近，彭澤第二水廠取水口和水功能區(qū)起點COD、NH₃-N、TP增量較小，未改變水質(zhì)分類值；長江省控自動水質(zhì)監(jiān)測斷面、長江省控人工水質(zhì)監(jiān)測斷面處COD、NH₃-N、TP增量較小，未改變水質(zhì)分類值。

3.2.2 豐水期

圖9～11為豐水期正常排放條件下的污染物分布及增量包絡線圖。由于豐水期水位高、流量大，因此污染因子稀釋快，故豐水期各水質(zhì)因子影響范圍遠小于枯水期。COD因子由于排放濃度高，因此污染物增量范圍最大，TP范圍最小，據(jù)統(tǒng)計，COD增量超過0.1 mg/L的水域長度約為100 m，寬度約為32 m，COD增量超過0.03 mg/L的水域長度約為650 m，寬度約為53 m；NH₃-N增量超過0.01 mg/L的水域長度約為201 m，寬度約為57 m，NH₃-N增量超過0.005 mg/L的水域長度約為627 m，寬度約為84 m；TP增量超過0.001 mg/L的水域長度約為91 m，寬度約為29 m。

表4為正常排放條件下的敏感目標濃度增量值，可以看出：各水質(zhì)因子基本聚集在排口下游附近，彭澤第二水廠取水口和水功能區(qū)起點COD、NH₃-N、TP基本無增量；長江省控自動水質(zhì)監(jiān)測斷面、長江省控人工水質(zhì)監(jiān)測斷面處COD、NH₃-N、TP亦基本無增量。

3.3 各污染因子增量敏感度分析

根據(jù)上述擴散結果，統(tǒng)計枯水期各污染因子增量包絡線影響距離見表5，豐水期見表6。

根據(jù)2.1節(jié)降解系數(shù)與流量之間關系，推測污染因子增量包絡線影響距離與增量為allometricl函數(shù)關系，即有：

y=A·x^B（3）

式中：A，B為待定系數(shù);y為各污染因子增量包絡線影響距離;x為各污染因子濃度增量。不同污染因子包絡線影響距離與濃度增量關系曲線如圖12所示，豐水期總磷由于數(shù)據(jù)量少沒有進行擬合。可以看出，該模型能較好地模擬包絡線影響距離與濃度增量之間的關系，待定系數(shù)A與工況、污染因子有關，豐水期小于枯水期，CODgt;氨氮gt;總磷，在0.22～21.29之間，COD枯水期最大，為21.29；B與工況有關，豐水期小于枯水期，在-1.65～-1.11之間。可見豐水期對污染物的稀釋作用較大，削減濃度所需要的距離更短，而枯水期污染物擴散范圍較大，相同的濃度削減需要更多的距離保證，因此在枯水期應重視海域的水環(huán)境保護，采取一定的措施削減污染物的擴散。

4 結論

（1）企業(yè)擴排后，COD、氨氮、總磷污排水影響范圍較小，僅限于排口上、下游兩側近岸水域，呈扁長狀貼岸分布形態(tài)，COD、氨氮、總磷增量影響范圍分布在現(xiàn)有排口沿岸方向枯水期最大擴散距離為2 300 m，豐水期為350 m，疊加增量后不影響所在水功能區(qū)水質(zhì)現(xiàn)狀和管理目標要求。

（2）污染因子增量包絡線影響距離與污染因子濃度增量呈現(xiàn)較好的allometricl函數(shù)關系，豐水期對污染物的稀釋作用較大，削減濃度所需要的距離更短，而枯水期污染物擴散范圍較大，相同的濃度削減需要更多的距離保證。

（3）為了兼顧資源環(huán)境保護與經(jīng)濟社會高質(zhì)量發(fā)展，應進一步梳理和優(yōu)化水功能區(qū)限制排污總量，并實時掌握該段排污和水功能區(qū)水質(zhì)情況，實行統(tǒng)一動態(tài)監(jiān)管與調(diào)配，保障水源安全與企業(yè)的合理排污需求。

參考文獻：

［1］ 畢雪，潘婷婷.蓄滯洪區(qū)安全建設工程水環(huán)境影響分析［J/OL］.人民長江，1-7［2024-10-08］.http：∥kns.cnki.net/kcms/detail/42.1202.TV.20240726.1418.008.html.

［2］ 肖新宗，郭芳，賈慶林，等.南水北調(diào)中線總干渠污染物降解及自凈能力研究［J］.人民長江，2024，55（5）：49-56.

［3］ 張云，蔡彬彬.污水處理廠達標尾水導流排江可行性研究：以南通市益民污水處理廠為例［J］.人民長江，2020，51（11）：53-57.

［4］ 韓非非，崔冬.長江口污水超標排放對水質(zhì)影響的數(shù)值模擬研究［J］.人民長江，2018，49（14）：17-23.

［5］ 尹華，昌鏡偉，章光新，等.新立城水庫水體富營養(yǎng)化成因及治理對策［J］.東北師大學報（自然科學版），2010，42（1）：152-156.

［6］ JUNG KW，HWANG MJ，PARK DS，et al.Combining fluidized metal-impregnated granular activated carbon in three-dimensional electrocoagulation system：feasibility and optimization test of color and COD removal from real cotton textile wastewater［J］.Separation and Purification Technology，2015，146：154-167.

［7］ 袁行知，許雪峰，俞亮亮，等.基于水動力水質(zhì)模型的平原河網(wǎng)排污模擬分析［J］.中國農(nóng)村水利水電，2022（12）：39-46.

［8］ 婁保鋒，余明星，黃波，等.1998～2019年三峽水庫總磷濃度與形態(tài)變化及通量匯算［J］.北京大學學報（自然科學版），2024，60（3）：548-562.

［9］ CUNHA C，CARRER J，OLIVEIRA M，et al.A study concerning the solution of advection-diffusion problems by the Boundary Element Method［J］.Engineering Analysis with Boundary Elements，2016，65：79-94.

［10］LAUNDER BE，SPALDING DB.The numerical computation of turbulent flows［M］.Amsterdam：Elsevier，1983.

［11］SAWFORD B.Turbulent relative dispersion［J］.Annual review of fluid mechanics，2001，33（1）：289-317.

［12］BEAR J，CHENG AH D.Modeling groundwater flow and contaminant transport［M］.Berlin：Springer，2010.

［13］李婉華，宋彥，冷星，等.小清河入海污染物擴散特征及其對近岸工程建設的響應［J］.海洋地質(zhì)前沿，2024，40（4）：47-59.

［14］劉浩，戴明新，彭士濤，等.渤海灣主要污染物環(huán)境容量的估算［J］.海洋通報，2011，30（4）：451-455.

［15］郭森，韓保新，楊靜，等.納污海域水環(huán)境容量計算與總量分配方法研究：以欽州灣為例［J］.環(huán)境科學與技術，2006（增1）：19-22.

［16］路川藤，羅小峰.基于非結構網(wǎng)格的高分辨率隱式算法研究及應用［J］.海洋通報，2015，34（1）：59-64.

（編輯：劉媛）