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

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

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

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

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

1 研究區(qū)域概況

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

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

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

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

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

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

2.3 污染物遷移控制方程

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

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

水動力彌散系數(shù)包括機械彌散系數(shù)與分子擴散系數(shù)：

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

2.4 初始條件與邊界條件

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

地下水運移的定解條件：

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

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

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

邊界條件：

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

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

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

3.1 有限差分數(shù)值解法

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

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

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

3.2 定解條件的確定

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

3.3 參數(shù)選取

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

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

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

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

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

3.4 結(jié)果分析

3.4.1 地下水流場分析

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

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

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

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

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

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

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

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

3.4.3 對居民點的影響

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

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

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

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

4 結(jié)論

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

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

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

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

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

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

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

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

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

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

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

[8]焦友軍,施小清,吳吉春,等.基于PHT3D的地下水中六價鈾吸附反應(yīng)運移數(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]薛禹群,吳吉春.地下水動力學(xué)[M].北京:中國水利水電出版社,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.