黃河干支流交匯河段污染物傳輸模擬研究

余真真，朱 海，毛頌平，閆 莉，張 雷，趙麗萍

（1.黃河水資源保護科學研究院，河南 鄭州 450004； 2.河海大學 水利水電學院，江蘇 南京 210098；3.黃河水利科學研究院，河南 鄭州 450003）

黃河是我國西北、華北地區重要的供水水源，黃河流域是我國糧食主產區、重要的能源和工業基地，在國家經濟社會發展和生態安全方面具有十分重要的地位［1－2］。 黃河流域水資源短缺，水環境承載能力低，以支流水環境污染為特征的水危機未得到根本扭轉［3－6］，水資源超載與耗水型產業布局，河道季節性斷流與廢水處理投入不足等問題交織。 目前黃河流域水環境質量穩中趨好，重污染河段集中在中游支流城鎮的下游，枯水期污水處理廠廢污水排入，難以實現水功能區水質目標。

黃河干流污染物包括入河排污口排入的污染物，同時包括大量通過支流排入干流的污染物。 支流水質差且含沙量高，輸污比重遠大于干流，支流水質污染和來水量銳減是流域水污染形勢嚴峻的主要原因［7］。據統計，入黃支流每年向黃河輸送化學需氧量（COD）約占黃河化學需氧量的70.4%，氨氮約占67.7%。 污染物通過支流入黃口匯入干流，導致干流水體存在污染風險，極大地影響著沿黃地區供水安全。 同時，干支流交匯河段水體摻混劇烈、水流條件特殊，流速分布不均，特殊的水力特性會對泥沙、污染物產生較強的滯留作用，污染物輸運過程極其復雜［8－9］。 倪晉仁等［10］針對交匯河段水力計算進行了探討；茅澤育等［11－12］實施了明渠交匯口三維水力特性試驗研究和污染物輸移數值計算；馮鏡潔等［13］開展了河流交匯分離區特性研究；魏文禮等［14］開展了不同交匯角度明渠交匯口三維水力特性研究。 然而，上述研究大都采用概化的順直河道進行模擬，未結合河流交匯口實際形態進行模擬。本文以黃河小北干流與汾河入黃河交匯河段為研究區域，采用黃河干支流實測地形，建立實尺度大渦模擬模型，對交匯河段污染物沿程傳輸過程進行研究，捕捉水體紊動誘發污染的三維混合與擴散過程，探索黃河水污染特性與規律。

1 黃河干支流水質狀況

（1）流域整體水質演變。 從20 世紀80 年代開始黃河流域具有較完整水質評價資料，根據各時期水質評價結果可知，全流域水質狀況經歷了良好、污染嚴重、污染得到初步遏制的歷程［15］，見圖1。 2000 年前后是全流域水質最差的時段，依據2003 年水質評價結果，流域年均劣Ⅴ類水河長占32.0%，而Ⅰ～Ⅲ類水河長僅占33.5%；甘陜、寧蒙、晉陜豫等省界水體污染相當嚴重，劣Ⅴ類水占37.9%；不符合飲用要求的重點飲用水源地達90.0%。 近年來，流域水質惡化趨勢基本得到遏制，全流域總體水質狀況逐步好轉。 與2003 年相比，2018 年全流域Ⅰ～Ⅲ類水河長比例提高到了73.5%，劣Ⅴ類水比例下降到了18.9%。 黃河流域水功能區水質達標率由2005 年的37.5%上升到了2018年的63.3%。

（2）干支流水質變化。 統計干流水質表明，2000年干流Ⅰ～Ⅲ類水河長占54.7%，Ⅳ～Ⅴ類水河長占40.5%，劣Ⅴ類水河長占4.8%。 2018 年干流Ⅰ～Ⅲ類水河長占比升至97.8%，Ⅳ～Ⅴ類水河長占比下降至2.2%，見圖2。 自2009 年起黃河干流連續多年未出現劣Ⅴ類水。 黃河干流重要河段水體主要污染物質量濃度大幅下降，上游蘭州、石嘴山、畫匠營，中游潼關、三門峽，下游花園口等重要斷面水體主要污染物質量濃度逐步下降。

2005 年，黃河支流Ⅰ～Ⅲ類水河長占比37.9%，劣Ⅴ類水河長占比高達41.7%。 其后，總體水質逐漸改善，至2018 年，支流Ⅰ～Ⅲ類水河長占比已上升至66.4%，劣Ⅴ類水河長占比下降為16.1%，見圖3。2000 年以來，以前污染嚴重的汾河西師、渭河吊橋、洛河七里鋪等支流控制斷面主要污染物質量濃度均大幅下降。 其中汾河西師斷面COD 質量濃度由2000 年的569.60 mg／L 下 降 至2018 年42.17 mg／L，下 降 了92.60%，氨氮質量濃度由2000 年16.32 mg／L 下降至2018 年0.90 mg／L，下降了94.49%。

圖3 2000—2018 年黃河支流水質變化

2 研究區域水環境特性

黃河中游禹門口至潼關河段是托克托至潼關北干流的下段部分，稱小北干流，全長132.5 km，為晉陜兩省的界河。 黃河出禹門口（也稱“龍門”）后，驟然放寬，河床由100 m 的峽谷展寬至4 km 以上，兩岸分布有大量灘地，經小北干流后在潼關河寬收縮為850 m，折向東流。 禹潼河段屬淤積性游蕩型河道，洪水具有峰高量大、含沙量高的特點。 泥沙大量淤積，河道寬淺，水流散亂，主流游蕩不定。 該河段是晉陜兩省交界的水質敏感區，處在承上啟下控制黃河水質的關鍵部位，龍門下游約55 km 處，流經山西省的汾河在黃河左岸匯入黃河干流。 汾河與黃河交匯處河道形態見圖4，其中x、y分別為研究區域相對坐標，x為河長、y為河寬。

圖4 汾河與黃河交匯處河道形態

汾河是黃河第二大支流，發源于寧武縣東寨鎮管涔山，在萬榮縣廟前村附近匯入黃河。 干流全長694 km，流域面積39 471 km2，多年平均徑流量18.47 億m3，實測多年平均輸沙量0.22 億t。 汾河屬于資源性缺水地區，隨著經濟社會的發展，水資源供需矛盾突出，存在地下水超采、河道斷流、水污染嚴重等問題，汾河入黃口西師斷面水質不能穩定達標。

圖5 表明，黃河干流龍門斷面2014—2018 年連續5 a 的COD、氨氮質量濃度滿足Ⅲ類水質目標要求，多年平均COD、氨氮質量濃度分別為12.95、0.25 mg／L；汾河西師斷面2014—2018 年的COD、氨氮質量濃度均未能達到Ⅳ類水質目標要求，多年平均COD、氨氮質量濃度分別為42.06、6.20 mg／L。 汾河入黃斷面COD、氨氮多年平均質量濃度約為黃河干流龍門斷面的3.25倍和24.8 倍。

圖5 黃河干流龍門、支流汾河西師斷面污染物質量濃度箱式圖

3 數學模型構建

（1）控制方程。 自然界中河流水流形態以紊流形式存在，紊流是十分復雜的多尺度不規則流動，該性質導致湍流脈動復雜的非線性輸運特性。 目前在洪水預報和工程設計中常采用雷諾平均數值模擬技術模擬水流運動，但其只能提供湍流的平均信息，計算精度較低［16］。 紊流除存在許多隨機性很強的小尺度渦運動外，還存在一些大尺度擬序結構，大渦對于質量、動量與能量的傳輸起主要作用，對流動的初始條件和邊界形狀有強烈依賴性［17］。 本文采用大渦數值模擬方法，對大尺度湍流直接進行數值模擬，對小尺度湍流脈動建立模型模擬較高雷諾數和較復雜的湍流運動，得到污染物在空間上的不對稱運動特性，獲得更多物質傳遞信息，實現湍流及其標量輸運過程的高精度模擬［18－19］。

不可壓常黏性系數的紊流運動大渦模擬控制方程為濾波后的Navier?Stokes 方程：

式中：x為空間坐標；u為流速；＝／2，為濾波后的速度應變率張量；ν為黏性系數；ρ為流體密度；t為時間；p為壓力；τij為亞格子應力張量；τkk為亞格子正應力張量之和；δij為克羅內克符號；νSGS為亞格子渦黏系數。

研究河段為黃河干流禹門口至潼關河段，屬于寬淺河道，河流水平尺度遠大于垂直尺度，將水平大渦模擬模型作為首選紊流模型。 對水平方向的亞格子渦黏系數和渦擴散系數建立模型：

式中：ks為尼古拉茲粗糙長度；γ為斯蒂芬－玻爾茲曼常數；σT為普朗特－施密特數；S?為應變率的平方和；U為水深平均水平速度矢量；C為謝才系數；H為水深；DSGS為渦擴散系數。

考慮污染物濃度的變化過程，其控制方程為濾波后的對流擴散方程：

（2）網格剖分。 為了精細模擬交匯河口污染物的輸運規律，在計算區域內，沿河流縱向和橫向分別劃分308、36 個計算網格，以保證網格盡可能平順過渡，沿水深方向劃分10 層網格，網格共計11.09 萬個。 干支流交匯處進行網格加密，在汾河入黃河口外延31×16個網格與黃河干流計算網格搭接。

（3）數值離散。 采用循環隱式進程的有限差分ADI 算法求解控制方程，每一時間步交替地沿x、y方向掃描，具體做法：從l步到l＋1／2 步，將x方向的動量方程和連續方程只沿x方向形成關于水位ζ和x方向流速u的隱式方程組，形成三對角矩陣，可用追趕法求解每一行（l＋1／2）Δt時刻的水位ζ和流速u，將求得的水位和流速代入y方向的動量方程顯式解出y方向流速v；從l＋1／2 步到l＋1 步，將y方向的動量方程和連續方程只沿y方向形成關于水位ζ和流速v的隱式方程組，同理可用追趕法求解每一列（l＋1）Δt時刻的水位ζ和流速v，顯式求解x方向的動量方程。

（4）邊界條件。 床面采用無滑移邊界條件：計算區域較大時，可忽略邊壁對水流的阻力，側向采用無滑移固壁邊界，其法向動量通量和標量通量均為0；水動力場的入口邊界給定流量值、濃度場入口邊界給定污染物濃度值，出流水動力邊界設置為水位，污染物濃度設置為自由出流，各標量梯度為0。

4 計算結果分析

4.1 模擬工況設置

結合2017 年實測地形數據，根據當年水動力條件，以COD、氨氮為代表水質因子進行分析。 水動力上邊界、下邊界分別為流量、水位過程，汾河口作為旁側入流邊界設定；開邊界同時給定COD、氨氮實測濃度邊界。 模擬工況邊界條件見表1。

表1 模擬工況邊界條件

4.2 水動力及污染物濃度平面變化

實尺度大渦模擬模型計算表明，汾河河口與黃河干流交匯區域，交匯口上游干流河槽過水斷面突然擴大、交匯后又突然縮窄，該復雜的交匯口平面幾何形態導致特殊的水動力特性，對污染物輸移產生了一定影響。 汾河與黃河干流流量比約為1 ∶15，經統計模擬范圍平均流速（矢量和）約為0.39 m／s，河道主泓流速較大，約為0.60 m／s，邊灘流速較小。 交匯口橫斷面水流流速不均勻系數為0.973，水流紊動能最大值為0.036 m2／s2。 干流深泓平面形態呈U 形，交匯口位于彎道的凹岸，匯合口下游左岸形成了長約1.75 km、寬約1 km 的低流速區，導致污染物滯留。 計算區域渦量等值線見圖6。 交匯口處高渦量區集中在黃河干流左岸，經過交匯口后，逐漸向右岸轉移。 渦量值為0～0.031／s，最大渦量值出現在x＝510 km 附近，主要是受地形影響，交匯口下游11 km 處河道突然縮窄，水流流速增大，加快了污染物的擴散。

圖6 渦量等值線（單位：s－1）

汾河與黃河交匯處COD、氨氮年均質量濃度平面輸移云圖見圖7。 由圖7 可知，支流汾河高濃度污染水體進入黃河后，黃河干流污染物質量濃度顯著增大，COD、氨氮質量濃度最大值分別為53.5、5.8 mg／L。 匯流后受地形影響，左岸形成污染物滯留區，COD、氨氮質量濃度分別為47.02、5.02 mg／L，隨后污染物呈帶狀逐漸偏向右岸，跨越省界到達陜西省境內，距離汾河入黃河口下游約10 km 處又偏向左岸隨水流向前輸移，期間干支流水體并未完全混合。

圖7 COD 和氨氮年均質量濃度輸移云圖

4.3 交匯區污染物三維混合擴散特性

（1）污染物三維混合特征。 經模擬計算得到汾河入黃河交匯口下游河段污染物三維空間分布。 研究河段屬于寬淺河道，污染物質量濃度分布不存在明顯的垂向梯度，但是具有較強的橫向及縱向不均勻性；河口斷面（x＝499.8 km）靠近交匯口一側（河道左岸）污染物質量濃度較高，且污染物摻混主要發生在左岸；隨著汾河高濃度污染水體與黃河干流水體互相摻混，下游x＝501 km 處高質量濃度污染物逐漸向右輸移；x＝502 km 處污染物擴散范圍更廣，已經達到斷面中心，分布更加均勻。

采用不均勻性指數Dev（x）定量分析交匯河口污染物的混合特征，以評估污染物混合程度。 不均勻性指數越小，混合越均勻。

式中：Cp為平均預測濃度；Cs（x，y）為坐標（x，y）處的濃度。

汾河河口及黃河下游沿主泓線斷面COD 和氨氮濃度分布不均勻性指數見圖8，可以看出，不均勻性指數沿程呈減小－平穩－減小趨勢。 COD 濃度Dev（x）最大值為0.44，氨氮濃度Dev（x）最大值為1.61，COD 濃度空間分布比氨氮更均勻，兩者均在x＝510 km 附近（河道斷面縮窄段）充分混合。

圖8 COD 和氨氮濃度不均勻性指數沿程變化情況

（2）污染物縱向擴散。 為了定量研究典型污染物向周邊水體的擴散規律，沿黃河干流泓道走向提取污染物沿程質量濃度，見圖9。 圖中X1 為汾河入黃交匯口位置，X2 為沿主泓線各污染物質量濃度達到平均預測質量濃度時的位置，表示該處污染物已充分混合。計算結果表明，交匯口上游黃河干流主泓沿線COD、氨氮質量濃度基本為恒定值，且與入口邊界質量濃度基本一致；支流高濃度污染水體的流入使得交匯口下游污染物質量濃度顯著抬升，在x＝505 km 處，COD、氨氮質量濃度均達到最大值；在干支流混合作用下，下游污染物質量濃度逐漸趨于恒定，各污染物質量濃度達到充分混合狀態時的位置距離交匯口11 km。

圖9 汾河入黃河交匯口下游COD、氨氮質量濃度縱向分布（沿干流主泓線）

（3）污染物橫向擴散。 在汾河入黃匯合口下游布置4 個監測斷面觀察污染物橫向擴散特性監測斷面見圖7，圖中曲線為1／2 質量濃度等值線范圍，即黃河干流與汾河污染物質量濃度之和的一半。 COD 及氨氮1／2 質量濃度等值線對應質量濃度分別為31.960、2.891 mg／L。 在匯合口附近，汾河高質量濃度COD、氨氮水體匯入，斷面左岸水體COD 和氨氮質量濃度高，右岸水體質量濃度基本與黃河干流上游水體污染物質量濃度一致，此時混合主要發生在左岸（斷面C1、C2）；隨著污染物的混合與輸移，左岸水體COD、氨氮質量濃度漸漸下降，斷面質量濃度漸漸趨于平穩（斷面C3）；在通過下游突然縮窄斷面（x＝510 km處）后，COD、氨氮在橫向擴散上達到充分摻混（斷面C4）。 經測算，汾河河口附近COD 和氨氮1／2 質量濃度等值線包絡范圍基本一致，沿流向長約3 200 m、橫向寬約750 m，COD 和氨氮具有相似的混合規律。

COD、氨氮的橫向斷面分布見圖10 至圖12，總體上污染物到達C4 斷面后均已發生充分摻混。 對于靠近河口的C1、C2 斷面，河道左、右岸COD、氨氮質量濃度梯度較大，且存在明顯的混合層（高質量濃度梯度區）。 以C2 斷面為代表進行統計，混合層范圍距離左岸200～1 000 m，長約800 m。

圖10 C1 斷面COD、氨氮質量濃度橫向分布

圖11 C2 斷面COD、氨氮質量濃度橫向分布

圖12 C4 斷面COD、氨氮質量濃度橫向分布

數值計算結果表明，污染物橫向擴散梯度與斷面紊動能分布相關。 以C2 斷面COD 橫向質量濃度梯度為例（見圖13），在距離左岸150 ～400 m 范圍內，COD質量濃度梯度達到峰值，最大為0.030 mg／（L·m）；相應地，C2 斷面紊動能在距離左岸150 ～400 m 范圍內逐漸增大，最大為0.005 m2／s2，表明污染物質量濃度橫向擴散梯度與水體紊動能分布顯著相關。

圖13 C2 斷面COD 橫向質量濃度梯度與紊動能

（4）污染物垂向擴散。 數值模擬結果表明，汾河河口C1 斷面及下游C2 斷面污染物在垂直方向上發生充分摻混，未見明顯三維特性。 在C1 斷面，高質量濃度COD、氨氮均集中分布于河道左岸1 km范圍內（見圖14）；在匯合口下游約3 km 的C2 斷面，污染物已逐漸混合至河道中心位置（見圖15）。 污染物沿河道垂直方向（z方向）基本沒有發生梯度變化。

圖14 C1 斷面COD、氨氮質量濃度垂向分布

圖15 C2 斷面COD、氨氮質量濃度垂向分布

5 結論

黃河流域水資源開發過度河段污染物入河量和水域納污能力矛盾突出，37%的河流水域承擔了流域超91%的入河污染負荷，目前黃河干流水質總體能滿足水質目標要求，但部分支流污染嚴重。 采用大渦模擬模型耦合物質對流擴散方程，對汾河入黃河干流交匯處河段污染物的傳遞過程進行數值模擬，在實尺度河流污染物模擬技術上有所創新，可準確捕捉污染物擴散混合的三維效應，更加全面揭示污染物的傳播過程，相對于物理模型試驗和數值模擬等小尺度模型，其空間尺度雷諾數的量級更加符合實際情況，得到的成果具有較強的實用價值。

天然河道污染物縱向擴散規律易受不規則地形影響，水流流速分布不均勻，汾河入黃河干流交匯口橫斷面流速不均勻系數為0.973，水流紊動能最大值為0.036 m2／s2，在匯合口下游左岸形成了較大范圍的低流速區，引起了污染物的滯留。 交匯口河段渦量值為0～0.031／s，最大渦量值出現在交匯口下游河道突然縮窄處，該區域水流流速增大，加快了污染物的擴散。 研究河段屬于寬淺河道，污染物質量濃度分布不存在明顯的垂向梯度，但具有較強的橫向及縱向不均勻性，受地形影響污染物在汾河入黃口下游形成約11 km 長的污染帶向下游輸移。 計算區域河段寬深比大，污染物垂向混合擴散效應極小，污染物擴散以二維為主，濃度橫向擴散梯度與水體紊動能正相關，紊動能越大，污染物濃度梯度越大。