基于數值模擬的土壤重金屬鉻污染預測分析

彭盼盼，伍靖偉，劉 聰，彭 翔，李 銀，劉 瑜

(1.武漢大學 水資源與水電工程科學國家重點實驗室，武漢 430072；2.天津市水利科學研究院，天津 300061)

隨著我國社會經濟的發展，土壤環境污染問題也日益突出。在造成我國土壤污染的各種污染物中，重金屬是重要的一種[1]，其中鉻已經成為主要的[2]土壤重金屬污染物之一。鉻元素在自然界中的主要賦存形態有三價鉻Cr(Ⅲ)和六價鉻Cr(Ⅵ)2種，前者是動物和人類必需的微量元素，后者對動植物和人類有很大的毒性，吸入后可能導致人體癌變或畸形。已有研究表明[3]，在土壤和地下水環境中，Cr(Ⅲ)和Cr(Ⅵ)在理化性質以及生物化學活性等方面存在著顯著的差異。Cr(Ⅲ)在土壤中具有較低的溶解度和流動性并且其自身生物毒性較弱；而土壤對Cr(Ⅵ)的吸附能力較弱，導致其在土體中的流動性較強。因此分析和預測土壤中Cr(Ⅵ)的運移轉化趨勢，對于改善作物種植環境、制定合理的環境保護和污染治理措施具有科學的指導意義。本文的目的是在野外監測鉻污染的基礎上運用數值模擬模型預測土壤Cr(Ⅵ的運移規律，為鉻污染防治提供科學依據。

1 研究區域概況及研究方法

1.1 研究區域概況

研究區位于天津市北辰區李嘴村，范圍覆蓋東經117°05′01.21″～117°06′45.25″，北緯39°13′55.38″～39°14′34.84″，北部邊界為新引河和永定新河支流，西部和南部邊界為北運河支流，東部邊界為外環西路及京津公路，總面積約1.87 km2。

研究區屬于海積沖積低平原地貌，整體地勢低平，局部由于人類活動影響有高差，海拔位于5～8 m。區內有農田、生活區、工業區、坑塘等，分別約占比40%、22%，37%，1%。該區域地下水埋深較淺，常年位于3 m以上，農田區地下水位稍高于工業區，地下水埋深在灌溉時期小于1 m。根據天津海濱工程勘察設計有限公司2015年5月完成的《引河橋李嘴村鉻污染調查項目土工試驗報告(工號：2015-063)》，研究區域埋深3 m深度范圍內土壤類型從上到下依次為人工填土層，厚度約為1.0 m；粉質黏土,厚度約為1.0 m；粉土，厚度約為1.0 m。

研究區鉻污染源來自于原天津市某化工廠鉻渣的堆放，該化工廠主要生產鉻鹽，停產后殘留了大量的工業廢渣。鉻渣露天堆放，受風吹和雨雪淋浸，所含的鉻隨降雨溶出滲入土壤和地下水中，造成污染環境。據相關部門的野外實測資料，目前該研究區污染源附近地下水中六價鉻的濃度含量高達700 mg/L。研究區域有大面積的農田，主要依靠抽取地下水進行灌溉，進而加劇了地下水中的污染物的運移，并且部分污染物滯留在土壤包氣帶，對土壤造成一定程度的污染。為保證居民的身體健康和糧食安全，有必要對土壤的鉻污染狀況進行預測。

1.2 研究方法

通過實地勘測，確定研究區地形地貌、鉻渣分布、農田種植結構、灌溉及居民用水情況以及污染物與周邊區域的聯系情況。在此基礎上，在研究區代表性位置布置觀測井對地下水位和地下水中污染物濃度進行動態監測，野外取樣對土壤鉻污染進行動態監測。通過實驗獲取土壤參數以及六價鉻在土壤中運移的相關物理參數，利用Hydrus-1D建立數值模擬模型，進而對鉻在土壤中的運移變化趨勢進行預測分析。

2 模型的建立和參數的設置

本研究運用Hydrus-1D軟件對土壤中Cr(Ⅵ)離子污染運移進行數值模擬模型的建立。Hydrus-1D[4]是由美國農業部、美國鹽堿土實驗室等機構開發的一款用于模擬一維飽和-非飽和土壤水分、熱量和溶質運移的專業軟件。

2.1 基本運動方程

Hydrus-1 D采用修正后的Richards方程[5]來描述水分在包氣帶中的運移過程，忽略包氣帶水水平和側向流動，只考慮垂向上的一維流動。取地面為坐標原點，z軸向上為正，Richards方程形式為：

(1)

式中：θ為體積含水率；h為土壤壓力水頭，cm；S為源匯項，h-1；K(h,z)是由以下方程給出的非飽和滲透系數：

K(h,z)=Ks(z)Kr(h,z)

(2)

式中：Kr是相對滲透系數，cm/s；Ks為飽和滲透系數，cm/s。

模型中采用經典對流-彌散方程描述一維溶質運移[6]，方程如下：

φ

(3)

式中：c為溶質濃度，g/cm3；s為吸附在土壤顆粒上的固態溶質濃度，g/cm3；D為飽和-非飽和水動力彌散系數；q為體積流動通量密度，cm/s；φ為匯源項(代表溶質發生反應的各種零級、一級其他反應)，g/cm3。

2.2 模型基本參數的確定

HYDRUs-1D模型所需的驅動數據包括氣象數據、土壤水力特性參數以及溶質運移和吸附等參數，各自獲取方法如下。

大氣數據主要用于估算作物需水量和農田灌溉用水量[7]，包括逐日最低氣溫、最高氣溫、濕度、日照時數、降雨量、風速等彭曼-蒙特斯公式需要的數據。本研究根據中國氣象科學數據共享服務網提供的數據進行本研究。

本研究采用Van-Genutchen水力特性模型表達土壤水力特性參數之間的關系。Van-Genutchen模型[8]如下:

(4)

K(h)=KsSle[1-(1-S1/m3)m]2

(5)

m=1-1/nn>1

(6)

(7)

式中：θs為土壤飽和體積含水率，cm3·cm-3；θr為土壤殘余體積含水率，cm3·cm-3；h為土壤負壓，cm；α為進氣值的倒數，cm-1；Ks為飽和水力傳導度，cm/d；Se為有效飽和度，無量綱；l為土壤孔隙連接參數，無量綱，模型中取l=0.5；m與n為土壤孔隙尺寸分布參數，無量綱。

根據現場勘測取樣和作物根系可能達到的最大深度，取土壤模擬深度為200 cm。剖分時采用一個1 cm為空間步長，共剖分201個節點，共設定5個觀測點，根據地質勘測提供的資料將土壤剖面劃分為2個土壤介質層，0～100 cm為粉質黏土層，100～200 cm為黏土層。根據地質勘探給出的地質資料和天然密度，利用Hydrus-1D上的神經網絡預測功能[8]，得到天然容重下實驗土料的土壤水分運動參數，見表1。

表1 土壤水分運動參數Tab.1 Soil hydrodynamic parameters

注：表中各參數符號意義與Van-Genutchen模型中相應參數一致。

土壤吸附參數由土柱實驗擬合得到，粉土的吸附分配系數取0.279 8 L/kg；粉質黏土的吸附分配系數取0.337 9 L/kg。

2.3 模型的建立

根據研究區的土地利用現狀，將該區域分為農田區、污染源區和生活區3個子區域分別建立模型進行污染物運移的數值模擬，見圖1。

圖1 研究區域分區Fig.2 Subranges of the study area

對于污染源區和生活區，由于模擬土壤剖面地表直接暴露于大氣中，接受大氣降雨補給，且研究區域有排水設施，降雨可以產生徑流，因此上邊界為有徑流的大氣邊界，通過對研究區地下水位的動態觀測可知地下水位波動不大，因此下邊界條件設為變水頭邊界。溶質運移的上邊界條件選擇濃度通量邊界，下邊界選擇零濃度梯度邊界。上邊界根據降雨量、降雨濃度(設為0)求得。

由于污染源區和生活區不進行灌溉，因此上邊界氣象條件只考慮降雨，此外該區基本為水泥地或屋頂覆蓋，地表蒸發量甚微，所以不考慮蒸發。初始條件，根據實際測得的土壤含水量和土壤鉻濃度設定。

對于農田區，該研究區域主要作物為玉米和白菜，因此在模擬時需要考慮灌溉和作物吸收對非飽和帶土壤中污染物運移變化的影響。此灌溉或降雨可產生積水，所以上邊界選取可積水的大氣邊界，下邊界設選取變化的水頭邊界。該區域的溶質運移模塊上邊界設置為濃度通量邊界，下邊界為零濃度梯度邊界。上邊界根據灌溉水量、灌溉水濃度(即地下水中的濃度)求得。

2.4 模型的標定和驗證

標定驗證的數據來自天津市水科院2015年開展的室內土柱污染物運移實驗。實驗用土來自研究區未受污染的粉土和粉質黏土，土樣經過風干后填裝土柱。土柱高30 cm，由馬氏瓶控制并保持5cm固定入滲水頭，溶液中污染物鉻的初始濃度為550 mg/L, 對土柱進行淋洗，并在土柱最下端口(30 cm)處采集水樣，測出滲出液中的鉻濃度。將模擬結果和實測結果進行對比，見圖2和圖3。

圖2 粉土滲出液中Cr(Ⅵ)濃度模擬值與實測值比較Fig.2 The comparison of Cr(Ⅵ)in leachate of silt between simulation and measured value

圖3 粉質黏土滲出液中污染物濃度模擬值與實測值比較Fig.3 The comparison of Cr(Ⅵ)in leachate of silt clay between simulation and measured value

從以上2圖可以直觀地看出，滲出液中Cr(Ⅵ)的濃度隨時間先增加，后趨于穩定狀態，說明2種土壤對Cr(Ⅵ)的吸附達到平衡。2種土壤土柱污染試驗的模擬值與實際觀測值差距很小，污染物濃度的總體變化趨勢一致，說明該模型可以較為準確的動態模擬土壤中污染物的運移變化，預測土壤中污染物濃度的變化趨勢。

3 土壤中Cr(Ⅵ)污染的預測分析

以Hydrus-1D的非飽和水流模型為基礎，建立數值模型，預測未來污染物Cr(Ⅵ)離子在淺層土壤中的運移變化趨勢。模擬時段從2015年3月開始到未來30 a，以天為計算單位。考慮到氣象變異較強烈，且對土壤水和地下水的影響較大，因此本研究采用長系列氣象資料進行模擬。采用1984-2014年的氣象數據作為模型驅動，計算作物ET[9]、灌溉水量、蒸發、淋洗、降雨入滲等。

3.1 污染源區土壤中Cr(Ⅵ)污染趨勢

由圖4可知，在土壤剖面上Cr(Ⅵ)分布并無明顯規律，5 a后Cr(Ⅵ)濃度分布趨于穩定狀態。由圖5可見在模擬初期0～3 a，不同深度的土層中Cr(Ⅵ)含量均有明顯的減少，3～30 a可能由于氣象條件、降雨入滲蒸發等因素的影響，各深度土層中Cr(Ⅵ)的含量均呈周期性變化。根據GB15618-2008土壤環境質量標準，工業工地土壤無機物污染的二級標準為30 mg/kg，由圖5可以看出，在污染源區在初始時段到1 a左右，土壤剖面污染物濃度超標，預測期5 a之后污染物濃度下降至土壤無機物的環境質量二級標準值以下。

圖4 污染源土壤剖面Cr(Ⅵ)濃度分布Fig.4 The distribution of Cr(Ⅵ)in soil profile of contaminant source area

圖5 污染源不同深度土層Cr(Ⅵ)濃度隨時間變化Fig. 5 The change of Cr(Ⅵ)in different depth of soil vary time in contaminant source area

3.2 生活區土壤中Cr(Ⅵ)污染趨勢

由圖6可見在土壤剖面上Cr(Ⅵ)分布并無明顯規律，5 a后Cr(Ⅵ)濃度分布趨于穩定狀態。由圖7可見，在模擬初期0～3 a，不同深度的土層中Cr(Ⅵ)含量均有明顯的減少，3～30 a間可能由于氣象條件、降雨入滲蒸發等因素的影響，各深度土層中Cr(Ⅵ)的含量均呈周期性變化。根據GB15618-2008土壤環境質量標準，居住地土壤無機物污染的二級標準為5 mg/kg，由預測結果可以看出，在生活區土壤中Cr(Ⅵ)的含量并沒有超出土壤無機物污染的二級標準值。

圖6 生活區土壤剖面Cr(Ⅵ)濃度分布Fig.6 The distribution of Cr(Ⅵ)in soil profile of living-area

圖7 生活區不同深度土層Cr(Ⅵ)濃度隨時間變化Fig.7 The change of Cr(Ⅵ)in different depth of soil vary time in living-area

由于鉻渣堆已經遷走，在污染源和生活區沒有鉻的補給來源；加之土壤中的有機質、鐵氧化物將Cr(Ⅵ)還原為Cr(Ⅲ)，成為沉淀或絡合物[3,10]；在降雨入滲淋洗作用下，部分Cr(Ⅵ)隨水分的運動進入飽和帶含水層中[2]，所以3～5 a內淺層土壤中Cr(Ⅵ)的含量明顯減少。一部分Cr(Ⅵ)被土壤顆粒吸附，并達到吸附平衡[11]，所以5 a之后含量呈周期性變化。

3.3 農田區土壤中Cr(Ⅵ)污染趨勢

由圖8可見在在土壤剖面0～100 cm，Cr(Ⅵ)隨深度的增加而增多，可能是受作物根系[12]分布的影響；100～200 cm處于作物根系層一下，分布并無明顯規律，1 a后Cr(Ⅵ)濃度分布趨于穩定狀態。由圖9可見，在模擬初期，不同深度的土層中Cr(Ⅵ)含量均有明顯的減少，1 a以后可能由于降雨入滲、農田灌溉、蒸發、作物吸收等因素的綜合影響，使農田區域土壤中Cr(Ⅵ)的含量一直處于比較低的狀態。由于現階段對土壤中六價鉻暫無同一的評價標準，本文按照GB15618-2008土壤環境質量標準，居住地土壤無機物污染的二級標準對農田區進行評價。由預測結果可見，農田區土壤中Cr(Ⅵ)的含量遠遠低于5 mg/kg，處于無污染狀態。

圖8 農田區土壤剖面Cr(Ⅵ)濃度變化Fig.8 The distribution of Cr(Ⅵ)in soil profile of farmland

圖9 農田區不同深度土層Cr(Ⅵ)濃度隨時間變化Fig.9 The change of Cr(Ⅵ)in different depth of soil vary time on farmland

4 結 語

(1)本研究在獲取了土壤參數和Cr(Ⅵ)在土壤中的吸附分配系數的基礎上，運用Hydrus-1D建立了Cr(Ⅵ)在淺層土壤中的運移轉化數值模擬模型。

(2)運用建立的運移轉化數值模擬模型對未來30 a Cr(Ⅵ)離子在研究區域淺層土壤中的運移轉化趨勢進行預測。

(3)由模型預測結果可知，預測期5 a之后，污染源處土壤中Cr(Ⅵ)的含量低于30 mg/kg，生活區土壤中Cr(Ⅵ)的含量低于5 mg/kg，根據GB15618-2008土壤環境質量標準，均低于土壤無機物污染的二級標準。

□

[1] 盧曉華. 基于數值模擬的企業地下水重金屬污染的環境影響預測評價 [J]. 安全與環境工程, 2014,21(1):93-7.

[2] 李喜林. 鉻渣堆場滲濾液對土壤-地下水系統污染規律研究[D]. 遼寧阜新：遼寧工程技術大學, 2012.

[3] 傅臣家, 劉洪祿, 吳文勇,等. 六價鉻在土壤中吸持和遷移的試驗研究 [J]. 灌溉排水學報, 2008,27(2):9-13.

[4] Genuchten BV. HYDRUs-1D code for simulating the one-dimensional movement of water, heat, and multiple solutes in variably saturated media[C]∥ proceedings of the International Ground Water Modeling Center, Colorado School of Mines, F, 2010.

[5] Celia Ma BE, Zarba Rl. A general mass conservation numerical solution for the unsaturated flow equation [J]. Water Resour Res, 1990,(26):1 483-1 496.

[6] Sudo H, Kawamata T, Matsuoka Y, et al. The HYDRUs-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media [J]. Dep. of Environmental Sciences Univ of California Riverside, 2008,68.

[7] Pereira LS, Gowing JW. Water and the environment: innovation issues in irrigation and drainage[J]. Water & the Environment Innovation Issues in Irrigation & Drainage, 1998,25(9):87.

[8] Radcliffe DE,imnek J. Soil physics with HYDRUS: modeling and applications[M]. CRC, 2010.

[9] Smith M, Group A. CLIMWAT for CROPWA: a climatic database for irrigation planning and management[M]. Food Agriculture Organization of the United Nations, 1993.

[10] 雷建森. 六價鉻在土壤中的吸附特性及風險評價研究[D]. 長春：吉林大學, 2015.

[11] 李喜林, 王來貴, 郝 喆,等. 粉質黏土和粉質砂土對鉻渣滲濾液中Cr(Ⅵ)吸附特性[J]. 環境工程學報, 2013,7(12):5 019-5 024.

[12] 潘 崠, 丁紹蘭. 鉻在土壤和植物中的吸附、遷移和累積規律研究進展[J]. 西部皮革, 2011,33(2):22-25.