劉歡， 孫洪廣

(河海大學(xué)力學(xué)與材料學(xué)院， 水文水資源與水利工程科學(xué)國家重點實驗室， 南京 210098)

20世紀(jì)及21世紀(jì)初，由于生產(chǎn)條件落后以及環(huán)保意識薄弱，許多工礦企業(yè)(特別是化工企業(yè))在生產(chǎn)過程中對附近土壤和地下水造成了嚴(yán)重的重金屬污染[1-2]。近年來，中國對生態(tài)環(huán)境問題越來越重視，老舊工業(yè)企業(yè)和廢棄場地的重金屬污染修復(fù)工程越來越多。作為土壤和地下水修復(fù)的基礎(chǔ)，數(shù)值模擬有助于掌握地下水溶質(zhì)運移規(guī)律，了解污染物遷移的趨勢，確定污染范圍及濃度分布。因此，數(shù)值模擬方法已被廣泛應(yīng)用于污染場地的定量評估與準(zhǔn)確預(yù)測，對水土污染管控和高效修復(fù)提供科學(xué)指導(dǎo)。

目前，關(guān)于土壤中重金屬遷移規(guī)律的研究已有大量報道。已有研究主要集中在兩方面：一是通過室內(nèi)土柱淋溶與砂箱實驗，結(jié)合理論分析考察重金屬在土壤中的遷移規(guī)律；二是通過數(shù)值方法對大型場地重金屬遷移進(jìn)行模擬與預(yù)測。江建明等[3]通過砂箱物理實驗與數(shù)值模擬發(fā)現(xiàn)裂隙對重金屬污染物遷移有導(dǎo)向作用。王玉雪等[4]通過室內(nèi)土柱實驗得到初始pH越低或初始濃度越高，重金屬Cr(Ⅵ)穿透時間與達(dá)到吸附飽和時間越早的結(jié)論。林青等[5]通過室內(nèi)置換實驗對銅、鎘、鋅、鉛在石英砂中的遷移規(guī)律進(jìn)行了研究，并通過數(shù)值模擬得出靜態(tài)吸附實驗得到的吸附參數(shù)不能直接應(yīng)用于重金屬運移模擬。Opekunova等[6]以某銅礦為研究對象，分析了預(yù)測銅、鋅、鎳及鎘等重金屬元素在土壤中的遷移。辛圓心等[7]針對獅子山尾礦，用COMSOL軟件考察了銅和鋅離子在土壤中的遷移規(guī)律。Lafratta等[8]通過對澳大利亞某鉛鋅冶煉廠的土壤進(jìn)行多參數(shù)研究，重現(xiàn)重金屬遷移過程，證實該冶煉廠仍是土壤重金屬污染的主要來源。Dong等[9]通過數(shù)值模擬重金屬在垃圾填埋場的遷移行為，得出在土壤與地下含水層最敏感的參數(shù)是阻滯因子與滲透系數(shù)。Zhang等[10]分析了鐵尾礦中重金屬的遷移特征及分布規(guī)律。孫洪廣等[11]研究表明，分?jǐn)?shù)階導(dǎo)數(shù)模型能很好的描述溶質(zhì)在土壤環(huán)境中的反常擴(kuò)散行為。

基于此，首先以南京某廢棄化工廠建立三維地下水重金屬溶質(zhì)運移模型，結(jié)合有限元方法，對該土壤中重金屬砷[As(V)]的濃度演化進(jìn)行模擬及預(yù)測；其次，基于該模型考察不同滲透系數(shù)、彌散度、分配系數(shù)以及孔隙率的條件下，重金屬As(V)在土壤中的遷移規(guī)律；最后，引入敏感性指數(shù)，比較各參數(shù)變化對重金屬As(V)溶質(zhì)遷移的影響。為類似地區(qū)的水土修復(fù)與地下水環(huán)境影響評價工作提供科學(xué)依據(jù)。

1 數(shù)學(xué)模型的建立

1.1 潛水層水流運動數(shù)學(xué)模型

(1)

式(1)中：h為地下水水頭；kxx、kyy、kzz分別為x、y、z方向的滲透系數(shù)；w為水流源匯項，包括降水入滲、河流入滲補(bǔ)給、井的抽水量等；Ss為給水系數(shù)；H0為含水層初始水位；H1為定水頭邊界條件；q為第二類邊界的單寬流量，流入為正，流出為負(fù)，隔水邊界為0；t為時間；q1為定流量邊界條件；Γ1為模型定水頭邊界；Γ2為垂向定流量邊界。

1.2 重金屬溶質(zhì)運移數(shù)學(xué)模型

i,j=x,y,z;xx=x;xy=y;xz=z

(2)

所研究的重金屬砷在土壤中的遷移，只考慮對流、彌散以及吸附作用，引入阻滯因子R，其表達(dá)式為

(3)

(4)

式中：C為土壤溶質(zhì)濃度；θ為土體孔隙率；ρ為土壤干體積質(zhì)量即容重；Q為溶質(zhì)源匯項；Dij為彌散系數(shù)；qi為達(dá)西流速；t為時間；Rn為土壤溶質(zhì)各個轉(zhuǎn)化過程包括生物吸收、化學(xué)反應(yīng)、降解、沉淀等過程；Kd為重金屬溶質(zhì)吸附相與溶解相的平衡分配系數(shù)。

2 幾何建模與邊界條件

主要研究中國南京某污染地塊的重金屬As(V)在地下水中的運移規(guī)律以及隨時間變化的濃度分布狀態(tài)。通過對該地塊的實地調(diào)查，建立重金屬在土壤中運移的三維模型。該模型可概化為250 m×160 m的矩形地塊(圖1)，潛水含水層約為10 m。東西南北四面皆臨近河流，地下水流大致為自東向西流動?？蓪|面和西面設(shè)定為恒定水頭，南北方向為沿程均勻下降的定水頭邊界。垂向邊界系統(tǒng)，上邊界為潛水含水層自由水面，系統(tǒng)與外部垂向水量交換都是通過該邊界完成的，如大氣降水補(bǔ)給，研究區(qū)域下部處理為零通量邊界。污染源為該廢棄產(chǎn)地的焚燒爐范圍，通過大氣降水淋溶進(jìn)入地下水。

為了考察各參數(shù)對重金屬As(V)在土壤中遷移的影響，以該模擬案例為基準(zhǔn)，將孔隙率、滲透系數(shù)、彌散度、吸附分配系數(shù)分別上下波動25%、50%作為4個水平值，連續(xù)通入10年污染物，研究分析各種參數(shù)變化對溶質(zhì)運移的影響。并采用因子變換法[12]，以敏感性指數(shù)為量化指標(biāo)對各參數(shù)進(jìn)行分析。具體參數(shù)如表1所示。

圖1 三維幾何建模圖Fig.1 3D geometric modeling drawing

表1 各參數(shù)水平值Table 1 Parameter level values

采用以下敏感性指數(shù)計算公式對模擬結(jié)果進(jìn)行敏感性評價，L值越大，表明參數(shù)改變對污染物遷移距離影響越大，即該參數(shù)的敏感性越高。計算公式為

(5)

式(5)中：L為敏感性指數(shù)，無量綱；n為參數(shù)水平值個數(shù)，取值為4；cj為對照組重金屬As(V)縱向遷移距離；c′j為其他水平值重金屬遷移距離。

3 計算結(jié)果及分析

3.1 重金屬砷溶質(zhì)遷移模擬結(jié)果

該化工廠生產(chǎn)含有重金屬As(V)產(chǎn)品滿7年后停止運營，隨后半年完成搬遷。

圖2 重金屬砷濃度分布隨時間變化Fig.2 Variation diagram of arsenic concentration distribution with time

圖2為土壤中重金屬砷的分布隨時間的變化，可以看出，前7年污染源重金屬As(V)持續(xù)泄露。污染暈中心以及污染羽中As(V)的質(zhì)量濃度逐年升高，預(yù)測時段分別為2、5、7年時，污染暈中心濃度達(dá)到2.25、5.55、7.12 mg/L；按照《地下水質(zhì)量標(biāo)準(zhǔn)》(GB/T 14848—2017)，重金屬砷標(biāo)準(zhǔn)值為0.05 mg/L，濃度分布如圖2所示，污染面積隨著時間推移不斷增加，超標(biāo)范圍也不斷增加。

7年后，由于該工廠搬遷，污染面源處不再流入污染物。在地下水對流作用下，污染范圍整體向下游遷移，在第20年、40年時，污染暈中心位置往下游分別遷移了2.7 m與6.5 m；彌散作用使得污染暈中心濃度向四周擴(kuò)散，污染暈中心濃度逐漸減少，污染范圍增加明顯，第40年時污染暈中心濃度已下降到0.33 mg/L。

同時可以注意到重金屬As(V)污染物在地下水中的遷移并不是向四周均勻擴(kuò)散的，污染羽縱向遷移距離更遠(yuǎn)，橫向遷移距離較短，在第7年時，污染羽縱向遷移距離達(dá)38.5 m，橫向遷移距離僅23.7 m。在垂向濃度主要集中在土壤表層，污染物濃度由表層向下減少，第20年超標(biāo)深度才達(dá)1.5 m，這是因為縱向是水力梯度的主要方向，該方向?qū)α饕约皬浬⑾禂?shù)值都比其他方向大。

3.2 模擬結(jié)果對不同參數(shù)的響應(yīng)

為研究分析各水文地質(zhì)參數(shù)對重金屬As(V)溶質(zhì)遷移的規(guī)律，在上述模型的基礎(chǔ)上，只考慮連續(xù)通入10年污染物的情況。并且以污染源中心及縱向濃度變化為主要研究對象，分析不同參數(shù)變化，重金屬As(V)溶質(zhì)遷移對其的響應(yīng)規(guī)律。

3.2.1 滲透系數(shù)對重金屬砷溶質(zhì)遷移的影響

當(dāng)滲透系數(shù)增加時，重金屬As(V)溶質(zhì)在縱向達(dá)到更大的遷移距離；另外，滲透系數(shù)也影響著溶質(zhì)在縱向分散上的傳輸，污染暈中心濃度隨著水力滲透系數(shù)的增加而減小。圖3為其他參數(shù)不變，滲透系數(shù)K分別取2.4×10-5、3.6×10-5、4.8×10-5、6.0×10-5、7.2×10-5m/s時，重金屬砷通入第10年的縱向濃度變化。滲透系數(shù)取值為2.4×10-5m/s時，污染暈中心濃度達(dá)到11.2 mg/L，相對參照組增幅達(dá)到29.2%，向下游遷移距離減少了6 m；滲透系數(shù)取值為7.2×10-5m/s時,污染暈中心濃度降低到7.1 mg/L，重金屬As(V)向下游遷移距離可達(dá)24.5 m，相比參照組污染源中心濃度降幅為19.1%，往下游遷移距離增加了4 m。這是因為在水頭差一定的情況下，一方面，滲透系數(shù)的增加直接導(dǎo)致地下水達(dá)西速度及滲流速度的增加，污染物隨地下水流運動，加快了溶質(zhì)向下游的遷移，從而污染暈中心濃度降低，向下游遷移距離增加；另一方面縱向彌散系數(shù)也相應(yīng)增大，隨著滲透系數(shù)增加，圖3中重金屬As(V)往水流上游遷移距離也呈增加趨勢。這也側(cè)面說明，污染物在吸附力度較強(qiáng)的土壤中，彌散作用對重金屬As(V)的遷移影響不可忽視。

圖3 不同滲透系數(shù)重金屬砷縱向濃度變化曲線Fig.3 Longitudinal concentration curves of heavy metal arsenic with different hydraulic conductivity

3.2.2 彌散度對重金屬砷溶質(zhì)遷移的影響

彌散度是用來描述含水層中污染物彌散作用的參數(shù)，彌散系數(shù)除以孔隙平均流速便是彌散度，其取值與土的孔隙結(jié)構(gòu)有關(guān)，且具有明顯的尺度效應(yīng)[13-14]。當(dāng)彌散度增加時，重金屬As(V)溶質(zhì)在土壤中遷移的變化趨勢與滲透系數(shù)增加一致，重金屬As(V)溶質(zhì)縱向遷移距離增加，污染源中心濃度減小。圖4為其他參數(shù)不變，αL分別取1.750、2.625、3.500、4.375、5.250 m時，重金屬溶質(zhì)通入第10年的縱向濃度變化?？v向彌散度取1.750 m時，污染源中心濃度為11.14 mg/L，與對照組彌散度取3.5 m相比，污染源中心濃度增加了28.4%，往下游遷移距離減少了4.5 m?？v向彌散度取5.250 m時，污染源中心濃度為7.13 mg/L，向下游遷移距離為23.5 m，與對照組相比，中心濃度減少17.8%，往下游遷移距離增加3.5 m。這是因為在地下水流速一定情況下，隨著彌散度增大，彌散系數(shù)也相應(yīng)增大，土壤中重金屬As(V)溶質(zhì)受到彌散作用越強(qiáng)，從而污染源溶質(zhì)向四周彌散越快，表現(xiàn)為污染源中心濃度越低，遷移距離越遠(yuǎn)，該數(shù)值模擬結(jié)果上游以及下游都是如此情況。

圖4 不同彌散度重金屬砷縱向濃度變化曲線Fig.4 Longitudinal concentration curve of heavy metal arsenic with different dispersion

3.2.3 分配系數(shù)對重金屬砷溶質(zhì)遷移的影響

圖5 不同分配系數(shù)重金屬砷縱向濃度變化曲線Fig.5 Longitudinal concentration curves of heavy metal arsenic with different partition coefficients

圖5為其他參數(shù)不變，分配系數(shù)分別取10、15、20、25、30 L/kg，重金屬溶質(zhì)通入第10年的縱向濃度變化。結(jié)果表明，重金屬As(V)在土壤中的遷移距離隨吸附分配系數(shù)的增加而減小，污染暈中心濃度與分配系數(shù)的大小也呈負(fù)相關(guān)變化。Kd為30 L/kg時，運移10年污染源中心濃度6.9 mg/L，向下游遷移17 m，相對參照組Kd=20 L/kg，污染源中心濃度降幅27.6%，向下游遷移距離減少了3 m；Kd為10 L/kg時，污染暈中心濃度為11.6 mg/L，向下游遷移28 m，與參照組相比，污染源中心濃度增幅34.4%，遷移距離增加了8 m。這是由于分配系數(shù)越大，土壤對污染物的吸附能力越強(qiáng)，地下水中越多重金屬As(V)溶質(zhì)被吸附到土壤中，導(dǎo)致地下水中As(V)濃度減少；土壤對重金屬As(V)污染物吸附增強(qiáng)，也表現(xiàn)為砷在土壤中受到的阻滯作用越強(qiáng)，溶質(zhì)在土壤中的遷移速度越慢，因而遷移距離越短[9]。

3.2.4 孔隙率對重金屬砷溶質(zhì)遷移的影響

孔隙率對溶質(zhì)遷移的影響在兩方面，一是孔隙率直接決定滲流速度，而滲流速度控制對流遷移；二是決定模型單元中儲存溶質(zhì)的孔隙體積大小。在給定滲透系數(shù)情況下，即達(dá)西流速不變，孔隙率對溶質(zhì)的運移作用主要表現(xiàn)在污染物濃度分散上面，孔隙率越大，孔隙實際平均流速越小，水動力彌散系數(shù)也越小，污染暈中心的濃度越高，污染物遷移距離越短[15]。在分配系數(shù)較小或者可以忽略的情況下的確如此，然而對于土壤吸附非常強(qiáng)的重金屬As(V)而言，污染物并非以水的滲流速度v做對流，其孔隙流速為v/R，對應(yīng)的彌散系數(shù)D也是如此，即D/R[16]。在0～1改變θ值，對污染物滲流速度以及彌散度幾乎沒有影響。圖6為其他參數(shù)不變，孔隙率分別為0.210、0.315、0.420、0.525、0.630時，重金屬砷通入第10年的縱向濃度變化。圖6中孔隙率θ在0.21～0.63變化，污染暈中心濃度相對于參照組變化幅度為-0.04～0.04 mg/L，遷移距離幾乎無差別。由此可見，在土壤對污染物吸附力度較強(qiáng)的土壤中，相對來說孔隙率的變化對污染物遷移的影響較弱。

圖6 不同孔隙率重金屬砷縱向濃度變化曲線Fig.6 Longitudinal concentration curve of heavy metal arsenic with different porosity

3.2.5 各參數(shù)對重金屬溶質(zhì)遷移影響的比較

各參數(shù)變化對污染物遷移距離的影響如表2所示，可以看出，在0～50%波動范圍，滲透系數(shù)、彌散度以及吸附分配系數(shù)對重金屬As(V)遷移距離的影響程度相當(dāng)，敏感指數(shù)分別為18.0%、16.3%、18.0%，孔隙率對重金屬As(V)的遷移影響極弱。這與不考慮吸附作用或者分配系數(shù)較小的溶質(zhì)運移規(guī)律存在較大偏差[17-18]，其原因是重金屬As(V)溶質(zhì)分配系數(shù)比較大，土壤對其的吸附能力比較強(qiáng)，導(dǎo)致在彌散度與孔隙率對重金屬As(V)溶質(zhì)遷移規(guī)律的影響上存在變化。一是重金屬溶質(zhì)遷移過程中，由于土壤的吸附作用被大量阻滯在土壤中，遷移距離受彌散作用影響較大，從而彌散度對重金屬的敏感性增強(qiáng)；二是重金屬溶質(zhì)在土壤中的實際對流遷移以及彌散系數(shù)都直接受阻滯因子大小的影響，而孔隙率的變化對其影響極小，所以在本文的模型中，孔隙率對溶質(zhì)的遷移作用影響可忽略。

表2 各參數(shù)不同幅度下重金屬砷遷移距離擬合偏差Table 2 Fitting deviation of migration distance of heavy metal arsenic under different ranges of each parameter

同時可以注意到，各參數(shù)的變化對重金屬As(V)在土壤中遷移距離的影響并不是呈線性變化，是因為其遷移距離受多種因素的影響，除了所探討的參數(shù)外，還有水壓差，降雨入滲系數(shù)等等。滲透系數(shù)及彌散度變化幅度在±25%以內(nèi)，可近似看作線性變化。

4 結(jié)論

基于南京某廢棄化工廠重金屬As(V)在土壤中的遷移，運用有限元方法進(jìn)行三維數(shù)值模擬，再在該模型基礎(chǔ)上，改變滲透系數(shù)、彌散度、分配系數(shù)及孔隙率的大小，對比分析重金屬溶質(zhì)在土壤中的遷移規(guī)律，并引入敏感性指數(shù)，比較各參數(shù)對重金屬溶質(zhì)遷移影響的強(qiáng)度。得出如下結(jié)論。

(1)該污染場地在連續(xù)通入重金屬As(V)污染物前7年，污染暈中心濃度持續(xù)上升到7.12 mg/L，污染范圍隨時間推移不斷擴(kuò)大，第7年縱向遷移距離38.5 m，遠(yuǎn)大于橫向遷移距離23.7 m。停止通入污染物后，在地下水對流作用下，污染范圍整體向下游移動，同時由于彌散作用，污染暈中心濃度持續(xù)下降，在第40年時污染暈中心往下游遷移了6.5 m，濃度下降到0.33 mg/L。

(2)重金屬砷As(V)在土壤中的遷移距離隨滲透系數(shù)、彌散度的增加而增加，隨分配系數(shù)增加而減小。并且由于土壤對重金屬As(V)的吸附能力強(qiáng)，在土壤阻滯作用下，重金屬As(V)在土壤中實際遷移速度較小，分配系數(shù)對重金屬的遷移影響較大，而孔隙率的變化卻影響極小，幾乎可忽略。

(3)通過對各參數(shù)的敏感性分析，得出在該模型中滲透系數(shù)、彌散度、分配系數(shù)對重金屬As(V)在土壤中的遷移分布影響程度相當(dāng)，敏感性指數(shù)分別為18.0%、16.3%、18.0%，孔隙度為0.006%。