萬(wàn)朔陽(yáng)， 吳 勇， 唐學(xué)芳， 鄧東平， 蘭 真， 韓莉璧

(1.成都理工大學(xué)環(huán)境與土木工程學(xué)院，成都 610059；2.地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，成都 610059)

隨著科學(xué)技術(shù)水平的提高,農(nóng)業(yè)技術(shù)也不斷得到改進(jìn), 人類活動(dòng)與土壤之間的聯(lián)系日益增加，即來(lái)對(duì)土壤污染也逐漸加劇，現(xiàn)如今因土壤鹽堿化, 化肥、農(nóng)藥污染、污水灌溉引起的地下水和土壤污染等問(wèn)題日益受到關(guān)注。探究土壤溶質(zhì)運(yùn)移規(guī)律逐漸成為解決此類問(wèn)題的關(guān)鍵所在，由于影響土壤溶質(zhì)運(yùn)移因素眾多僅靠室內(nèi)試驗(yàn)難以掌控其實(shí)際運(yùn)移規(guī)律。故通過(guò)對(duì)受污染場(chǎng)地土壤建立數(shù)學(xué)模型來(lái)處理這些問(wèn)題已經(jīng)成為眾多學(xué)者的熱門研究對(duì)象[1-2]。早在20世紀(jì)90年代Simunek等[3]就已利用Hydrus-1D對(duì)非飽和土的水力和溶質(zhì)輸運(yùn)參數(shù)進(jìn)行估算，驗(yàn)算出各參數(shù)的控制方程，為分析土壤溶質(zhì)運(yùn)移特征提供了重要的理論基礎(chǔ)。郭瑞等[4]在此基礎(chǔ)上以水鹽運(yùn)移理論和環(huán)境條件為依據(jù)，建立了適宜于不同條件下的水鹽運(yùn)移模型。王銳等[5]在分析入滲壓力水頭對(duì)土壤水分入滲影響的基礎(chǔ)上，建立了不同入滲水頭條件下土壤水分運(yùn)動(dòng)數(shù)學(xué)模型，模擬值與實(shí)驗(yàn)值最大誤差不超過(guò)6%。劉目興等[6]對(duì)比研究了黏質(zhì)土壤在不同初始含水率條件下對(duì)初始入滲率的影響，結(jié)果表明土壤初始含水率與初始入滲率成反比，并且通過(guò)構(gòu)建Horton模型的對(duì)結(jié)論進(jìn)行擬合模擬，擬合效果較好。李茜等[7]提出構(gòu)建土壤水與地下水溶質(zhì)運(yùn)移聯(lián)合模型，將兩者相結(jié)合，更好的描述區(qū)域溶質(zhì)遷移的規(guī)律。由上述可知，對(duì)土壤建立溶質(zhì)運(yùn)移模型不僅能對(duì)運(yùn)移過(guò)程的分析提供理論依據(jù)，同時(shí)也能較好預(yù)測(cè)土壤污染物遷移規(guī)律，對(duì)土壤修復(fù)及防治有著重要意義。

選取四川西壩鎮(zhèn)某高產(chǎn)農(nóng)田為研究區(qū)域，對(duì)所取土樣進(jìn)行相關(guān)室內(nèi)實(shí)驗(yàn)參數(shù)分析的基礎(chǔ)上，運(yùn)用Hydrus-1D建立水流溶質(zhì)運(yùn)移模型，對(duì)恒定水流條件下土壤室內(nèi)淋濾實(shí)驗(yàn)結(jié)果進(jìn)行數(shù)據(jù)擬合；并且借助地下水建模系統(tǒng)(groundwater modeling systems,GMS)構(gòu)建空間分布模型，分析各項(xiàng)重金屬在研究區(qū)范圍內(nèi)的空間分布特征，對(duì)研究區(qū)內(nèi)土壤重金屬遷移規(guī)律及分布特征差異作出總體的評(píng)述。

1 研究區(qū)概況

研究區(qū)位于四川省樂山市五通橋區(qū)南西(方位253°)，直距約4.0 km處，面積2.4 km2(圖1)；此次共取72組土樣對(duì)研究區(qū)內(nèi)的高產(chǎn)農(nóng)田重金屬遷移及分布特征展開分析，各取樣點(diǎn)如圖2所示。

2 材料與方法

2.1 土樣實(shí)驗(yàn)方法

對(duì)研究區(qū)所取72組土樣采用高溫密閉消解法對(duì)土壤樣品進(jìn)行消解。采用ICP-MC對(duì)土壤消解產(chǎn)物進(jìn)行重金屬測(cè)試分析。重金屬檢測(cè)項(xiàng)目主要包括：Cd、Pb、Cr、Cu、Hg、As等，主要選取Cd、Cu、Zn、As、Cr等離子對(duì)區(qū)域內(nèi)土壤重金屬分布及遷移特征展開討論，實(shí)驗(yàn)數(shù)據(jù)統(tǒng)計(jì)見表1，均在成都理工大學(xué)地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室測(cè)得。

圖1 研究區(qū)位置及地貌

圖2 研究區(qū)采樣點(diǎn)分布

表1 土壤重金屬數(shù)據(jù)統(tǒng)計(jì)表

2.2 室內(nèi)淋濾實(shí)驗(yàn)方法

土柱實(shí)驗(yàn)所用土樣取至研究區(qū)某塊高產(chǎn)農(nóng)田深度為1 m的探坑內(nèi)，沿地表向下1 m，每隔10 cm取若干組土樣；實(shí)驗(yàn)土柱是用厚5 mm、長(zhǎng)20 cm，直徑為10 cm有機(jī)玻璃管制作而成，如圖3和圖4所示，柱頂連接供水管、柱底有滲漏孔，土柱共分為四層，按照探坑不同深度取樣依次分層填樣，第一層(0～5 cm)填充探坑0～20 cm所取土樣，依次填充，第四層(15～20 cm)為探坑80～100 cm所取土樣。實(shí)驗(yàn)采用馬氏瓶穩(wěn)定供水，淋濾液為超純水，淋濾時(shí)間為11天，每隔一天取一次滲漏水樣，記錄其中Cd、Cu、Zn、Cr離子含量變化值。

圖3 淋濾實(shí)驗(yàn)裝置

3 結(jié)果及分析

3.1 建立溶質(zhì)運(yùn)移模型

Hydrus-1D軟件[8]是一套用于模擬飽和-非飽和多孔介質(zhì)中水分運(yùn)移和溶質(zhì)運(yùn)移的數(shù)值模型，在對(duì)非飽和帶土壤水溶質(zhì)運(yùn)移模擬方面，具有操作簡(jiǎn)單，擬合效果較好等優(yōu)勢(shì)，因此廣泛被運(yùn)用于此類問(wèn)題的研究。通過(guò)確立溶質(zhì)基本控制方程 、邊界條件、土壤水力特征參數(shù)等過(guò)程，來(lái)構(gòu)建溶質(zhì)運(yùn)移模型。

3.1.1 模型基本控制方程

采用Hydrus-1D軟件，對(duì)研究區(qū)土壤垂向上建立一維水流溶質(zhì)運(yùn)移模型，此次水流模型概化為多層均質(zhì)各項(xiàng)同性的一維垂向穩(wěn)定流[9]，土壤水分運(yùn)動(dòng)主要受Richards方程控制，土壤溶質(zhì)運(yùn)移主要受對(duì)流-彌散方程控制[10]。

(1)

(2)

式中：θ為土壤體積含水率，cm3/cm3；h為土壤壓力水頭，cm；S為根系吸水項(xiàng)或其他源匯項(xiàng)，cm3/cm3/T(T為時(shí)間單位常量，可代表h、min、s)；K為非飽和土壤水力梯度傳導(dǎo)，cm/T；KAij為各項(xiàng)異性張量的分量，(-)；Xi為空間坐標(biāo)，cm；Z為垂項(xiàng)坐標(biāo)(向上為正)，cm；t為時(shí)間，T；c為溶質(zhì)濃度，g/cm3；D為水動(dòng)力彌散系數(shù)，cm2/h；q為達(dá)西流速，cm/h。

3.1.2 初始條件與邊界條件

本次模擬條件在不考慮表面蒸發(fā)及植物根系吸收作用下進(jìn)行，將初始條件設(shè)定為土柱中各液相濃度分布，對(duì)于水流模型上邊界為通量邊界，下邊界為自由滲漏邊界；對(duì)于溶質(zhì)運(yùn)移模型上邊界為大氣通量邊界，下邊界為零梯度邊界[11]。

(1)土壤水分運(yùn)移的初始條件與邊界條件。

初始條件：

h(x,0)=h0(x)， 0≤x≤L

(3)

上邊界條件：

(4)

下邊界條件：

h(x,t)=hb(t)

(5)

式中：h0(x)表示土壤剖面初始土壤水負(fù)壓,cm;hb(t)表示下邊界壓力水頭,cm；qs表示地表水入滲量。

(2)土壤溶質(zhì)運(yùn)移的初始條件與邊界條件。

初始條件：

c(x,0)=c0(x),t=0,0≤x≤L

(6)

上邊界條件：

(7)

下邊界條件：

(8)

式中：c0為初始剖面土壤液相濃度,g/cm3；cs為上邊界流量溶質(zhì)濃度，此次試驗(yàn)采用超純水故濃度cs=0。

3.1.3 土壤水力特征參數(shù)

土壤水力特征參數(shù)屬于土壤的基本物理性質(zhì),其各項(xiàng)參數(shù)屬性直接影響土壤水在包氣帶中的遷移轉(zhuǎn)化關(guān)系,對(duì)分析土壤溶質(zhì)運(yùn)移起重要作用[12]；通常采用van-Genuchten模型描述土壤水力特征參數(shù)[13]。為了將模型更好地與室內(nèi)實(shí)驗(yàn)結(jié)果相擬合，使模擬過(guò)程更具真實(shí)性，基于實(shí)驗(yàn)土壤不同粒徑分選占比量和土壤容重?cái)?shù)據(jù)，利用Hydrus-1D內(nèi)置的預(yù)測(cè)及反演土壤水力特征參數(shù)功能，并結(jié)合野外場(chǎng)地實(shí)驗(yàn)值對(duì)預(yù)測(cè)值進(jìn)行校正，最終推測(cè)出本研究區(qū)內(nèi)土壤的水力特征參數(shù)(表2)。模擬時(shí)間設(shè)為11 d，土柱分為2層，上層為砂質(zhì)壤土、下層為黃褐土粉砂壤土，虛擬觀測(cè)點(diǎn)位于土柱低端滲水口處，根據(jù)土壤自身所含重金屬含量和相關(guān)土壤特征性質(zhì)，假定土柱液相初始濃度分別為1.2 μg/L(Cd)、8 μg/L(Cu)、50 μg/L(Zn)、6 μg/L(Cr)。將實(shí)驗(yàn)數(shù)據(jù)導(dǎo)入inverse-solution模塊進(jìn)行反向演算，最終得到模擬結(jié)果與實(shí)驗(yàn)值的擬合曲線，如圖5～圖8所示。

表2 土壤水力特性參數(shù)

注：θs是土壤飽和含水率；θr是土壤殘余含水率；Ks是飽和導(dǎo)水率，cm/h；α為與進(jìn)氣吸力相關(guān)的參數(shù)，l/cm；n為土壤水分特征曲線的坡度大小影響系數(shù)；L為經(jīng)驗(yàn)擬合參數(shù)，一般為0.5。

圖5 Cd離子觀測(cè)值與模擬值

圖6 Zn離子觀測(cè)值與模擬值

圖7 Cu離子觀測(cè)值與模擬值

圖8 Cr離子觀測(cè)值模擬值

3.1.4 誤差分析

為了檢驗(yàn)?zāi)Ｐ偷哪M效果，運(yùn)用Hydrus-1d將實(shí)測(cè)數(shù)據(jù)與模擬數(shù)據(jù)進(jìn)行反演運(yùn)算，得到了能代表模型擬合程度的各項(xiàng)參數(shù)如表3所示。由表3可以看出，在各項(xiàng)離子擬合參數(shù)中R2均大于0.92，表示模型能夠模擬淋濾實(shí)驗(yàn)中各重金屬離子轉(zhuǎn)移量；均方根誤差(root-mean-square error，RMSE)是表示實(shí)際值與模擬值的偏差程度[14]，3種離子RMSE范圍為0.08～0.1，表示模型擬合程度較好；平均絕對(duì)誤差(mean absolute error, MAE)和平均誤差(mean error,ME)能更好的反映預(yù)測(cè)值誤差的實(shí)際情況，兩者參數(shù)均小于0.1，確保了模型的真實(shí)性。

表3 數(shù)據(jù)誤差參數(shù)

3.2 建立空間分布模型

由于影響土壤重金屬遷移因素眾多，利用常規(guī)手段對(duì)遷移特征進(jìn)行系統(tǒng)分析，雖能較好表述重金屬在一維或二維空間上的遷移特征，但在考慮外界環(huán)境因素對(duì)重金屬遷移影響方面相對(duì)欠缺。為了更加直觀的反映影響重金屬遷移因素，此次結(jié)合水文地質(zhì)中常用的空間分析方法對(duì)土壤重金屬遷移特征進(jìn)行分析?？臻g分布模型通常運(yùn)用于對(duì)區(qū)域單因子特征評(píng)價(jià)及來(lái)源解析；地下水建模系統(tǒng)(groundwater modeling system, GMS)是應(yīng)用廣泛的地下水水文模擬系統(tǒng)，在對(duì)污染物分布特征及遷移預(yù)測(cè)方面具有強(qiáng)大的空間計(jì)算能力[15]。基于前期對(duì)研究區(qū)野外地質(zhì)勘查，將研究區(qū)等高線提取到GMS中作出三維立體概化地質(zhì)模型(圖9)，由圖9可知，整個(gè)區(qū)域地形北高南低，地勢(shì)較低處為建筑帶及灌溉農(nóng)田，此次取樣點(diǎn)密度主要集中在農(nóng)田地帶；圖10為研究區(qū)地下水等水位線平面圖，流向至西北向東南，地下水水位與地形高程相關(guān)性較高。

為了進(jìn)一步分析研究區(qū)各重金屬的遷移及來(lái)源影響因素，在參照淋濾實(shí)驗(yàn)結(jié)果的基礎(chǔ)上，選取研究區(qū)所取土樣中的主要四種重金屬(Cu、As、Cr、Cd)離子作為研究對(duì)象，借助GMS中的泛克里金插值法對(duì)研究區(qū)內(nèi)重金屬空間分布特征展開分析，分析結(jié)果如圖11～圖14所示。

從圖11、圖12可知研究區(qū)域內(nèi)土壤中的Cr、Cu離子濃度沿西北向東南方向逐漸增加，呈規(guī)律分布，地勢(shì)較低的農(nóng)田地帶濃度相對(duì)偏高，造成此現(xiàn)象的主要原因可能是受長(zhǎng)期人為灌溉影響和淺層地下水水流方向上的路徑遷移所造成的，兩種離子受污染源影響不明顯；由圖13可知，As離子集中分布在研究區(qū)中上部及南部方向上，主要受區(qū)域內(nèi)污染源影響；圖14中Cd離子分布較為分散，局部呈現(xiàn)無(wú)明顯規(guī)律的富集現(xiàn)象，主要也集中在污染源附近，受水流遷移作用相對(duì)較弱。結(jié)合上述4種離子的空間分布特征，Cu、Cr離子與As、Cd離子之間存在著明顯的差異，而造成的這種差異的主要因素與區(qū)域內(nèi)土壤對(duì)各離子的吸附量有關(guān)，Cu、Cr離子在土壤中的含量明顯高于As、Cd離子，從淋濾試驗(yàn)也可以看出在土壤中含量越高的重金屬離子其在單位時(shí)間內(nèi)淋出量也較大，故在降雨作用和淺層地下水遷移作用下其濃度變化易向水流方向逐漸增大；相反對(duì)于土壤含量較小的重金屬元素，其分布受外界因素影響較大。

圖9 研究區(qū)三維地質(zhì)模

圖10 地下水等水位線圖

圖11 Cu離子空間分布圖

圖12 Cr離子空間分布圖

圖13 As離子空間分布圖

圖14 Cd離子空間分布

4 結(jié)論

(1)利用Hydrus-1D對(duì)研究區(qū)域建立水流溶質(zhì)運(yùn)移模型能較好的預(yù)測(cè)土壤淋濾實(shí)驗(yàn)中各重金屬遷移變化量，且擬合參數(shù)均在置信區(qū)間范圍內(nèi)。

(2)空間分布特征顯示土壤初始重金屬含量中占比較高的重金屬元素其遷移過(guò)程受淺層地下水水流作用較大；對(duì)于含量較低的重金屬元素，其更易受外界污染源影響，造成局部濃度偏高，分布特征無(wú)明顯規(guī)律。