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

（1.西安理工大學(xué) 西北旱區(qū)生態(tài)水利國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710048；2.河北省邯鄲水文勘測(cè)研究中心，河北 邯鄲 056011；3.四川水利職業(yè)技術(shù)學(xué)院，四川 成都 611830；4.河南智河工程技術(shù)有限公司，河南 鄭州 450003；5.生態(tài)環(huán)境部 土壤與農(nóng)業(yè)農(nóng)村生態(tài)環(huán)境監(jiān)管技術(shù)中心，北京 100020）

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

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

1 材料與方法

1.1 研究區(qū)概況

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

圖1 研究區(qū)礦渣堆分布及監(jiān)測(cè)點(diǎn)布設(shè)

1.2 流量及污染物數(shù)據(jù)

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

1.3 數(shù)值模型及求解方法

模型的地表水動(dòng)力模塊控制方程為平面二維淺水方程（簡(jiǎn)稱SWEs）［12－14］，主要針對(duì)具有自由液面且以平面運(yùn)動(dòng)為主的水流，只考慮水平方向流速，忽略垂向流動(dòng)，且忽略運(yùn)動(dòng)黏性項(xiàng)、紊流黏性項(xiàng)、風(fēng)應(yīng)力和科氏力。二維淺水方程的守恒格式可用如下矢量形式來(lái)表示：

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

水質(zhì)模塊控制方程為對(duì)流擴(kuò)散方程：

式中：C為污染物垂向平均濃度；Dx、Dy分別為x、y方向的擴(kuò)散系數(shù)，因?qū)α鬏^強(qiáng)，故忽略擴(kuò)散作用；qin為點(diǎn)源排放的流量強(qiáng)度；Cin為點(diǎn)源污染物垂向平均濃度；S為污染物濃度對(duì)時(shí)間的全導(dǎo)數(shù)dC／dt，即轉(zhuǎn)化項(xiàng)，如污染物生化反應(yīng)、生物作用等。

僅考慮衰減與沖刷反應(yīng)時(shí)，轉(zhuǎn)化項(xiàng)S可表示為

式中：K1為衰減系數(shù)，K2為沖刷系數(shù)。

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

1.4 污染物通量計(jì)算

污染物通量是指水體中的污染物在單位時(shí)間內(nèi)通過(guò)某一斷面的總質(zhì)量，是水文、地質(zhì)、化學(xué)及生物等綜合作用的結(jié)果［18－20］。污染物通量計(jì)算公式為

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

2 模型構(gòu)建

地形數(shù)據(jù)的分辨率為2 m，共有269 437 個(gè)網(wǎng)格，模型計(jì)算采用開邊界，僅在支流監(jiān)測(cè)點(diǎn)位置設(shè)置入流流量和污染物初始濃度。研究區(qū)處于山區(qū)，且河道長(zhǎng)期有水，可認(rèn)為河道水流無(wú)下滲；參照文獻(xiàn)［21－22］，根據(jù)當(dāng)?shù)貙?shí)際情況，取曼寧系數(shù)為0.04。由于礦洞水僅有極少部分能匯入地表水，且流量較小，因此將其與各支流的發(fā)源地進(jìn)行合并建模。研究區(qū)以礦渣堆為主的污染源釋放的污染物會(huì)進(jìn)入地表水和地下水，再隨地表水和地下水向下游遷移，使下游地表水和地下水中污染物超過(guò)相應(yīng)水質(zhì)標(biāo)準(zhǔn)。由于礦渣堆積方量不同，因此采用GIS 技術(shù)確定各支流礦渣堆積的范圍，分別率定水流經(jīng)過(guò)各支流礦渣堆及其上下游時(shí)污染物的沖刷系數(shù)和衰減系數(shù)，當(dāng)?shù)乇硭砸欢魉偻ㄟ^(guò)礦渣堆時(shí)，模型將計(jì)算出各支流礦渣堆析出的污染物量。

為分析礦渣堆中硫酸根離子和鐵離子對(duì)流域水質(zhì)的影響，采用基于GPU 加速技術(shù)的水動(dòng)力水質(zhì)耦合模型對(duì)污染物的遷移轉(zhuǎn)化過(guò)程進(jìn)行模擬。為驗(yàn)證模型模擬結(jié)果的準(zhǔn)確性，根據(jù)郭家溝、寶泉寺、孫家溝、清石溝及竹園溝的水質(zhì)監(jiān)測(cè)數(shù)據(jù)，對(duì)模型進(jìn)行參數(shù)率定，結(jié)果見表1。

表1 水質(zhì)參數(shù)率定結(jié)果 10－3／s

3 結(jié)果分析

水質(zhì)模擬結(jié)果見表2，可以看出各監(jiān)測(cè)點(diǎn)污染物模擬結(jié)果相對(duì)誤差均在10%以內(nèi)，表明所構(gòu)建模型能較好地模擬小文峪河流域硫酸根離子和鐵離子在河網(wǎng)中的遷移轉(zhuǎn)化過(guò)程。

表2 硫酸根離子和鐵離子模擬結(jié)果和實(shí)測(cè)結(jié)果對(duì)比

3.1 污染物空間分布特征

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

圖2 硫酸根離子質(zhì)量濃度空間分布

圖3 鐵離子質(zhì)量濃度空間分布

3.2 水量及污染物通量分析

各支流對(duì)小文峪河干流水量貢獻(xiàn)不一，由大到小依次為寶泉寺、郭家溝、竹園溝、清石溝和孫家溝，占比分別為21.97%、21.21%、20.46%、18.94%、17.42%。整體水量差距相對(duì)不大，但由于各支流下游污染物質(zhì)量濃度不同，因此各支流硫酸根離子和鐵離子負(fù)荷對(duì)主河道的貢獻(xiàn)率差距較大。各支流水流和日污染物通量見圖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。硫酸根離子和鐵離子入河通量在不同支流變化趨勢(shì)較為一致，最大值均出現(xiàn)在寶泉寺。受水量和礦渣堆積方量不同影響，污染物通量在各支流呈現(xiàn)不同變化，以硫酸根離子為例，其污染物入河通量的最大值出現(xiàn)在寶泉寺，最小值出現(xiàn)在孫家溝，最大值為最小值的14.4倍，而寶泉寺水量?jī)H為孫家溝的1.26倍，表明污染物通量的差異主要與礦渣堆積方量有關(guān)。

圖4 污染物通量及水量

3.3 空間分布特征

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

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

4 結(jié)論

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

綜上，基于GPU 加速技術(shù)的水動(dòng)力水質(zhì)耦合模型對(duì)礦山區(qū)河網(wǎng)水環(huán)境進(jìn)行模擬效果較好，在小文峪河流域水環(huán)境治理方案制定過(guò)程中，應(yīng)首先考慮寶泉寺和郭家溝兩條支流的整治。