礦渣堆積對河網水質影響的數值模擬分析

李文鈺，侯精明，劉占衍，張 松，欒廣學，杜穎恩，韓占濤

（1.西安理工大學 西北旱區生態水利國家重點實驗室，陜西 西安 710048；2.河北省邯鄲水文勘測研究中心，河北 邯鄲 056011；3.四川水利職業技術學院，四川 成都 611830；4.河南智河工程技術有限公司，河南 鄭州 450003；5.生態環境部 土壤與農業農村生態環境監管技術中心，北京 100020）

隨著人類社會的發展，對自然資源的需求不斷提升，這一過程中礦產資源的開發利用造成大量礦渣堆積，威脅生態安全［1－3］。如富含多種礦產資源的文峪河流域，因早年廢棄礦渣未得到有效處理而致使其在河道堆積，上游來水或降水對其沖刷、淋溶和浸泡后產生了大量高濃度污染物，對下游河道水環境及居民飲用水安全造成了極大威脅。為遏止礦渣堆積對生態環境的不利影響［4－8］，提出高效精準的水環境治理方案，需要定量分析礦渣堆積后流域水系污染物的分布特征及遷移總量。

近年來國內外學者常使用的水環境模型有SWMM、SWAT、MIKE 模型等，如：曾曉嵐 等［9］使用SWMM 模型對滇池東岸花卉大棚種植區降雨徑流水質、水量進行了模擬，并分析了氮磷的單位面積負荷率及輸出特征；劉騫等［10］應用SWAT 模型構建了岷江流域分布式水文和污染負荷模型，并計算了各項減排措施及氣象驅動條件對流域水環境改善的貢獻率；劉晨輝等［11］采用MIKE 模型模擬了COD、BOD5、氨氮等在長江的遷移和擴散情況。以上模型對特定區域具有較好的模擬效果，但SWMM 等一維模型無法展示污染物的對流擴散過程；SWAT 等水文模型參數較多，建模困難，且根據其本底特征確定污染物輸移路徑后無法在模擬過程中依照水動力條件自動調整路徑，適用于大尺度流域水文模擬，無法對中小尺度流域的水動力和水質變化過程進行高分辨率模擬；MIKE 等水動力水質模型雖然可以進行高分辨率水環境模擬，但其模型重點在于水動力模塊，且對于百萬數量級以上網格，模型運行速度較慢。基于此，本文采用GPU 加速技術構建了小文峪河流域水環境數值模型，并充分利用現有數據模擬分析小文峪河流域污染物空間分布及污染物來源，以期為小文峪河流域水環境治理提供參考。

1 材料與方法

1.1 研究區概況

文峪河位于湖北省十堰市竹山縣西北部，上游主要支流有小文峪河和界嶺河，下游匯入堵河。文峪河流域屬于亞熱帶季風氣候區，雨量充沛，熱量充足。水量多寡取決于氣候條件，呈季節性變化。小文峪河流域地表水主要污染源為礦洞和礦渣堆，礦渣堆沿自然坡面、河道溝谷及近溝坡面溜坡放置，降水時水流從礦渣堆表面和內部流過，形成礦渣堆淋溶水，嚴重影響小文峪河的水環境。小文峪河從上游到下游主要支流有竹園溝、清石溝、孫家溝、寶泉寺及郭家溝（見圖1），各支流礦渣堆面積分別為8 402、5 997、951、16 638、20 824 m2，礦渣方量分別為50 423、19 679、3 377、96 788、84 500 m3。本研究選取超標污染物硫酸根離子和鐵離子為目標污染物，對研究區污染物時空分布進行高分辨率模擬并量化分析小文峪河污染物來源。

圖1 研究區礦渣堆分布及監測點布設

1.2 流量及污染物數據

本研究采用2021 年野外單場次監測數據，監測點15 個，見圖1。監測指標主要有硫酸根離子、鐵離子、流量等。監測時間為旱季，流量較為穩定，因此以該場次監測數據為模型率定數據。監測點主要分布在礦渣堆上下游及支流與主河道交匯處，因此以支流礦渣堆上游監測水質和流量數據為模型邊界條件，認為無礦渣堆積的支流上游來水污染物濃度為0。支流上游來水對部分礦渣進行沖刷與礦洞水混合后經徑流過程到達下游，發生復雜的物理變化，采用一級沖刷過程來綜合概化地表水流經礦渣堆并與礦洞水混合到達下游的污染物濃度變化規律。一級支流在向主河道匯流過程中會有二級支流的匯入，由于缺乏二級支流監測數據，因此采用一級衰減過程來描述污染物濃度下降過程。小文峪河主河道上游流量和水質邊界條件根據監測數據設置，即不考慮研究區域上游污染物空間分布及水量匯集特征。

1.3 數值模型及求解方法

模型的地表水動力模塊控制方程為平面二維淺水方程（簡稱SWEs）［12－14］，主要針對具有自由液面且以平面運動為主的水流，只考慮水平方向流速，忽略垂向流動，且忽略運動黏性項、紊流黏性項、風應力和科氏力。二維淺水方程的守恒格式可用如下矢量形式來表示：

式中：t為時間；qx、qy分別為x、y方向的單寬流量；h為水深；q為變量矢量，包括水深及x、y方向的單寬流量；F、G分別為x、y方向的通量矢量；g為重力加速度；S為源項，包括降水源項、底坡源項及摩阻源項；u、v分別為x、y方向的流速；i為入滲源項；zb為河床底面高程；Cf為謝才系數，Cf＝gn2／h1／3，其中n為曼寧系數。

水質模塊控制方程為對流擴散方程：

式中：C為污染物垂向平均濃度；Dx、Dy分別為x、y方向的擴散系數，因對流較強，故忽略擴散作用；qin為點源排放的流量強度；Cin為點源污染物垂向平均濃度；S為污染物濃度對時間的全導數dC／dt，即轉化項，如污染物生化反應、生物作用等。

僅考慮衰減與沖刷反應時，轉化項S可表示為

式中：K1為衰減系數，K2為沖刷系數。

采用Godunov 格式有限體積法求解二維淺水方程，選用二階MUSCL 方法進行變量重構［15－16］，通過HLLC 近似黎曼求解器求解控制單元界面上的質量通量與動量通量。為提高計算效率，模型結合CUDA 架構實現GPU 加速［17］，同時保持精確度和高穩定性。

1.4 污染物通量計算

污染物通量是指水體中的污染物在單位時間內通過某一斷面的總質量，是水文、地質、化學及生物等綜合作用的結果［18－20］。污染物通量計算公式為

式中：L為污染物通量；Q為平均流量；m為轉化系數，值為0.086 4［18］。

2 模型構建

地形數據的分辨率為2 m，共有269 437 個網格，模型計算采用開邊界，僅在支流監測點位置設置入流流量和污染物初始濃度。研究區處于山區，且河道長期有水，可認為河道水流無下滲；參照文獻［21－22］，根據當地實際情況，取曼寧系數為0.04。由于礦洞水僅有極少部分能匯入地表水，且流量較小，因此將其與各支流的發源地進行合并建模。研究區以礦渣堆為主的污染源釋放的污染物會進入地表水和地下水，再隨地表水和地下水向下游遷移，使下游地表水和地下水中污染物超過相應水質標準。由于礦渣堆積方量不同，因此采用GIS 技術確定各支流礦渣堆積的范圍，分別率定水流經過各支流礦渣堆及其上下游時污染物的沖刷系數和衰減系數，當地表水以一定流速通過礦渣堆時，模型將計算出各支流礦渣堆析出的污染物量。

為分析礦渣堆中硫酸根離子和鐵離子對流域水質的影響，采用基于GPU 加速技術的水動力水質耦合模型對污染物的遷移轉化過程進行模擬。為驗證模型模擬結果的準確性，根據郭家溝、寶泉寺、孫家溝、清石溝及竹園溝的水質監測數據，對模型進行參數率定，結果見表1。

表1 水質參數率定結果 10－3／s

3 結果分析

水質模擬結果見表2，可以看出各監測點污染物模擬結果相對誤差均在10%以內，表明所構建模型能較好地模擬小文峪河流域硫酸根離子和鐵離子在河網中的遷移轉化過程。

表2 硫酸根離子和鐵離子模擬結果和實測結果對比

3.1 污染物空間分布特征

研究區硫酸根離子和鐵離子空間分布見圖2、圖3，可以看出礦渣堆上游污染物質量濃度均較低，經過礦渣堆途中污染物質量濃度逐漸上升，在礦渣堆下游污染物質量濃度達到峰值。竹園溝、清石溝、孫家溝、寶泉寺、郭家溝5 條支流硫酸根離子質量濃度分別從26.0、23.5、134.0、19.5、20.3 mg／L 上升至222.0、121.1、170.0、453.0、409.0 mg／L，鐵離子質量濃度分別從0.13、0.07、6.94、0.06、0.15 mg／L 上升至10.60、28.90、10.50、83.10、59.60 mg／L。礦渣堆被沖刷過程中，各支流污染物質量濃度總體變化趨勢大致相同，在入流點到礦渣堆上游污染物質量濃度呈下降趨勢，在礦渣堆上游至礦渣堆下游污染物質量濃度呈上升趨勢。其中竹園溝上游礦渣堆積范圍較廣，水流到達該支流與主河道交匯處時，污染物質量濃度達到最大，之后受上游水流匯入、污染物衰減過程及河網本底特征影響，污染物質量濃度呈減小趨勢。清石溝礦渣堆面積較小，在水流經礦渣堆到達主河道之前，兩種污染物質量濃度經稀釋和底泥吸附等作用有所下降。支流孫家溝地表水的硫酸根離子和鐵離子質量濃度均小于上游來水的，水流匯合后對主河道污染物具有一定稀釋作用，使主河道污染物質量濃度下降；寶泉寺和郭家溝礦渣堆直達下游與主河道交匯處，故在交匯前硫酸根離子和鐵離子質量濃度均較大，且寶泉寺和郭家溝支流流量較大，因此水流進入主河道后，硫酸根離子和鐵離子質量濃度均發生了明顯變化，致使主河道下游污染物質量濃度增大。礦渣堆下游是河網污染物集中整治的關鍵區域。

圖2 硫酸根離子質量濃度空間分布

圖3 鐵離子質量濃度空間分布

3.2 水量及污染物通量分析

各支流對小文峪河干流水量貢獻不一，由大到小依次為寶泉寺、郭家溝、竹園溝、清石溝和孫家溝，占比分別為21.97%、21.21%、20.46%、18.94%、17.42%。整體水量差距相對不大，但由于各支流下游污染物質量濃度不同，因此各支流硫酸根離子和鐵離子負荷對主河道的貢獻率差距較大。各支流水流和日污染物通量見圖4。各支流輸入主河道的污染物通量由大到小依次為寶泉寺、郭家溝、竹園溝、清石溝、孫家溝，其中硫酸根離子通量分別為1.135、0.989、0.518、0.157、0.079 t／d，鐵離子通量分別為0.204、0.139、0.026、0.019、0.007 t／d。硫酸根離子和鐵離子入河通量在不同支流變化趨勢較為一致，最大值均出現在寶泉寺。受水量和礦渣堆積方量不同影響，污染物通量在各支流呈現不同變化，以硫酸根離子為例，其污染物入河通量的最大值出現在寶泉寺，最小值出現在孫家溝，最大值為最小值的14.4倍，而寶泉寺水量僅為孫家溝的1.26倍，表明污染物通量的差異主要與礦渣堆積方量有關。

圖4 污染物通量及水量

3.3 空間分布特征

各支流硫酸根離子、鐵離子入河通量大小順序均為寶泉寺＞郭家溝＞竹園溝＞清石溝＞孫家溝，可見入河污染物主要來源于寶泉寺和郭家溝，二者占比合計在70%以上，其他支流的占比相對較小。硫酸根離子污染負荷主要來自郭家溝和寶泉寺，其對小文峪河主河道的貢獻率分別為34.36%、39.42%；除郭家溝和寶泉寺外，竹園溝對河網中硫酸根離子的貢獻率也較大，為17.99%。在各支流中，寶泉寺對小文峪河主河道鐵離子的貢獻率最大，為51.64%，其次為郭家溝，為35.10%。由各支流對小文峪河主河道污染物的貢獻率可知，整體上，寶泉寺的污染負荷大于郭家溝的，其次為竹園溝、清石溝和孫家溝；流域內污染負荷主要來自下游支流寶泉寺和郭家溝，其產生的硫酸根離子、鐵離子污染負荷分別占整個研究區的73.78% 和86.74%。

各支流污染物通量及貢獻率不同的主要原因是，礦渣堆積方量不同，其中寶泉寺和郭家溝礦渣堆占地面積分別為16 638、20 824 m2，其他3 條支流礦渣堆占地面積之和為15 350 m2；寶泉寺和郭家溝流量較大，水流對礦渣堆的沖刷較為嚴重，導致污染物質量濃度增大，污染負荷增大。因此對小文峪河流域的水環境治理方案制定過程中，應首先考慮寶泉寺和郭家溝兩條支流的整治。

4 結論

以十堰市小文峪河流域為例，基于二維水動力水質耦合模型，充分利用各支流礦渣堆屬性、流域河網本底特征及地表水、礦洞水、淋溶水水質和水量監測數據等，構建了基于GPU 加速技術的礦山區河網水環境高分辨率數值模型，對小文峪河流域硫酸根離子和鐵離子的空間分布進行了可視化模擬，并量化分析了小文峪河污染負荷來源。模型模擬結果與污染物實測值相比，相對誤差小于10%，表明模型可以對礦山區河網水環境變化過程進行模擬。小文峪河主河道的水量貢獻率由大到小分別為寶泉寺、郭家溝、竹園溝、清石溝和孫家溝，占比分別為21.97%、21.21%、20.46%、18.94%、17.42%。硫酸根離子和鐵離子對小文峪河貢獻率由大到小均為寶泉寺＞郭家溝＞竹園溝＞清石溝＞孫家溝，其中寶泉寺和郭家溝硫酸根離子貢獻率分別為39.42%和34.36%，鐵離子貢獻率分別為51.64%和35.10%。

綜上，基于GPU 加速技術的水動力水質耦合模型對礦山區河網水環境進行模擬效果較好，在小文峪河流域水環境治理方案制定過程中，應首先考慮寶泉寺和郭家溝兩條支流的整治。