基于数值模拟的土壤重金属铬污染预测分析

彭盼盼，伍靖伟，刘 聪，彭 翔，李 银，刘 瑜

(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.