重金屬砷在飽和土壤中遷移過程模擬

劉歡， 孫洪廣

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

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

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

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

1 數學模型的建立

1.1 潛水層水流運動數學模型

(1)

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

1.2 重金屬溶質運移數學模型

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

(2)

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

(3)

(4)

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

2 幾何建模與邊界條件

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

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

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

表1 各參數水平值Table 1 Parameter level values

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

(5)

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

3 計算結果及分析

3.1 重金屬砷溶質遷移模擬結果

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

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

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

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

同時可以注意到重金屬As(V)污染物在地下水中的遷移并不是向四周均勻擴散的，污染羽縱向遷移距離更遠，橫向遷移距離較短，在第7年時，污染羽縱向遷移距離達38.5 m，橫向遷移距離僅23.7 m。在垂向濃度主要集中在土壤表層，污染物濃度由表層向下減少，第20年超標深度才達1.5 m，這是因為縱向是水力梯度的主要方向，該方向對流以及彌散系數值都比其他方向大。

3.2 模擬結果對不同參數的響應

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

3.2.1 滲透系數對重金屬砷溶質遷移的影響

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

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

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

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

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

3.2.3 分配系數對重金屬砷溶質遷移的影響

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

圖5為其他參數不變，分配系數分別取10、15、20、25、30 L/kg，重金屬溶質通入第10年的縱向濃度變化。結果表明，重金屬As(V)在土壤中的遷移距離隨吸附分配系數的增加而減小，污染暈中心濃度與分配系數的大小也呈負相關變化。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。這是由于分配系數越大，土壤對污染物的吸附能力越強，地下水中越多重金屬As(V)溶質被吸附到土壤中，導致地下水中As(V)濃度減少；土壤對重金屬As(V)污染物吸附增強，也表現為砷在土壤中受到的阻滯作用越強，溶質在土壤中的遷移速度越慢，因而遷移距離越短[9]。

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

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

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

3.2.5 各參數對重金屬溶質遷移影響的比較

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

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

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

4 結論

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

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

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

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