孫 杰，郝 喆，陳紅丹

(遼寧大學(xué) 環(huán)境學(xué)院，遼寧 沈陽 110036)

0 引 言

尾礦庫用于堆放礦山開采后產(chǎn)生的尾礦，是礦山開采中的重要設(shè)施[1]。它也是礦山開采中產(chǎn)生的人造污染源[2]，若管理不當(dāng)，可能會造成尾礦砂外泄、污染物擴(kuò)散等后果，嚴(yán)重威脅周圍環(huán)境、資源及下游人民的生命財產(chǎn)安全[3]。尾礦壩是尾礦庫中最主要的構(gòu)筑物，由初期壩和后期尾礦材料形成的壩體共同構(gòu)成[4]。尾礦廢水中含有重金屬污染物及其他有毒有害物質(zhì)，其經(jīng)過遷移作用會對周圍的水和土壤造成污染[5-7]，因此，尾礦庫引發(fā)的重金屬污染問題愈加受到關(guān)注。

目前，許多學(xué)者在尾礦庫污染物運(yùn)移的數(shù)值模擬方面已有一定的研究進(jìn)展。張建軍[1]等應(yīng)用GEO-SLOPE軟件中的SEEP/W和CTRAN/W模塊分析了尾礦砂各層滲透系數(shù)比、干灘長度等對滲流場的影響，并且模擬了考慮和不考慮吸附兩種情況下100天、1 000天壩體內(nèi)銅離子的遷移過程；陳紅丹[8]等模擬了降雨、蒸發(fā)與地下水耦合作用下的大型尾礦庫重金屬Cu2+的遷移過程；周澤超[9]等利用TOUGH2軟件根據(jù)建立的地下水流動模型與地下水污染物遷移模型，模擬計算出了各時刻銅離子在壩體中的分布情況；Yu Liren等[10]模擬研究了阿納馬市附近的索利莫斯河的流動與污染物的運(yùn)移；B.O.Tsydenov[11]通過數(shù)值模擬分析了某海域秋季熱壩發(fā)展過程中，污染物在不同風(fēng)況下的空間遷移規(guī)律及分布特點(diǎn)；Ben Abdelghani Farouk[12]等運(yùn)用Hydro-Geosphere軟件模擬了污染物在露天礦山中隨時間推移的遷移情況。上述研究多基于現(xiàn)有尾礦庫研究重金屬污染物的對流-彌散作用，未見在尾礦庫堆積過程中，尾礦壩加高對于重金屬遷移范圍的影響研究。本文以本鋼集團(tuán)歪頭山尾礦庫為例，利用GeoStudio軟件中的SEEP/W和CTRAN/W模塊建立尾礦庫數(shù)值模型，研究尾礦壩加高對于重金屬Cu2+在尾礦庫中遷移范圍的影響，為大型尾礦庫環(huán)境地質(zhì)評價和污染預(yù)測防治提供依據(jù)。

1 研究區(qū)域及方法

1.1 研究區(qū)域

研究區(qū)域位于遼寧省本溪市溪湖區(qū)，該尾礦庫屬于二等尾礦庫。庫區(qū)占地面積1.62 km2，現(xiàn)庫容量約為1億m3。尾礦庫初期壩由透水堆石壩構(gòu)成，初期壩壩頂標(biāo)高為191.9 m，壩底標(biāo)高180.0 m，相對高度11.9 m，頂寬4.0 m，壩長410.0 m，內(nèi)坡比1∶2，外坡比1∶2.5。尾礦庫的主壩尾礦堆積，壩頂標(biāo)高為280.0 m，壩庫壩底標(biāo)高180 m，主壩高度差約為100 m，壩頂壩長約為1 600 m。

1.2 研究方法

Geostudio系統(tǒng)軟件是由加拿大巖土軟件GEO-SLOPE公司開發(fā)的面向巖土、采礦、交通、水利、地質(zhì)、環(huán)境工程等領(lǐng)域的仿真分析軟件，是全球最知名的巖土工程分析軟件之一。它包含以下八種專業(yè)分析模塊：SLOPE/W(邊坡穩(wěn)定性分析)、SEEP/W(地下水滲流分析)、SIGMA/W(巖土應(yīng)力變形分析)、QUAKE/W(地震動力響應(yīng)分析)、TEMP/W(地?zé)岱治?、CTRAN/W(地下水污染物傳輸分析)、VADOSE/W(綜合滲流蒸發(fā)區(qū)和土壤表層分析)。當(dāng)用戶需要進(jìn)行多模塊工作時，只需將模型的邊界條件和材料參數(shù)等定義在一個模塊中，使用其他模塊時不需重復(fù)輸入相關(guān)參數(shù)，因此該軟件建模高效、操作簡單[13]。本研究中所使用的Geostudio版本為2012版。

SEEP/W模塊主要用于分析多孔滲水材料，它可以分析飽和穩(wěn)態(tài)及飽和-不飽和瞬態(tài)問題，幾乎可以對所有地下水問題進(jìn)行建模分析，其計算原理是通過內(nèi)插函數(shù)來表證該單元上每一節(jié)點(diǎn)處的水頭值，單元內(nèi)的水頭分布用線性差值函數(shù)來近似。CTRAN/W模塊主要用于分析污染物通過土層和巖石等滲流介質(zhì)外界環(huán)境中的水體，經(jīng)過地面和地下水的非飽和區(qū)域進(jìn)入地下水體的問題，其計算原理是伽遼金有限元法，該方法從剩余加權(quán)法出發(fā)對對流-彌散方程進(jìn)行數(shù)值離散，是求解微分方程近似解的一種方法。CTRAN/W和SEEP/W相結(jié)合，可以用來分析污染物的運(yùn)移。SEEP/W模塊計算水的流速、水體積含量和水流量，CTRAN/W模塊則用這些參數(shù)計算污染物的遷移。SEEP/W和CTRAN/W模塊均屬于有限元分析軟件。

運(yùn)用GeoStudio中的SEEP/W和CTRAN/W兩個模塊開展尾礦壩加高對重金屬污染物Cu2+在尾礦庫中遷移規(guī)律影響的數(shù)值模擬研究。在SEEP/W模塊中設(shè)置尾礦庫模型材料參數(shù)、邊界條件等參數(shù)，通過穩(wěn)態(tài)滲流計算得出尾礦庫孔隙水壓力等勢圖、體積含水率等結(jié)果；利用CTRAN/W模塊建立對流-彌散和粒子示蹤模型，通過計算重金屬Cu2+的對流-彌散過程，分析對流-彌散過程計算出的重金屬濃度等勢圖、粒子路徑遷移圖等結(jié)果，得到Cu2+在尾礦庫中遷移的情況，最終通過分析得出尾礦壩加高過程中Cu2+的遷移規(guī)律。

2 基礎(chǔ)原理及模型建立

2.1 Darcy定律

Darcy定律是滲流理論中的基本定律[14]，是形容飽和土中水的滲流速度與水力坡降之間的聯(lián)系的。其最初是從飽和土中得到的，同樣也適用于非飽和土，不同之處在于滲流系數(shù)是否為常量[15]。表達(dá)式為：

Q=ki

(1)

式(1)中：Q——流速，m/d；k——滲流系數(shù)，m/d；i——滲透坡降，無量綱。

2.2 滲流基本方程

滲流計算方程[16]：

(2)

式(2)中：H——總水頭，m；kx——x方向的滲流系數(shù)，m/d；ky——y方向的滲流系數(shù)，m/d；Q——邊界流量，m3/d；θ——體積水容量，L/m3；t——時間，d。

2.3 對流-彌散方程

污染物的遷移機(jī)理受到對流-彌散、吸附-沉淀、解吸、絡(luò)合、生物降解等作用。根據(jù)尾礦庫各參數(shù)性質(zhì)，本文僅考慮重金屬的對流-彌散作用。對流-彌散方程為[17]：

(3)

式(3)中：Dx——橫向彌散系數(shù)，m2/d；Dy——縱向彌散系數(shù)，m2/d；C——污染物濃度，g/m3；ux——橫向滲流速度，m/d；uy——縱向滲流速度，m/d。

2.4 模型建立

本文以歪頭山尾礦庫為例，模擬在尾礦壩加高過程中重金屬Cu2+的遷移過程，通過對比在不同尾礦壩壩高下重金屬的遷移范圍及其濃度分布情況，研究尾礦壩加高對重金屬污染物遷移的影響規(guī)律。如圖2所示，本次研究共選取三種不同壩高，工況1、工況2、工況3的頂部高程分別為236、263、280 m，三個概化計算模型剖面的地面高程均為180 m，計算區(qū)域的水平距離均為595 m。模型網(wǎng)格剖分為三角形或四邊形，在模擬中已做網(wǎng)格無關(guān)化檢驗，網(wǎng)格尺寸大小的劃分對重金屬遷移的影響可忽略不計，本文尾礦庫模型全局單元尺寸取5 m，工況1、工況2、工況3剖分單元分別為1 212個、1 450個、1 545個，節(jié)點(diǎn)分別為1 320個、1 551個、1 644個。

圖2 尾礦庫概化計算模型Fig.2 Generalized calculation model of the tailings pond

2.5 模型邊界條件

在進(jìn)行穩(wěn)態(tài)滲流、對流-彌散模擬時，設(shè)置排水邊界、零流量邊界、上游水頭、下游水頭等水力邊界條件和污染物邊界條件。如圖2所示，在穩(wěn)態(tài)滲流模擬中，AB為零流量邊界，BC為上游水頭，DE為下游水頭，EF為排水邊界。在進(jìn)行對流-彌散模擬時，BC為重金屬Cu2+污染物邊界，初始濃度為50 g/m3。

2.6 材料參數(shù)

在概化計算模型中，包含尾細(xì)砂、尾粉質(zhì)粘土、基巖等三種不同土層，其中前兩層為飽和/不飽和模型，基巖為飽和模型，各介質(zhì)水文地質(zhì)參數(shù)如表1所示，土水特征曲線及各土層滲透系數(shù)-基質(zhì)吸力曲線詳見圖3和圖4。

表1 各介質(zhì)水文地質(zhì)參數(shù)

圖3 各介質(zhì)土水特征曲線Fig.3 Soil and water characteristic curves of each medium

圖4 各介質(zhì)滲流系數(shù)-基質(zhì)吸力特征曲線Fig.4 Seepage coefficient-matric suction characteristic curve of each medium

3 數(shù)值模擬計算結(jié)果分析

3.1 滲流分析

圖5為不同壩高下的尾礦庫孔隙水壓力等勢圖，圖中所示剖面1和剖面2交點(diǎn)高程坐標(biāo)分別為(211.93 m，240.31 m)、(146.35 m，223.01 m)。圖6為兩個剖面在工況1高度、工況2高度和工況3高度下的垂直方向孔隙水壓力等勢分布情況。

根據(jù)圖5和圖6可知，在同一尾礦庫模型中，垂直方向上孔隙水壓力隨著尾礦庫深度的增加而增大，尾礦庫底部的孔隙水壓力最大，尾礦庫頂部的孔隙水壓力最小。且孔隙水壓力隨高程變化的曲線大致呈一條直線，所以孔隙水壓力隨高程變化的關(guān)系為一次函數(shù)關(guān)系，經(jīng)過數(shù)據(jù)分析可得函數(shù)斜率為-9.77，即高程每增加1 m，孔隙水壓力減少9.77 kPa。

圖5 穩(wěn)態(tài)滲流作用下尾礦庫孔隙水壓力等勢分布圖Fig.5 Equipotential distribution diagrams of pore water pressure of the tailing pond under steady seepage

圖6 垂直剖面孔隙水壓力分布Fig.6 Pressure distribution of pore water in vertical section

通過圖6可知，在剖面1處，當(dāng)高程為199.27 m 時，工況1、工況2、工況3的孔隙水壓力分別為29.43、180.69、244.08 kPa；當(dāng)高程為240.20 m 時，工況1、工況2、工況3的孔隙水壓力分別為-378.70、-217.80、-148.49 kPa。在剖面2處，當(dāng)高程為169.97 m時，工況1、工況2、工況3的孔隙水壓力分別為251.81、362.86、408.23 kPa；當(dāng)高程為223.00 m時，工況1、工況2、工況3的孔隙水壓力分別為-272.86、-148.59、-96.36 kPa。由此可知，孔隙水壓力會受到尾礦壩高度的影響，當(dāng)尾礦壩高度增加時，相同位置的孔隙水壓力增大。

3.2 重金屬Cu2+遷移規(guī)律分析

根據(jù)實(shí)際水質(zhì)檢測結(jié)果，在不考慮吸附和離子間相互作用的情況下，研究區(qū)域的重金屬Cu2+初始濃度為50 g/m3。通過對尾礦庫進(jìn)行對流-彌散和粒子示蹤模型計算，分析研究尾礦壩加高對于重金屬污染物的遷移規(guī)律。根據(jù)GeoStudio數(shù)值模擬2 920 d計算結(jié)果，得出重金屬Cu2+在不同壩高下的濃度分布情況和粒子遷移路徑，如圖7和8所示。圖7中粉色陰影表示Cu2+濃度高于40 g/m3，藍(lán)色陰影表示Cu2+濃度低于10 g/m3。

通過觀察圖7可以看出，在同一尾礦庫模型中，污染物濃度沿初始濃度邊界向初期壩方向逐漸降低。通過對比發(fā)現(xiàn)，不同壩高的尾礦庫重金屬污染物Cu2+遷移的質(zhì)量濃度等勢分布不同，尾礦壩加高會對污染物遷移產(chǎn)生影響。在剖面1和剖面2處可以看出，當(dāng)尾礦庫深度相同時，尾礦壩越高，重金屬Cu2+濃度越大。由圖8粒子遷移路徑可知，重金屬Cu2+是隨地下水沿尾礦壩方向遷移的，且隨著尾礦壩的加高，污染物的遷移距離越大。

圖7 重金屬Cu2+遷移2 920 d濃度等勢圖Fig.7 Concentration equipotential diagrams of heavy metal Cu2+ after migration for 2 920 days

根據(jù)圖9可知，在同一尾礦庫模型中，Cu2+的質(zhì)量濃度隨著高程的增加而減小，尾礦庫底部濃度最高，尾礦庫表層濃度最低且接近于0。在剖面1處，當(dāng)高程為169.97 m時，工況1、工況2、工況3的Cu2+質(zhì)量濃度分別為6.59、28.65、35.14 g/m3；當(dāng)高程為199.25 m時，工況1、工況2、工況3的Cu2+質(zhì)量濃度分別為3.63、26.27、33.97 g/m3。在剖面2處，當(dāng)高程為169.97 m時，工況1、工況2、工況3的Cu2+質(zhì)量濃度分別為2.72、20.72、27.83 g/m3；當(dāng)高程為200.90 m時，工況1、工況2、工況3的Cu2+質(zhì)量濃度分別為0.59、16.20、24.08 g/m3。由此可見，在同一深度，經(jīng)過相同的時間，尾礦壩加高會使重金屬污染物的質(zhì)量濃度增大。通過圖10可知，隨著時間的推移，在同一高程位置處，污染物Cu2+濃度呈上升趨勢。

圖8 重金屬Cu2+遷移2 920 d路徑分布Fig.8 Path distribution of heavy metal Cu2+ after migration for 2 920 days

圖9 垂直剖面重金屬Cu2+濃度分布(2 920 d)Fig.9 Concentration distribution of heavy metal Cu2+ in vertical profile(2 920 days)

圖10 不同遷移時間下Cu2+濃度分布(剖面1)Fig.10 Concentration distribution of Cu2+ after migration for different time(profile1)

圖11 各粒子移距變化(2 920 d)Fig.11 Displacement variation of each particle(2 920 days)

由圖11和圖12可以看出，尾礦壩加高會使污染物顆粒的遷移距離和平均遷移速度增大，且尾礦壩越高，顆粒的遷移距離和平均遷移速度越大。當(dāng)尾礦壩加高到工況2時，1～6號顆粒在2 920 d 內(nèi)的移動距離和平均遷移速度比工況1時分別增加了190.35%、84.14%、73.29%、74.95%、75.64%、76.21%；當(dāng)尾礦壩加高到工況3時，1～6號顆粒在2 920 d內(nèi)的移動距離和平均遷移速度比工況2時分別增加了57.92%、21.85%、23.57%、24.10%、24.64%、25.20%。

圖12 各粒子平均遷移速度變化(2 920 d)Fig.12 Variation of average migration velocity ofeach particle(2 920 days)

4 結(jié) 論

通過以歪頭山尾礦庫為例進(jìn)行數(shù)值模擬分析，利用GeoStudio軟件中的SEEP/W和CTRAN/W兩個模塊建立數(shù)值模型，研究尾礦壩加高對于重金屬Cu2+在尾礦庫中遷移范圍的影響，得出以下結(jié)論：

(1)在同一尾礦庫模型中，垂直方向上孔隙水壓力隨著尾礦庫深度的增加而增大，深度每增加1 m，孔隙水壓力增加9.77 kPa，尾礦庫底部的孔隙水壓力最大，尾礦庫頂部的孔隙水壓力最小。此外，孔隙水壓力會受到尾礦壩高度的影響，當(dāng)尾礦壩加高時，相同位置的孔隙水壓力增大。

(2)在同一尾礦庫模型中，重金屬污染物濃度沿初始濃度邊界向初期壩方向逐漸降低。當(dāng)尾礦庫深度相同時，尾礦壩加高會使重金屬Cu2+濃度增大。重金屬Cu2+是隨地下水沿尾礦壩方向遷移的，且隨著尾礦壩的加高，污染物的遷移距離越大。

(3)在同一尾礦庫模型中，Cu2+濃度隨著尾礦庫的深度增加而增大，尾礦庫底部濃度最高，尾礦庫表層重金屬濃度最低且接近于0，且隨著時間的推移，尾礦庫同一深度處，污染物Cu2+濃度呈上升趨勢。重金屬污染物Cu2+經(jīng)過相同時間的遷移下，尾礦壩加高會使重金屬污染物的質(zhì)量濃度增大。

(4)尾礦壩加高會使污染物顆粒的遷移距離和平均遷移速度增大，且尾礦壩越高，顆粒的遷移距離和平均遷移速度越大。