重金屬在土壤-地下水中的遷移模型研究
——以山西、陜西黃土高原塬區(qū)為例

張?zhí)飙^，王海芳

(中北大學(xué) 環(huán)境與安全工程學(xué)院，山西 太原 030051)

0 引 言

根據(jù)《2019中國生態(tài)環(huán)境狀況公報》顯示，全國2 830處淺層地下水水質(zhì)監(jiān)測井中，Ⅰ～Ⅲ類水質(zhì)監(jiān)測井占23.7%，Ⅳ類占30.0%，Ⅴ類占46.2%，可見我國地下水質(zhì)狀況較差，地下水污染問題已成為人們關(guān)注的焦點問題[1].

重金屬在降雨等因素作用下會遷移至地下水，地下水作為人類的主要水源之一，在飲用含有重金屬的地下水后會直接影響人體健康，嚴重時會造成人體患病[2].地下水環(huán)境一旦受到污染，很難使其恢復(fù)至原來的狀態(tài)[3].所以，研究重金屬在土壤非飽和帶的遷移轉(zhuǎn)化，預(yù)測污染物進入地下水的濃度具有重要的理論與現(xiàn)實意義[4].現(xiàn)階段研究人員解決環(huán)境污染問題的一個重要技術(shù)方法是將所研究的環(huán)境系統(tǒng)行為抽象為一個數(shù)學(xué)模型，這是進行定量研究工作的基礎(chǔ)[5].

黃土高原塬區(qū)土壤表面重金屬向地下水的遷移中，降雨的入滲水分對土壤重金屬淋失具有顯著影響[6].國內(nèi)外已經(jīng)研究出很多降雨入滲模型，如Green-Ampt、Smith-Parlange、Horton、Kostiakov模型等[7].Kostiakov模型計算簡單，能夠有效描述短期范圍的降雨入滲過程.方正三[8]對Kostiakov公式和Horton公式進行了修正，用于研究在暴雨條件下黃土高原地區(qū)土壤滲透與時間的關(guān)系.研究土壤包氣帶污染問題的主要模型包括經(jīng)驗?zāi)Ｐ蚚9-10]、對流-彌散模型[11-14]和隨機模型[15-16].其中，對流彌散方程是最常用、最基本的描述溶質(zhì)運移的數(shù)學(xué)模型[17].Sun等[18]對 4種土壤進行了小土柱滲透試驗和溶析儀實驗，以對流-擴散方程為基礎(chǔ)建立了銻在土壤中的遷移模型，對流-彌散模型可以解釋銻在土壤中的遷移過程.對污染物遷移至地下水的研究方面，鐘茂生等[19]采用三相平衡模型和SESOIL模型耦合地下水稀釋模型對北京市不同水文地質(zhì)條件下的 69種有機污染物和農(nóng)藥推導(dǎo)了基于保護地下水的土壤通用篩選值.基于此，本文將入滲模型、對流-彌散模型和地下水稀釋模型進行耦合，描述污染物基于降雨入滲條件向地下水的遷移過程.最后利用MATLAB編程對該模型進行了驗證.

1 重金屬污染物在土壤中運移的影響因素

重金屬污染物易溶于水，土壤中的微生物很難將其降解，重金屬污染物一旦進入土壤就會發(fā)生持久性的污染.重金屬在土壤中的遷移會受到自身理化性質(zhì)，土壤對流、擴散和彌散作用，土壤質(zhì)地，土壤水分等因素的影響.土壤質(zhì)地的粘土含量越高，土壤中存在的自由水分含量和土壤空隙都會相應(yīng)降低，重金屬在土壤中運移就會受到限制[20]；而土壤含水量越高，對流和擴散作用就越大.對于黃土高原地區(qū)，降雨量少，土壤表面蒸發(fā)量大，重金屬進入土壤后的遷移主要受降雨入滲水分的影響[21].

2 塬區(qū)水文地質(zhì)條件

黃土高原塬區(qū)含水層巖性為離石組粉土質(zhì)黃土中夾帶有多層古土壤及鈣質(zhì)結(jié)核層，結(jié)構(gòu)疏松且空隙、裂隙發(fā)育，其下部三門組粉質(zhì)粘土，結(jié)構(gòu)相對較密實，空隙、裂隙不甚發(fā)育，構(gòu)成黃土塬區(qū)的隔水底板[22].

黃土高原塬區(qū)潛水主要分布在較大的黃土塬體中，在塬面寬度大于400 m時均存在潛水，深度一般在50 m～80 m.

黃土高原塬區(qū)潛水補給的最主要來源為大氣降水.黃土塬區(qū)寬闊有利于降雨下滲，黃土顆粒細密.野外滲水試驗顯示黃土的垂直滲透系數(shù)為2.1 m/d～7.8 m/d，透水性較好[22].

3 降雨入滲條件下污染物運移模型的構(gòu)建

3.1 概念模型

黃土高原地區(qū)降雨稀少，降雨多集中在每年的夏季，蒸發(fā)強烈，研究區(qū)土壤表面到潛水面的厚度為H，土壤中污染物向下遷移主要依靠降雨入滲的垂向遷移，存在的微弱側(cè)向流動可忽略，故可將非飽和帶概化為垂向一維流.假設(shè)土壤包氣帶是均質(zhì)各向同性的，污染物運移模型為對流-彌散方程，模型上邊界設(shè)定為定濃度邊界，下邊界設(shè)定為零濃度邊界.

污染物遷移模型見圖1.

圖1 研究區(qū)污染物遷移概念模型圖

3.2 降雨入滲模型

在降雨條件下的降雨入滲量采取Kostiakov 三參數(shù)入滲經(jīng)驗?zāi)Ｐ停撃Ｐ褪呛芏鄧业貐^(qū)使用的入滲模型，研究比較成熟且入滲參數(shù)的物理意義明確[23]，模型如下

I(t)=ξtα+Kst，

(1)

式中：I(t)為t時刻的累積入滲量，cm；ξ為入滲系數(shù)，表示入滲開始后第一個單位時間末扣除Ks后的累積入滲量，cm；α為入滲指數(shù)，表示土壤入滲速度的衰減速度；Ks為土壤相對穩(wěn)定入滲率，也稱為飽和導(dǎo)水率[24-25]，是在單位土壤勢梯度下飽和土壤的入滲速度或非飽和土壤入滲達到相對穩(wěn)定階段的入滲速度，cm/min.

降雨入滲條件下，由降雨強度及根系吸附作用得到的土壤包氣帶空隙水實際平均流速v為[26]

(2)

式中：ξ為降雨入滲系數(shù)；Jw為降雨入滲速度；CET為植物蒸騰系數(shù)；θ為包氣帶平均體積含水率；q為包氣帶滲流速度.

由于黃土高原地區(qū)干旱少雨，植被稀少，忽略植被蒸騰作用，省去由于植物蒸騰作用消耗的水量項，將式(1)代入式(2)得

(3)

由于黃土高原塬區(qū)土壤地下水補給主要為大氣降雨，污染物進入表土層后，隨著降雨入滲的土壤水分向地下遷移.由于各地區(qū)土壤質(zhì)地不同，降雨入滲量也不相同，從美國EPA獲取到各種土壤質(zhì)地的飽和導(dǎo)水率值，見表1，根據(jù)研究區(qū)土壤質(zhì)地條件即可知研究區(qū)土壤飽和導(dǎo)水率，利用式(2)可求得土壤孔隙水流速.

表1 土壤質(zhì)地對應(yīng)的飽和導(dǎo)水率值

3.3 污染物在土壤非飽和帶遷移模型

3.3.1 模型推導(dǎo)過程

污染物在土壤中的運移主要考慮對流、分子擴散和機械彌散3個物理過程以及土壤對污染物的吸附作用.

對流引起的溶質(zhì)通量與土壤水流通量和水中溶質(zhì)濃度有關(guān)，溶質(zhì)對流通量可表示為

Jc=qC,

(4)

式中：Jc為溶質(zhì)的對流通量，mol·m-2·s-1；q為水通量，m/s；C為單位體積土壤水中溶質(zhì)的量，mol/m3、kg/m3或g/L，即溶質(zhì)濃度.

由于土壤水流通量可表示為

q=vθ，

(5)

所以，式(4)可寫

Jc=vθC,

(6)

式中：v為平均空隙流速，m/s；θ為土壤體積含水量，m3/m3.

由于機械彌散和擴散在土壤中會引起溶質(zhì)濃度的混合和分散，且微觀流速不易測定，彌散和擴散結(jié)果不易區(qū)分，所以，將兩者聯(lián)合作為水動力彌散.用菲克第一定律表示的擴散作用為

(7)

式中：Js為溶質(zhì)的擴散通量，mol·m-2·s-1或kg·m-2·s-1；Ds為溶質(zhì)的有效擴散系數(shù)，m2/s；dC/dx為濃度梯度.土壤溶質(zhì)擴散可表示為

(8)

式中：Ds為擴散系數(shù).

機械彌散為

(9)

式中：Dh為機械彌散系數(shù)，m2/s.

將式(8)與式(9)聯(lián)立可得到水動力彌散通量為

(10)

式中：D為水動力彌散系數(shù).

污染物等溶質(zhì)在進入水土環(huán)境后，土壤固液相快速發(fā)生吸附、解吸作用，在很短時間或瞬間就能達到平衡狀態(tài).假設(shè)土壤中所吸附的溶質(zhì)量和土壤溶質(zhì)濃度是線性關(guān)系，采用Henry吸附模型(線性吸附)表示土壤與污染物之間的吸附

S=KdC,

(11)

式中：S為單位土壤吸附污染物的量；Kd為污染物分配系數(shù).

土壤溶質(zhì)運移主要是對流和水動力彌散(機械彌散和擴散)作用的結(jié)果，聯(lián)立式(4)和式(10)可得出溶質(zhì)通量為

(12)

假設(shè)土壤三維空間內(nèi)一單元六面體ABCD-EFGH(見圖2)為一單位容積土壤體，其邊長為Δx、Δy、Δz.在Δt時間內(nèi)，污染物進出單元體的變化符合質(zhì)量守恒原理，無源匯項存在，即進出該單元體的污染物的量的差等于Δt時段內(nèi)該單元體內(nèi)污染物質(zhì)量的變化.

圖2 污染物進入六面體單元示意圖

設(shè)進入六面體上界面ABCD面溶質(zhì)通量為Jx，則Δt時間內(nèi)由上界面進入的污染物的量為

mx=JxΔyΔzΔt，

(13)

流出六面體下界面EFGH面的溶質(zhì)通量為

(14)

Δt時間內(nèi)流出下界面EFGH的溶質(zhì)量為

(15)

在Δt時間內(nèi)沿x軸方向的污染物流入與流出單元體的污染物質(zhì)量差值為

(16)

同理，在Δt時間內(nèi)，沿y軸和z軸方向的污染物流入與流出質(zhì)量之差為

(17)

(18)

所以，在x，y，z3個方向上溶質(zhì)流入量與流出量的總差量為

(19)

根據(jù)質(zhì)量守恒定律，Δt時間內(nèi)該單元體中溶質(zhì)量的變化為流入與流出單元體中的溶質(zhì)質(zhì)量的差值，即為

(20)

將式(20)兩邊除以ΔxΔyΔzΔt得

(21)

利用愛因斯坦求和約定表示為

(22)

在一維條件下

(23)

將式(12)代入式(23)得一維條件下的對流彌散方程為

(24)

一維條件下考慮污染物在土壤的吸附交換項，則式(24)變?yōu)?/p>

(25)

將式(5)和式(11)代入式(25)，考慮對流、水動力彌散和吸附作用，污染物在黃土高原塬區(qū)土壤非飽和帶中遷移的一維水動力-彌散方程為

(26)

式中：ρ為土壤質(zhì)量密度；θ為土壤的體積含水率；C為土壤中污染物濃度；Dx為垂直方向上水動力彌散系數(shù).

將阻滯因子引入方程，降雨因素以源匯項考慮，則式(26)變?yōu)?/p>

(27)

式中：Rd為阻滯因子；Q為降雨導(dǎo)致的源匯項.

3.3.2 模型邊界條件及解析解

研究區(qū)污染源主要為連續(xù)或短期釋放源，并且考慮降雨和污染物在土壤中的吸附，且輸入濃度為一個常量，根據(jù)水文地質(zhì)概念模型可將邊界條件定為如下表示：

初始條件

C(x,t)|t=0=Ci,-∞≤x≤+∞；

(28)

邊界條件

(29)

(30)

根據(jù)初始條件和邊界條件，式(27)的解析解為

C(x,t)=

(31)

(32)

B(x,t)=

(33)

3.4 污染物進入飽和地下潛水模型

污染物隨土壤水分進入潛水時，與地下潛水混合后被稀釋，污染物濃度會隨之降低.污染物進入地下水與地下水混合稀釋后的濃度預(yù)測模型主要采用箱式模型來完成.

Cgw=Cp/DAF,

(34)

式中：Cgw為混合稀釋后污染物的質(zhì)量濃度，mg/L；Cp為污染物接觸潛水面時的濃度，mg/L，通過對流-彌散方程進行確定；DAF為地下水稀釋因子.

地下水稀釋因子DAF也可以根據(jù)當?shù)氐牡叵滤牡刭|(zhì)條件計算求得[27]，公式如下

(35)

D=(0.011 2×L2)1/2+

(36)

式中：K為含水層導(dǎo)水率，m/a；i為水力梯度，m/m；D為混合區(qū)深度(污染物與地下水混合的深度)，m；I為地下水滲透速率；L為污染源到監(jiān)測點的水平距離；Hgw為含水層厚度.

4 模型的優(yōu)勢與對比分析

目前，對于重金屬在非飽和-飽和土壤中遷移的整體性模型研究還不多.姜利國等[28]考慮了重金屬在地下水中的對流、彌散、吸附以及微生物降解作用，在非飽和-動理論及溶質(zhì)運移理論基礎(chǔ)上建立了非飽和-飽和區(qū)域內(nèi)水流方程與重金屬污染物運移方程耦合的數(shù)學(xué)模型.李錫夔[29]提出了一個考慮了污染物運移的對流、機械逸散、分子彌散、吸附、蛻變、不動水效應(yīng)的模擬飽和-非飽和土壤中污染物運移過程的數(shù)值模型.利用算例驗證了在一二維條件下污染物的運移狀況.劉培斌等[30]建立并驗證了在排水條件下田間一維飽和-非飽和土壤中氮素運移與轉(zhuǎn)化的耦合模型，模型中考慮了有機質(zhì)的礦化、氮素的吸附、硝化、反硝化、氨氣揮發(fā)及作物根系吸氮等氮素轉(zhuǎn)化作用過程，同時也考慮了土壤溫度和濕度對氮素轉(zhuǎn)化的影響.Grift B V D等[31]研究出一種模型方法來評估污染物對地下水和地表水負荷的影響，耦合非飽和帶淋溶模型和三維地下水運移模型，采用質(zhì)量平衡方法對3個冶煉廠附近的3個不同集水區(qū)進行了污染物預(yù)測模擬.

依據(jù)以上論述，目前國內(nèi)外在該方面的研究主要是實驗室結(jié)合土柱實驗或?qū)S區(qū)小范圍內(nèi)重金屬在土壤-地下水中的遷移模型研究，無法推廣使用.本文所建立的模型針對不同的暴露情景，選取黃土高原塬區(qū)為研究對象，將降雨入滲模型、對流彌散方程與地下水稀釋模型進行耦合，該耦合模型充分考慮降雨、水文、地質(zhì)及污染物特征等因素.在污染物進入土壤后用一維對流-彌散方程進行描述，給定具體的邊界條件即可計算得到污染物遷移至潛水面的污染物濃度.污染物進入地下水后用地下水稀釋模型描述，充分考慮研究區(qū)的水文地質(zhì)條件，進而計算污染物與地下水混合稀釋后的濃度.將土壤飽和導(dǎo)水率與土壤質(zhì)地結(jié)合在一起確定飽和導(dǎo)水率.前人研究的模型的計算方法都比較復(fù)雜，考慮的因素太多不易操作，很少考慮水文地質(zhì)及降雨等因素.

5 數(shù)學(xué)模型驗證

假設(shè)土壤中污染物的環(huán)境背景值為零，地表除污染源以外無任何流體進入模擬區(qū)域.參考文獻[32-34]，設(shè)置相應(yīng)參數(shù)為：土壤質(zhì)量密度ρ=1.6 g/cm3，土壤的體積含水率θ=0.4 cm3/cm3；垂直方向上水動力彌散系數(shù)Dx=400 cm2/d；土壤水平均孔隙流速v=20 cm/d；液相和吸附相間溶質(zhì)分配系數(shù)為1.17 cm3/g；模擬深度為潛水埋深x=50 m；模擬時間為365 d.

初始條件與邊界條件：

C(x,0)=0，x>0；

C(0,t)=8.96 mg/L, 0

C(0,t)=C(x,10),t≥10；

C(∞,t)=0,t>0.

利用MATLAB編程模擬污染源存在10 d時的重金屬污染物運移過程，分別獲得10 d，20 d，40 d時土壤中重金屬污染物質(zhì)量濃度分布及50 m深度處運移時間與重金屬濃度之間的關(guān)系圖，如圖3～圖6 所示.

圖3 10 d時重金屬與包氣帶深度的關(guān)系曲線

圖4 20 d時重金屬與包氣帶深度的關(guān)系曲線

圖5 40 d時重金屬與包氣帶深度的關(guān)系曲線

圖6 50 m深度處重金屬與運移時間的的關(guān)系曲線

從模擬結(jié)果看出，當污染源是瞬時污染源時，在釋放重金屬污染物階段，重金屬濃度分布隨著距離的增加而減小，直至在空間分布上趨于穩(wěn)定.這是因在靠近污染源處，污染源濃度較大，經(jīng)過一段時間的彌散，濃度會逐漸減小，因彌散速度較小，在短時間內(nèi)，重金屬只能彌散很短距離，較深土壤中污染物濃度較少.當重金屬污染源不存在后，重金屬在土壤各深度處的分布先迅速增加然后減小，這是因為在沒有新的污染源進入時，重金屬濃度逐漸減小，重金屬遷移過程中濃度梯度會越來越小，使某一點的重金屬濃度有一個積累的過程，然后又進行釋放，這與宋新山等[4]的研究基本吻合.

6 結(jié) 論

研究重金屬在飽和-非飽和土壤中遷移的整體性模型不多，本研究將降雨入滲與重金屬所在飽和非飽和土壤作為一個整體考慮，通過將入滲模型、對流-彌散模型和地下水稀釋模型進行耦合，達到計算基于保護地下水的土壤風(fēng)險值的目的.與前人研究的飽和-非飽和帶遷移模型進行對比分析，本模型考慮了降雨以及當?shù)厮牡刭|(zhì)條件，模型操作簡單，普適性強.利用MATLAB軟件編程對模型在給定初始條件及邊界條件下進行了驗證，與前人研究結(jié)果基本吻合，能夠預(yù)測污染物在包氣帶的遷移規(guī)律以及重金屬到達潛水面的濃度.本次分析只以我國黃土高原塬區(qū)為研究對象，后期研究應(yīng)該對我國各個地區(qū)地下水進行污染調(diào)查，選擇典型樣地，并結(jié)合當?shù)厮牡刭|(zhì)條件、氣候條件、地下水類型和埋深條件建立場地暴露概念模型，為后期重金屬的地下水質(zhì)量評估提供理論支撐.