張?zhí)飙^，王海芳

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

0 引 言

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

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

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

1 重金屬污染物在土壤中運(yùn)移的影響因素

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

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

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

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

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

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

3.1 概念模型

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

污染物遷移模型見(jiàn)圖1.

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

3.2 降雨入滲模型

在降雨條件下的降雨入滲量采取Kostiakov 三參數(shù)入滲經(jīng)驗(yàn)?zāi)Ｐ?，該模型是很多?guó)家地區(qū)使用的入滲模型，研究比較成熟且入滲參數(shù)的物理意義明確[23]，模型如下

I(t)=ξtα+Kst，

(1)

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

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

(2)

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

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

(3)

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

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

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

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

污染物在土壤中的運(yùn)移主要考慮對(duì)流、分子擴(kuò)散和機(jī)械彌散3個(gè)物理過(guò)程以及土壤對(duì)污染物的吸附作用.

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

Jc=qC,

(4)

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

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

q=vθ，

(5)

所以，式(4)可寫(xiě)

Jc=vθC,

(6)

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

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

(7)

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

(8)

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

機(jī)械彌散為

(9)

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

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

(10)

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

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

S=KdC,

(11)

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

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

(12)

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

圖2 污染物進(jìn)入六面體單元示意圖

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

mx=JxΔyΔzΔt，

(13)

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

(14)

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

(15)

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

(16)

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

(17)

(18)

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

(19)

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

(20)

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

(21)

利用愛(ài)因斯坦求和約定表示為

(22)

在一維條件下

(23)

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

(24)

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

(25)

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

(26)

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

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

(27)

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

3.3.2 模型邊界條件及解析解

研究區(qū)污染源主要為連續(xù)或短期釋放源，并且考慮降雨和污染物在土壤中的吸附，且輸入濃度為一個(gè)常量，根據(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 污染物進(jìn)入飽和地下潛水模型

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

Cgw=Cp/DAF,

(34)

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

地下水稀釋因子DAF也可以根據(jù)當(dāng)?shù)氐牡叵滤牡刭|(zhì)條件計(jì)算求得[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)測(cè)點(diǎn)的水平距離；Hgw為含水層厚度.

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

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

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

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

假設(shè)土壤中污染物的環(huán)境背景值為零，地表除污染源以外無(wú)任何流體進(jìn)入模擬區(qū)域.參考文獻(xiàn)[32-34]，設(shè)置相應(yīng)參數(shù)為：土壤質(zhì)量密度ρ=1.6 g/cm3，土壤的體積含水率θ=0.4 cm3/cm3；垂直方向上水動(dòng)力彌散系數(shù)Dx=400 cm2/d；土壤水平均孔隙流速v=20 cm/d；液相和吸附相間溶質(zhì)分配系數(shù)為1.17 cm3/g；模擬深度為潛水埋深x=50 m；模擬時(shí)間為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時(shí)的重金屬污染物運(yùn)移過(guò)程，分別獲得10 d，20 d，40 d時(shí)土壤中重金屬污染物質(zhì)量濃度分布及50 m深度處運(yùn)移時(shí)間與重金屬濃度之間的關(guān)系圖，如圖3～圖6 所示.

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

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

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

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

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

6 結(jié) 論

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