鉛鋅礦廢石淋濾液在地下水潛水含水層中的遷移模擬研究
——基于有限差分法

伍陽，黃鑫，張靜，楊利均

(1.四川景星環(huán)境科技有限公司,四川成都 610037；2.四川農(nóng)業(yè)大學(xué)，四川成都 611130)

地下水埋藏在土壤、巖石的孔隙、裂隙和溶隙中，以各種不同形式的形態(tài)存在[1]，全世界的地下水估計(jì)總量多達(dá)1.5億立方公里，幾乎占地球總水量的十分之一。目前，我國(guó)約有70%人口以地下水為主要飲用水源，全國(guó)95%以上的農(nóng)村人飲用地下水，全國(guó)40%的耕地使用地下水灌溉[2-3]。它為全球的絕大部分人口提供著水源且支持著農(nóng)業(yè)和工業(yè)活動(dòng)，是一種十分重要的自然資源。然而地下水水質(zhì)在水環(huán)境中十分脆弱，一旦受到污染，不僅危害時(shí)間長(zhǎng)，并且治理難度很大，因此對(duì)地下水污染的預(yù)防重要性遠(yuǎn)大于治理。

金屬礦產(chǎn)的開發(fā)帶來巨大經(jīng)濟(jì)與能源效益，但也會(huì)出現(xiàn)一系列對(duì)環(huán)境的負(fù)面影響。其中，廢石場(chǎng)淋濾液對(duì)周圍地下水環(huán)境的影響尤為突出。鉛鋅尾礦在降雨作用下形成含有重金屬元素的淋濾液，下滲到地下水環(huán)境中隨著地下水遷移。它在污染地下水環(huán)境的同時(shí)，帶來飲用水源污染風(fēng)險(xiǎn)和對(duì)居民健康的危害。通過數(shù)值模擬手段定量評(píng)估污染物在地下水環(huán)境中的遷移規(guī)律，對(duì)加強(qiáng)礦區(qū)監(jiān)測(cè)管理，以及周圍水環(huán)境的保護(hù)有著重要意義。盧曉華[4]通過建立地下水溶質(zhì)運(yùn)移模型，對(duì)廢水處理車間發(fā)生泄漏情況下地下水重金屬污染情景進(jìn)行了數(shù)值模擬與環(huán)境影響預(yù)測(cè)評(píng)價(jià)，對(duì)環(huán)保部門從源頭上防治和管理潛在的地下水重金屬污染、保護(hù)地下水環(huán)境具有重要的意義；馮書順等[5]采用數(shù)值方法對(duì)礦坑井涌水量進(jìn)行了預(yù)測(cè)和模擬；鄧鼎興[6]采用有限元法對(duì)礦山地下水污染的分區(qū)治理效果進(jìn)行了數(shù)值模擬檢驗(yàn)，結(jié)果表明數(shù)值模擬預(yù)測(cè)結(jié)果與實(shí)際監(jiān)測(cè)值較為接近，可為類似礦山地下水污染治理提供借鑒。憑借數(shù)值方法，在研究二維以及三維的地下水滲流和污染物的基本遷移的基礎(chǔ)之上，國(guó)內(nèi)外還出現(xiàn)了對(duì)于考慮污染物的化學(xué)反應(yīng)、生物降解、吸附、衰減等多過程的污染物遷移轉(zhuǎn)化研究[7-9]，以及地下水與地表水耦合[10]、地面沉降與地下水滲流耦合等多場(chǎng)耦合模型[11]。

本文以西部某鉛鋅礦為研究區(qū)域，將研究區(qū)域合理概化后，建立地下水水流模型與溶質(zhì)運(yùn)移模型。以有限差分方法上游加權(quán)法為基礎(chǔ)，應(yīng)用MATLAB軟件模擬了廢石場(chǎng)在非正常情況下特征污染物Pb離子的遷移結(jié)果。模擬結(jié)果可以為庫區(qū)地下水環(huán)境評(píng)價(jià)和治理提供科學(xué)依據(jù)。

1 研究區(qū)域概況

鉛鋅礦所在的東、北及西側(cè)都以山脊為界，南側(cè)以斯弄多支溝為界，地形起伏，總體呈現(xiàn)西高東低，從西北到東南向海拔高度逐漸減小。該地區(qū)受西南支氣流、西南季風(fēng)等大氣環(huán)流的控制而具有高原亞寒帶半干旱季風(fēng)氣候的特點(diǎn)，具有日照充足，空氣干燥稀薄，多大風(fēng)，夏秋多雨，干濕季節(jié)分明的主要?dú)夂蛱卣?。最大年降水量?05.7mm，年平均降水量為352.38mm，最大日降雨量24.2mm，降水主要集中在6～9月，占全年降水的87%；年蒸發(fā)量1938.4mm。

鉛鋅礦位于該區(qū)域支溝內(nèi)，礦區(qū)地形坡度大，切割較強(qiáng)烈，溝谷發(fā)育且縱坡降大，第四系覆蓋嚴(yán)重，水文地質(zhì)邊界條件簡(jiǎn)單，構(gòu)成晶屑凝灰?guī)r、角礫凝灰?guī)r、溶結(jié)凝灰?guī)r與巖漿巖裂隙含水為主的水文地質(zhì)單元，屬以裂隙含水層為主的礦床。松散堆積物呈條帶狀分布在溝谷平坦部位，形成匯水谷地；溝谷斜坡上部基巖裸露，中、下部覆蓋殘坡積層。巖石類型以晶屑凝灰?guī)r、角礫凝灰?guī)r、溶結(jié)凝灰?guī)r為主，少量地帶出露花崗斑巖、閃長(zhǎng)玢巖。因此，區(qū)內(nèi)地下水主要賦存在基巖裂隙中，溝谷為其主要排泄區(qū)，地下水與地表水分水嶺基本一致，地下水徑流途徑短，分散排泄，水量貧乏。地表水補(bǔ)給源以大氣降水和冰雪融水為主，降水多沿地表徑流，少量沿表土、基巖裂隙入滲，最終在排泄基準(zhǔn)面溢出。

主要水化學(xué)類型包括Ca-Mg-HCO3，pH值介于8.19～8.38之間，地下水礦化度介于414.45～422.62mg/L，屬中性偏弱堿性，質(zhì)軟到低礦化淡水，水質(zhì)狀況良好。

2 滲流與溶質(zhì)遷移數(shù)學(xué)模型的建立

2.1 水文地質(zhì)概化模型

圖1 水文模型概念圖Fig.1 Concept of hydrological model

研究區(qū)內(nèi)，在天然情況下地下水由西北向東南流動(dòng)，與河流1平行的兩側(cè)分水嶺的垂向上幾乎無地下水交換，故兩側(cè)山脊線可概化為零流量邊界。西北處的為海拔較高的地形起伏較大的區(qū)域，故選取西北的天然分水嶺作為隔水邊界線。河流1與河流2所處第四系松散層，加之居民處的水井開采對(duì)邊界水位幾乎無影響，故概化為定水頭邊界。研究區(qū)頂部常年接受降水補(bǔ)給，可概化為潛水面邊界。研究區(qū)底部主要為為礦體底板，巖層的透水性弱，故概化為隔水邊界[12]。

2.2 潛水含水層滲流控制方程

因?yàn)楸镜貐^(qū)巖層傾角平緩，可以把地下水當(dāng)作平面二維問題。此外，由于受地形起伏控制研究區(qū)地下水流場(chǎng)形態(tài)較為復(fù)雜，地下水以泉水及溢出面形式出露，為此本次計(jì)算采用穩(wěn)定流建模。根據(jù)達(dá)西定律，對(duì)其建立起二維穩(wěn)定流形式的控制方程[13]：

其中Kx為在橫向上的滲透系數(shù)，Ky為在縱向上的滲透系數(shù)，其單位為m/d。

2.3 污染物遷移控制方程

污染物的遷移主要受到兩方面的驅(qū)動(dòng)：一是濃度梯度引起的溶質(zhì)分子的擴(kuò)散，其擴(kuò)散通量與濃度差的關(guān)系滿足費(fèi)克定律；二是由于地下水流動(dòng)引起的機(jī)械彌散。溶質(zhì)運(yùn)移控制方程[14]：

其中θ為含水層有效孔隙率；Dx為橫向水動(dòng)力彌散系數(shù)；Dy為縱向水動(dòng)力彌散系數(shù)；Rd為阻滯因子；Vx代表在橫向上的地下水流速；Vy代表在縱向上的地下水流速。流速Vx或Vy由達(dá)西定律得到：

水動(dòng)力彌散系數(shù)包括機(jī)械彌散系數(shù)與分子擴(kuò)散系數(shù)：

其中V為地下水實(shí)際流速，Vx與Vy分別為地下水速度在橫向和縱向上的速度分量，Dm為分子擴(kuò)散系數(shù)。

2.4 初始條件與邊界條件

在本研究中，地下水滲流場(chǎng)模型中海需要確定的邊界條件有第一類邊界與第二類邊界。

地下水運(yùn)移的定解條件：

（1）第一類邊界條件：地下水滲流方程中的第一類邊界為定水頭邊界，本研究區(qū)域中的河流可視為定水頭邊界：

（2）第二類邊界條件：地下水滲流方程中第二類邊界為定流量邊界，本研究區(qū)域中的分水嶺為隔水邊界，即零流量邊界:

溶質(zhì)運(yùn)移的定解條件包括了初始條件和邊界條件：初始條件：

邊界條件：

（3）第一類邊界條件：廢石浸出液入滲為連續(xù)注入式，地下水污染溶質(zhì)遷移模擬公式參考《多孔介質(zhì)滲流與污染物遷移數(shù)學(xué)模型》（仵彥卿編著）一書中平面連續(xù)注入示蹤劑，一維流動(dòng)二維彌散模型，污染源連續(xù)注入含水層，污染物濃度在各個(gè)時(shí)刻均為常數(shù)：

（4）第二類邊界條件：污染物溶質(zhì)運(yùn)移的第二類邊界為定擴(kuò)散通量邊界，此邊界與地下水滲流的第二類邊界線重合，則此邊界為零擴(kuò)散通量邊界：

3 數(shù)值模擬及結(jié)果

3.1 有限差分?jǐn)?shù)值解法

有限差分法是一種將微分方程的控制方程及對(duì)應(yīng)定界條件離散，利用迭代法求解出微分方程數(shù)值解的方法。地下水流動(dòng)模型與污染物運(yùn)移模型是耦合的系統(tǒng)模型，但在穩(wěn)定流的情況下，可先求出穩(wěn)定流條件下的流場(chǎng)分布，再根據(jù)研究區(qū)域內(nèi)的各個(gè)點(diǎn)水頭求出速度場(chǎng)，最后根據(jù)離散化后的溶質(zhì)運(yùn)移差分方程迭代求出污染物隨時(shí)空變化的濃度值。

本文研究區(qū)域?yàn)楦鶕?jù)自然邊界線劃分的不規(guī)則區(qū)域。采用矩形網(wǎng)格剖分方法，共剖分了12334個(gè)節(jié)點(diǎn)，每個(gè)矩形網(wǎng)格大小為30×30m2。由于研究區(qū)域被劃分為正方形網(wǎng)格且含水層各向同性，即Δx=Δy，Kx=Ky對(duì)于地下水滲流方程，其離散后差分格式為[14-15]：

采用上游加權(quán)法，得到地下水溶質(zhì)運(yùn)移方程的顯式差分方程為：

3.2 定解條件的確定

如圖1所示，地下水水流模型的隔水邊界在整個(gè)研究區(qū)域的北側(cè)、西側(cè)和南側(cè)，定水頭邊界為內(nèi)部河流與東側(cè)邊界處河流邊界。地下水水流模型的隔水邊界為污染物溶質(zhì)運(yùn)移的零通量邊界。根據(jù)2016年地勘資料對(duì)研究區(qū)域的污染物檢測(cè)值作為初始濃度即CO=0.02mg/L。

3.3 參數(shù)選取

在服務(wù)期滿(11a)前，認(rèn)為污染源為連續(xù)注入式，即淋濾液通過地表和包氣帶到達(dá)地下水系統(tǒng)后的濃度為常數(shù)。其源濃度Cq計(jì)算按公式計(jì)算[14]：

其中sourceC為污染物下滲漏量[14]:

其中Cp為污染物實(shí)驗(yàn)最大浸出液濃度(mg/L)，α降水入滲補(bǔ)給系數(shù)，A為廢石場(chǎng)滲水面積(m2)，I為降水量(mm)。

廢石場(chǎng)中的廢水滲透穿過庫底后，進(jìn)入第四系的潛水含水層，漂石含量小于2%，卵石含量40%～55%，卵石最大直徑10cm；混雜粘土與砂礫，不均勻。根據(jù)項(xiàng)目水文地質(zhì)勘查報(bào)告，評(píng)價(jià)范圍內(nèi)坡積層孔隙水的滲透系數(shù)為5.4m/d。廢石場(chǎng)所在區(qū)域的潛水含水層都是第四系坡積物、洪積物或沖積物，有效孔隙度取值0.10。采取類比的方法，確定含水層的縱向彌散度為0.3(m2/d)。根據(jù)經(jīng)驗(yàn)公式取橫向彌散度為縱向彌散度的0.1倍，即橫向彌散度為0.03(m2/d)[16]。將研究區(qū)域已有水文資料中的參數(shù)作為初始值賦給模型，通過實(shí)測(cè)的地下水水文位數(shù)據(jù)與模擬的水文數(shù)據(jù)進(jìn)行對(duì)比驗(yàn)證，采用經(jīng)過識(shí)別與驗(yàn)證水位后的參數(shù)作為模型參數(shù)來預(yù)測(cè)地下水流場(chǎng)與溶質(zhì)運(yùn)移。驗(yàn)證后的參數(shù)取值見表1。

表1 模擬水文地質(zhì)參數(shù)及取值Tab.1 Chemical characteristics of the tested swine manure before and after leaching

3.4 結(jié)果分析

3.4.1 地下水流場(chǎng)分析

以海拔高度表示研究區(qū)域內(nèi)各點(diǎn)的水位高度，穩(wěn)定流情況下的地下水的水位等值線圖見圖2。

研究區(qū)的地下水流場(chǎng)受該地區(qū)的地形控制，一致地呈現(xiàn)出西高東低，由西北向東南遞減的規(guī)律。內(nèi)部河流的前半段為河流溝谷，地下水除自西北向東南向流動(dòng)，還會(huì)向內(nèi)部河流形成的溝谷匯集流動(dòng)，實(shí)現(xiàn)地下水補(bǔ)給地表河流。內(nèi)部河流后半段的沖積平層，即圖2由紅色虛線圍成三角區(qū)域，形成了一個(gè)微小的分水嶺，它的分流作用使水流不僅僅沿著內(nèi)部河流向邊界河流匯集。此外，邊界處的河流流向?yàn)樽阅舷虮?，匯集內(nèi)部河流后向下游流動(dòng)，在研究區(qū)域內(nèi)水流最終以邊界處的河流下游為出口。

圖2 水頭等值線圖Fig.2 Contour line of water head

研究區(qū)域內(nèi)的水力梯度不大，變化不劇烈，但由于所處區(qū)域的巖層滲透系數(shù)較大，會(huì)由于較強(qiáng)的對(duì)流作用使得污染物能夠較為快速地隨著地下水遷移。另一方面，水流主要流向邊界河流，最終沿邊界河流的下游方向遷移，使得該區(qū)域地下水污染物的不易累積。

3.4.2 污染物遷移規(guī)律分析

分別取在1000天、11年（服務(wù)期滿）、16年（服務(wù)期滿后5年）、30年后這幾個(gè)時(shí)間節(jié)點(diǎn)，得到各時(shí)間節(jié)點(diǎn)的污染物濃度分布見圖3。

在服務(wù)期滿前，廢石場(chǎng)淋濾液被簡(jiǎn)化為濃度恒定的污染源。污染物在地下水系統(tǒng)中主要受對(duì)流作用遷移，沿著地下水主要流動(dòng)方向由西北向東南方向遷移，同時(shí)具有向內(nèi)部河流形成的河谷匯集的趨勢(shì)。根據(jù)對(duì)流場(chǎng)的分析可知，由于靠近邊界河流處的扇形沖積平層的分水嶺作用，污染物在隨著地下水遷移到內(nèi)部河流下游時(shí)將被分散到河流兩側(cè)。同時(shí)在彌散作用下，在河流2的上游地帶部分污染物擴(kuò)散到了河流2的對(duì)側(cè)。由于流速較大同時(shí)縱向彌散度較小，所形成的污染羽形狀較為扁平。在離污染源最遠(yuǎn)處的濃度最小，在污染羽中心且靠近源強(qiáng)處的濃度最大，呈現(xiàn)出中心高四周低的規(guī)律。

在服務(wù)期滿后，廢石場(chǎng)下滲淋濾液量減少為零。污染物繼續(xù)沿著地下水主要流動(dòng)方向遷移，在服務(wù)期滿前后污染物遷移到達(dá)居民點(diǎn)，并且在河流1的下游的污染帶寬度逐漸擴(kuò)大。污染羽仍然呈現(xiàn)出中心濃度高而邊緣濃度低的特征。污染物匯入河流1后，隨著河流1從南向北向下游遷移，使得研究區(qū)域內(nèi)的污染物不斷減小，所形成的污染羽長(zhǎng)度逐漸減小，中心也逐步向河流2的下游移動(dòng)。在30年后，大部分污染物已經(jīng)流出研究區(qū)域，居民點(diǎn)也不再處于污染羽的范圍內(nèi)。污染物在研究區(qū)域內(nèi)總體快速遷移而不滯留。

3.4.3 對(duì)居民點(diǎn)的影響

圖3 污染物在各時(shí)期的濃度分布圖Fig.3 Concentration of Pollution in different times

居民點(diǎn)內(nèi)與整個(gè)研究區(qū)域內(nèi)污染物的最大濃度隨天數(shù)變化規(guī)律曲線如圖4所示。總體來看，研究區(qū)和居民點(diǎn)的最大污染物濃度變化規(guī)律較為一致，最大濃度隨天數(shù)變化曲線形狀呈現(xiàn)為寬峰。最初研究區(qū)域內(nèi)的最大濃度為廢石場(chǎng)下滲處的源強(qiáng)濃度，隨著時(shí)間增加最大濃度也增加，在5500天之后區(qū)域內(nèi)的最大濃度達(dá)到最大的0.2mg/L，并保持穩(wěn)定。在8700天左右，區(qū)域內(nèi)的最大濃度開始下降，到10000天左右降低到0.05mg/L度，隨著時(shí)間增加最大濃度也增加，在5500天之后區(qū)域內(nèi)的最大濃度達(dá)到最大的0.2mg/L，并保持穩(wěn)定。在8700天左右，區(qū)域內(nèi)的最大濃度開始下降，到10000天左右降低到0.05mg/L以下。由于污染物到達(dá)居民點(diǎn)需要一定的時(shí)間，所以0-3900天居民點(diǎn)內(nèi)未檢測(cè)到污染物。3900天后居民點(diǎn)內(nèi)最大濃度開始隨時(shí)間上升，在6000天左右達(dá)到最大值。居民點(diǎn)內(nèi)保持該最大濃度到8700天后開始下降，到11000天左右，居民點(diǎn)內(nèi)無法檢出污染物。

圖4 研究區(qū)與居民區(qū)的最大鉛離子濃度變化圖Fig.4 Maximum concentration of pollution in research area and residence area

居民點(diǎn)鉛離子的最大檢出濃度為0.0386mg/L。根據(jù)GB/T14848-93Ⅲ類水體標(biāo)準(zhǔn)，Pb≤0.05mg/L。居民點(diǎn)內(nèi)的鉛離子濃度未超標(biāo)。以0.003mg/L的檢出限值為標(biāo)準(zhǔn)，污染物到達(dá)居民點(diǎn)需要3910天左右。而在研究區(qū)域內(nèi)，大部分時(shí)間鉛離子濃度超過了0.05mg/L，且最大濃度達(dá)到了0.2mg/L。

4 結(jié)論

本文根據(jù)西藏某鉛鋅礦的水文地質(zhì)條件，建立了地下水滲流與廢石淋濾液在地下水系統(tǒng)中運(yùn)移的數(shù)學(xué)模型，數(shù)值結(jié)果合理可靠，對(duì)于定量化研究地下水污染防治提供理論依據(jù)。

（1）研究區(qū)域內(nèi)的水位與高程基本吻合，地下水流場(chǎng)受地形控制，自西北向東南的流動(dòng)為地下水的主要滲流方向，地下水沿河流1向下游流動(dòng)，與河流2匯集后向河流2下游流動(dòng)，有利于污染物的擴(kuò)散而不累積。

（2）污染物沿著地下水主要流動(dòng)方向由西北向東南方向遷移，同時(shí)具有向內(nèi)部河流形成的河谷匯集的趨勢(shì)，當(dāng)污染物到達(dá)河流2后將沿河流2下游方向遷移離開研究區(qū)域。整體來看，運(yùn)移過程中所形成的污染羽形狀較為扁平。在離污染源最遠(yuǎn)處的濃度最小，在污染羽中心且靠近源強(qiáng)處的濃度最大，呈現(xiàn)出中心高四周低的規(guī)律。

（3）污染物到達(dá)居民點(diǎn)需要約3910天，污染物的最大檢出濃度為0.0386mg/L，未超過地下水Ⅲ類水質(zhì)量標(biāo)準(zhǔn)限值，因此模擬結(jié)果表明淋濾液對(duì)居民點(diǎn)沒有造成風(fēng)險(xiǎn)和危害。但在5000-8000天的時(shí)段內(nèi)，需要加強(qiáng)對(duì)敏感點(diǎn)的監(jiān)測(cè)，提前預(yù)防可能出現(xiàn)的風(fēng)險(xiǎn)并進(jìn)行實(shí)時(shí)相應(yīng)防治措施。

[1]焦甜,王軍霞,唐仲華,等.地下水溶質(zhì)運(yùn)移方程有限差分格式的實(shí)證研究[J].安全與環(huán)境工程,2016,23(3):17-23.

[2]孟展.西藏斯弄多鉛鋅銀礦床礦物學(xué)特征研究及地質(zhì)意義[D].成都理工大學(xué),2016:23-39.

[3]陳華君,胡若鵬,鄧競(jìng)超.淋濾條件下鉛鋅尾礦中重金屬對(duì)環(huán)境影響研究[J].云南冶金,2016,45(1):56-63.

[4]盧曉華.基于數(shù)值模擬的企業(yè)地下水重金屬污染的環(huán)境影響預(yù)測(cè)評(píng)價(jià)[J].安全與環(huán)境工程,2014,21(1):93-97.

[5]馮書順,王國(guó)瑞,馬自強(qiáng),等.基于Visual Modflow的礦井涌水量預(yù)測(cè)模擬研究[J].煤炭技術(shù),2016,35(2):239-242.

[6]鄧鼎興.某礦山地下水污染分區(qū)治理效果的數(shù)值模擬[J].安全與環(huán)境工程,2013,20(4):27-31.

[7]張曉葉.基于可滲透反應(yīng)格柵的地下水中2，4-二氯酚溶質(zhì)運(yùn)移模擬研究[D].北京工業(yè)大學(xué),2015,1-2.

[8]焦友軍,施小清,吳吉春,等.基于PHT3D的地下水中六價(jià)鈾吸附反應(yīng)運(yùn)移數(shù)值模擬[J].地下水,2015(2):8-10.

[9]Salmon S U,Prommer H,Park J,et al.A general reactive transportmodelingframeworkforsimulatingandinterpreting groundwater14Cageandδ13C[J].WaterResourcesResearch,2015,51(1):359-376.

[10]孫驍磊.銀北井渠結(jié)合灌區(qū)地表水地下水耦合模擬及優(yōu)化配置研究[D].寧夏大學(xué),2016,2-3.

[11]Ahad J M,Ganeshram C L,Bryant L M.,et al.Sources of n-alkanes in an urbanized estuary:Insights from molecular distributions and compound-specific stabl e a nd radiocarbonisotopes[J].MarineChemistry,2011,126(1),239-249.

[12]張人權(quán),等.水文地質(zhì)學(xué)基礎(chǔ)(第六版)[M].北京:地質(zhì)出版社,2011:1-184.

[13]薛禹群,吳吉春.地下水動(dòng)力學(xué)[M].北京:中國(guó)水利水電出版社,2009:1.

[14]仵彥卿.多孔介質(zhì)滲流與污染物遷徙數(shù)學(xué)模型[M].北京:科學(xué)出版社,2012.

[15]薛禹群,謝春紅.地下水?dāng)?shù)值模擬[M].北京:科學(xué)出版社,2007.

[16]GelharLW,Welty C,Rehfeldt K R.Reply to Comments on“A critical review of data on field-scale dispersion in aquifers”by S.P.Neuman[J].Water Resources Research,1993,29(6):1867-1869.