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

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

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

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

1 黃河干支流水質(zhì)狀況

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

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

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

圖3 2000—2018 年黃河支流水質(zhì)變化

2 研究區(qū)域水環(huán)境特性

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

圖4 汾河與黃河交匯處河道形態(tài)

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

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

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

3 數(shù)學(xué)模型構(gòu)建

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

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

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

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

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

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

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

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

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

4 計算結(jié)果分析

4.1 模擬工況設(shè)置

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

表1 模擬工況邊界條件

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

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

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

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

圖7 COD 和氨氮年均質(zhì)量濃度輸移云圖

4.3 交匯區(qū)污染物三維混合擴(kuò)散特性

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

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

式中：Cp為平均預(yù)測濃度；Cs（x，y）為坐標(biāo)（x，y）處的濃度。

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

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

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

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

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

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

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

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

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

數(shù)值計算結(jié)果表明，污染物橫向擴(kuò)散梯度與斷面紊動能分布相關(guān)。 以C2 斷面COD 橫向質(zhì)量濃度梯度為例（見圖13），在距離左岸150 ～400 m 范圍內(nèi)，COD質(zhì)量濃度梯度達(dá)到峰值，最大為0.030 mg／（L·m）；相應(yīng)地，C2 斷面紊動能在距離左岸150 ～400 m 范圍內(nèi)逐漸增大，最大為0.005 m2／s2，表明污染物質(zhì)量濃度橫向擴(kuò)散梯度與水體紊動能分布顯著相關(guān)。

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

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

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

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

5 結(jié)論

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

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