万朔阳， 吴 勇， 唐学芳， 邓东平， 兰 真， 韩莉璧

(1.成都理工大学环境与土木工程学院，成都 610059；2.地质灾害防治与地质环境保护国家重点实验室，成都 610059)

随着科学技术水平的提高,农业技术也不断得到改进, 人类活动与土壤之间的联系日益增加，即来对土壤污染也逐渐加剧，现如今因土壤盐碱化, 化肥、农药污染、污水灌溉引起的地下水和土壤污染等问题日益受到关注。探究土壤溶质运移规律逐渐成为解决此类问题的关键所在，由于影响土壤溶质运移因素众多仅靠室内试验难以掌控其实际运移规律。故通过对受污染场地土壤建立数学模型来处理这些问题已经成为众多学者的热门研究对象[1-2]。早在20世纪90年代Simunek等[3]就已利用Hydrus-1D对非饱和土的水力和溶质输运参数进行估算，验算出各参数的控制方程，为分析土壤溶质运移特征提供了重要的理论基础。郭瑞等[4]在此基础上以水盐运移理论和环境条件为依据，建立了适宜于不同条件下的水盐运移模型。王锐等[5]在分析入渗压力水头对土壤水分入渗影响的基础上，建立了不同入渗水头条件下土壤水分运动数学模型，模拟值与实验值最大误差不超过6%。刘目兴等[6]对比研究了黏质土壤在不同初始含水率条件下对初始入渗率的影响，结果表明土壤初始含水率与初始入渗率成反比，并且通过构建Horton模型的对结论进行拟合模拟，拟合效果较好。李茜等[7]提出构建土壤水与地下水溶质运移联合模型，将两者相结合，更好的描述区域溶质迁移的规律。由上述可知，对土壤建立溶质运移模型不仅能对运移过程的分析提供理论依据，同时也能较好预测土壤污染物迁移规律，对土壤修复及防治有着重要意义。

选取四川西坝镇某高产农田为研究区域，对所取土样进行相关室内实验参数分析的基础上，运用Hydrus-1D建立水流溶质运移模型，对恒定水流条件下土壤室内淋滤实验结果进行数据拟合；并且借助地下水建模系统(groundwater modeling systems,GMS)构建空间分布模型，分析各项重金属在研究区范围内的空间分布特征，对研究区内土壤重金属迁移规律及分布特征差异作出总体的评述。

1 研究区概况

研究区位于四川省乐山市五通桥区南西(方位253°)，直距约4.0 km处，面积2.4 km2(图1)；此次共取72组土样对研究区内的高产农田重金属迁移及分布特征展开分析，各取样点如图2所示。

2 材料与方法

2.1 土样实验方法

对研究区所取72组土样采用高温密闭消解法对土壤样品进行消解。采用ICP-MC对土壤消解产物进行重金属测试分析。重金属检测项目主要包括：Cd、Pb、Cr、Cu、Hg、As等，主要选取Cd、Cu、Zn、As、Cr等离子对区域内土壤重金属分布及迁移特征展开讨论，实验数据统计见表1，均在成都理工大学地质灾害防治与地质环境保护国家重点实验室测得。

图1 研究区位置及地貌

图2 研究区采样点分布

表1 土壤重金属数据统计表

2.2 室内淋滤实验方法

土柱实验所用土样取至研究区某块高产农田深度为1 m的探坑内，沿地表向下1 m，每隔10 cm取若干组土样；实验土柱是用厚5 mm、长20 cm，直径为10 cm有机玻璃管制作而成，如图3和图4所示，柱顶连接供水管、柱底有渗漏孔，土柱共分为四层，按照探坑不同深度取样依次分层填样，第一层(0～5 cm)填充探坑0～20 cm所取土样，依次填充，第四层(15～20 cm)为探坑80～100 cm所取土样。实验采用马氏瓶稳定供水，淋滤液为超纯水，淋滤时间为11天，每隔一天取一次渗漏水样，记录其中Cd、Cu、Zn、Cr离子含量变化值。

图3 淋滤实验装置

3 结果及分析

3.1 建立溶质运移模型

Hydrus-1D软件[8]是一套用于模拟饱和-非饱和多孔介质中水分运移和溶质运移的数值模型，在对非饱和带土壤水溶质运移模拟方面，具有操作简单，拟合效果较好等优势，因此广泛被运用于此类问题的研究。通过确立溶质基本控制方程 、边界条件、土壤水力特征参数等过程，来构建溶质运移模型。

3.1.1 模型基本控制方程

采用Hydrus-1D软件，对研究区土壤垂向上建立一维水流溶质运移模型，此次水流模型概化为多层均质各项同性的一维垂向稳定流[9]，土壤水分运动主要受Richards方程控制，土壤溶质运移主要受对流-弥散方程控制[10]。

(1)

(2)

式中：θ为土壤体积含水率，cm3/cm3；h为土壤压力水头，cm；S为根系吸水项或其他源汇项，cm3/cm3/T(T为时间单位常量，可代表h、min、s)；K为非饱和土壤水力梯度传导，cm/T；KAij为各项异性张量的分量，(-)；Xi为空间坐标，cm；Z为垂项坐标(向上为正)，cm；t为时间，T；c为溶质浓度，g/cm3；D为水动力弥散系数，cm2/h；q为达西流速，cm/h。

3.1.2 初始条件与边界条件

本次模拟条件在不考虑表面蒸发及植物根系吸收作用下进行，将初始条件设定为土柱中各液相浓度分布，对于水流模型上边界为通量边界，下边界为自由渗漏边界；对于溶质运移模型上边界为大气通量边界，下边界为零梯度边界[11]。

(1)土壤水分运移的初始条件与边界条件。

初始条件：

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

(3)

上边界条件：

(4)

下边界条件：

h(x,t)=hb(t)

(5)

式中：h0(x)表示土壤剖面初始土壤水负压,cm;hb(t)表示下边界压力水头,cm；qs表示地表水入渗量。

(2)土壤溶质运移的初始条件与边界条件。

初始条件：

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

(6)

上边界条件：

(7)

下边界条件：

(8)

式中：c0为初始剖面土壤液相浓度,g/cm3；cs为上边界流量溶质浓度，此次试验采用超纯水故浓度cs=0。

3.1.3 土壤水力特征参数

土壤水力特征参数属于土壤的基本物理性质,其各项参数属性直接影响土壤水在包气带中的迁移转化关系,对分析土壤溶质运移起重要作用[12]；通常采用van-Genuchten模型描述土壤水力特征参数[13]。为了将模型更好地与室内实验结果相拟合，使模拟过程更具真实性，基于实验土壤不同粒径分选占比量和土壤容重数据，利用Hydrus-1D内置的预测及反演土壤水力特征参数功能，并结合野外场地实验值对预测值进行校正，最终推测出本研究区内土壤的水力特征参数(表2)。模拟时间设为11 d，土柱分为2层，上层为砂质壤土、下层为黄褐土粉砂壤土，虚拟观测点位于土柱低端渗水口处，根据土壤自身所含重金属含量和相关土壤特征性质，假定土柱液相初始浓度分别为1.2 μg/L(Cd)、8 μg/L(Cu)、50 μg/L(Zn)、6 μg/L(Cr)。将实验数据导入inverse-solution模块进行反向演算，最终得到模拟结果与实验值的拟合曲线，如图5～图8所示。

表2 土壤水力特性参数

注：θs是土壤饱和含水率；θr是土壤残余含水率；Ks是饱和导水率，cm/h；α为与进气吸力相关的参数，l/cm；n为土壤水分特征曲线的坡度大小影响系数；L为经验拟合参数，一般为0.5。

图5 Cd离子观测值与模拟值

图6 Zn离子观测值与模拟值

图7 Cu离子观测值与模拟值

图8 Cr离子观测值模拟值

3.1.4 误差分析

为了检验模型的模拟效果，运用Hydrus-1d将实测数据与模拟数据进行反演运算，得到了能代表模型拟合程度的各项参数如表3所示。由表3可以看出，在各项离子拟合参数中R2均大于0.92，表示模型能够模拟淋滤实验中各重金属离子转移量；均方根误差(root-mean-square error，RMSE)是表示实际值与模拟值的偏差程度[14]，3种离子RMSE范围为0.08～0.1，表示模型拟合程度较好；平均绝对误差(mean absolute error, MAE)和平均误差(mean error,ME)能更好的反映预测值误差的实际情况，两者参数均小于0.1，确保了模型的真实性。

表3 数据误差参数

3.2 建立空间分布模型

由于影响土壤重金属迁移因素众多，利用常规手段对迁移特征进行系统分析，虽能较好表述重金属在一维或二维空间上的迁移特征，但在考虑外界环境因素对重金属迁移影响方面相对欠缺。为了更加直观的反映影响重金属迁移因素，此次结合水文地质中常用的空间分析方法对土壤重金属迁移特征进行分析。空间分布模型通常运用于对区域单因子特征评价及来源解析；地下水建模系统(groundwater modeling system, GMS)是应用广泛的地下水水文模拟系统，在对污染物分布特征及迁移预测方面具有强大的空间计算能力[15]。基于前期对研究区野外地质勘查，将研究区等高线提取到GMS中作出三维立体概化地质模型(图9)，由图9可知，整个区域地形北高南低，地势较低处为建筑带及灌溉农田，此次取样点密度主要集中在农田地带；图10为研究区地下水等水位线平面图，流向至西北向东南，地下水水位与地形高程相关性较高。

为了进一步分析研究区各重金属的迁移及来源影响因素，在参照淋滤实验结果的基础上，选取研究区所取土样中的主要四种重金属(Cu、As、Cr、Cd)离子作为研究对象，借助GMS中的泛克里金插值法对研究区内重金属空间分布特征展开分析，分析结果如图11～图14所示。

从图11、图12可知研究区域内土壤中的Cr、Cu离子浓度沿西北向东南方向逐渐增加，呈规律分布，地势较低的农田地带浓度相对偏高，造成此现象的主要原因可能是受长期人为灌溉影响和浅层地下水水流方向上的路径迁移所造成的，两种离子受污染源影响不明显；由图13可知，As离子集中分布在研究区中上部及南部方向上，主要受区域内污染源影响；图14中Cd离子分布较为分散，局部呈现无明显规律的富集现象，主要也集中在污染源附近，受水流迁移作用相对较弱。结合上述4种离子的空间分布特征，Cu、Cr离子与As、Cd离子之间存在着明显的差异，而造成的这种差异的主要因素与区域内土壤对各离子的吸附量有关，Cu、Cr离子在土壤中的含量明显高于As、Cd离子，从淋滤试验也可以看出在土壤中含量越高的重金属离子其在单位时间内淋出量也较大，故在降雨作用和浅层地下水迁移作用下其浓度变化易向水流方向逐渐增大；相反对于土壤含量较小的重金属元素，其分布受外界因素影响较大。

图9 研究区三维地质模

图10 地下水等水位线图

图11 Cu离子空间分布图

图12 Cr离子空间分布图

图13 As离子空间分布图

图14 Cd离子空间分布

4 结论

(1)利用Hydrus-1D对研究区域建立水流溶质运移模型能较好的预测土壤淋滤实验中各重金属迁移变化量，且拟合参数均在置信区间范围内。

(2)空间分布特征显示土壤初始重金属含量中占比较高的重金属元素其迁移过程受浅层地下水水流作用较大；对于含量较低的重金属元素，其更易受外界污染源影响，造成局部浓度偏高，分布特征无明显规律。