李文鈺,侯精明,劉占衍,張 松,欒廣學(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)境治理提供參考。
文峪河位于湖北省十堰市竹山縣西北部,上游主要支流有小文峪河和界嶺河,下游匯入堵河。文峪河流域?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è)
本研究采用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ū)域上游污染物空間分布及水量匯集特征。
模型的地表水動(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)定性。
污染物通量是指水體中的污染物在單位時(shí)間內(nèi)通過(guò)某一斷面的總質(zhì)量,是水文、地質(zhì)、化學(xué)及生物等綜合作用的結(jié)果[18-20]。污染物通量計(jì)算公式為
式中:L為污染物通量;Q為平均流量;m為轉(zhuǎn)化系數(shù),值為0.086 4[18]。
地形數(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
水質(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ì)比
研究區(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ì)量濃度空間分布
各支流對(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 污染物通量及水量
各支流硫酸根離子、鐵離子入河通量大小順序均為寶泉寺>郭家溝>竹園溝>清石溝>孫家溝,可見入河污染物主要來(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)首先考慮寶泉寺和郭家溝兩條支流的整治。
以十堰市小文峪河流域?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)首先考慮寶泉寺和郭家溝兩條支流的整治。