刘欢， 孙洪广

(河海大学力学与材料学院， 水文水资源与水利工程科学国家重点实验室， 南京 210098)

20世纪及21世纪初，由于生产条件落后以及环保意识薄弱，许多工矿企业(特别是化工企业)在生产过程中对附近土壤和地下水造成了严重的重金属污染[1-2]。近年来，中国对生态环境问题越来越重视，老旧工业企业和废弃场地的重金属污染修复工程越来越多。作为土壤和地下水修复的基础，数值模拟有助于掌握地下水溶质运移规律，了解污染物迁移的趋势，确定污染范围及浓度分布。因此，数值模拟方法已被广泛应用于污染场地的定量评估与准确预测，对水土污染管控和高效修复提供科学指导。

目前，关于土壤中重金属迁移规律的研究已有大量报道。已有研究主要集中在两方面：一是通过室内土柱淋溶与砂箱实验，结合理论分析考察重金属在土壤中的迁移规律；二是通过数值方法对大型场地重金属迁移进行模拟与预测。江建明等[3]通过砂箱物理实验与数值模拟发现裂隙对重金属污染物迁移有导向作用。王玉雪等[4]通过室内土柱实验得到初始pH越低或初始浓度越高，重金属Cr(Ⅵ)穿透时间与达到吸附饱和时间越早的结论。林青等[5]通过室内置换实验对铜、镉、锌、铅在石英砂中的迁移规律进行了研究，并通过数值模拟得出静态吸附实验得到的吸附参数不能直接应用于重金属运移模拟。Opekunova等[6]以某铜矿为研究对象，分析了预测铜、锌、镍及镉等重金属元素在土壤中的迁移。辛圆心等[7]针对狮子山尾矿，用COMSOL软件考察了铜和锌离子在土壤中的迁移规律。Lafratta等[8]通过对澳大利亚某铅锌冶炼厂的土壤进行多参数研究，重现重金属迁移过程，证实该冶炼厂仍是土壤重金属污染的主要来源。Dong等[9]通过数值模拟重金属在垃圾填埋场的迁移行为，得出在土壤与地下含水层最敏感的参数是阻滞因子与渗透系数。Zhang等[10]分析了铁尾矿中重金属的迁移特征及分布规律。孙洪广等[11]研究表明，分数阶导数模型能很好的描述溶质在土壤环境中的反常扩散行为。

基于此，首先以南京某废弃化工厂建立三维地下水重金属溶质运移模型，结合有限元方法，对该土壤中重金属砷[As(V)]的浓度演化进行模拟及预测；其次，基于该模型考察不同渗透系数、弥散度、分配系数以及孔隙率的条件下，重金属As(V)在土壤中的迁移规律；最后，引入敏感性指数，比较各参数变化对重金属As(V)溶质迁移的影响。为类似地区的水土修复与地下水环境影响评价工作提供科学依据。

1 数学模型的建立

1.1 潜水层水流运动数学模型

(1)

式(1)中：h为地下水水头；kxx、kyy、kzz分别为x、y、z方向的渗透系数；w为水流源汇项，包括降水入渗、河流入渗补给、井的抽水量等；Ss为给水系数；H0为含水层初始水位；H1为定水头边界条件；q为第二类边界的单宽流量，流入为正，流出为负，隔水边界为0；t为时间；q1为定流量边界条件；Γ1为模型定水头边界；Γ2为垂向定流量边界。

1.2 重金属溶质运移数学模型

i,j=x,y,z;xx=x;xy=y;xz=z

(2)

所研究的重金属砷在土壤中的迁移，只考虑对流、弥散以及吸附作用，引入阻滞因子R，其表达式为

(3)

(4)

式中：C为土壤溶质浓度；θ为土体孔隙率；ρ为土壤干体积质量即容重；Q为溶质源汇项；Dij为弥散系数；qi为达西流速；t为时间；Rn为土壤溶质各个转化过程包括生物吸收、化学反应、降解、沉淀等过程；Kd为重金属溶质吸附相与溶解相的平衡分配系数。

2 几何建模与边界条件

主要研究中国南京某污染地块的重金属As(V)在地下水中的运移规律以及随时间变化的浓度分布状态。通过对该地块的实地调查，建立重金属在土壤中运移的三维模型。该模型可概化为250 m×160 m的矩形地块(图1)，潜水含水层约为10 m。东西南北四面皆临近河流，地下水流大致为自东向西流动。可将东面和西面设定为恒定水头，南北方向为沿程均匀下降的定水头边界。垂向边界系统，上边界为潜水含水层自由水面，系统与外部垂向水量交换都是通过该边界完成的，如大气降水补给，研究区域下部处理为零通量边界。污染源为该废弃产地的焚烧炉范围，通过大气降水淋溶进入地下水。

为了考察各参数对重金属As(V)在土壤中迁移的影响，以该模拟案例为基准，将孔隙率、渗透系数、弥散度、吸附分配系数分别上下波动25%、50%作为4个水平值，连续通入10年污染物，研究分析各种参数变化对溶质运移的影响。并采用因子变换法[12]，以敏感性指数为量化指标对各参数进行分析。具体参数如表1所示。

图1 三维几何建模图Fig.1 3D geometric modeling drawing

表1 各参数水平值Table 1 Parameter level values

采用以下敏感性指数计算公式对模拟结果进行敏感性评价，L值越大，表明参数改变对污染物迁移距离影响越大，即该参数的敏感性越高。计算公式为

(5)

式(5)中：L为敏感性指数，无量纲；n为参数水平值个数，取值为4；cj为对照组重金属As(V)纵向迁移距离；c′j为其他水平值重金属迁移距离。

3 计算结果及分析

3.1 重金属砷溶质迁移模拟结果

该化工厂生产含有重金属As(V)产品满7年后停止运营，随后半年完成搬迁。

图2 重金属砷浓度分布随时间变化Fig.2 Variation diagram of arsenic concentration distribution with time

图2为土壤中重金属砷的分布随时间的变化，可以看出，前7年污染源重金属As(V)持续泄露。污染晕中心以及污染羽中As(V)的质量浓度逐年升高，预测时段分别为2、5、7年时，污染晕中心浓度达到2.25、5.55、7.12 mg/L；按照《地下水质量标准》(GB/T 14848—2017)，重金属砷标准值为0.05 mg/L，浓度分布如图2所示，污染面积随着时间推移不断增加，超标范围也不断增加。

7年后，由于该工厂搬迁，污染面源处不再流入污染物。在地下水对流作用下，污染范围整体向下游迁移，在第20年、40年时，污染晕中心位置往下游分别迁移了2.7 m与6.5 m；弥散作用使得污染晕中心浓度向四周扩散，污染晕中心浓度逐渐减少，污染范围增加明显，第40年时污染晕中心浓度已下降到0.33 mg/L。

同时可以注意到重金属As(V)污染物在地下水中的迁移并不是向四周均匀扩散的，污染羽纵向迁移距离更远，横向迁移距离较短，在第7年时，污染羽纵向迁移距离达38.5 m，横向迁移距离仅23.7 m。在垂向浓度主要集中在土壤表层，污染物浓度由表层向下减少，第20年超标深度才达1.5 m，这是因为纵向是水力梯度的主要方向，该方向对流以及弥散系数值都比其他方向大。

3.2 模拟结果对不同参数的响应

为研究分析各水文地质参数对重金属As(V)溶质迁移的规律，在上述模型的基础上，只考虑连续通入10年污染物的情况。并且以污染源中心及纵向浓度变化为主要研究对象，分析不同参数变化，重金属As(V)溶质迁移对其的响应规律。

3.2.1 渗透系数对重金属砷溶质迁移的影响

当渗透系数增加时，重金属As(V)溶质在纵向达到更大的迁移距离；另外，渗透系数也影响着溶质在纵向分散上的传输，污染晕中心浓度随着水力渗透系数的增加而减小。图3为其他参数不变，渗透系数K分别取2.4×10-5、3.6×10-5、4.8×10-5、6.0×10-5、7.2×10-5m/s时，重金属砷通入第10年的纵向浓度变化。渗透系数取值为2.4×10-5m/s时，污染晕中心浓度达到11.2 mg/L，相对参照组增幅达到29.2%，向下游迁移距离减少了6 m；渗透系数取值为7.2×10-5m/s时,污染晕中心浓度降低到7.1 mg/L，重金属As(V)向下游迁移距离可达24.5 m，相比参照组污染源中心浓度降幅为19.1%，往下游迁移距离增加了4 m。这是因为在水头差一定的情况下，一方面，渗透系数的增加直接导致地下水达西速度及渗流速度的增加，污染物随地下水流运动，加快了溶质向下游的迁移，从而污染晕中心浓度降低，向下游迁移距离增加；另一方面纵向弥散系数也相应增大，随着渗透系数增加，图3中重金属As(V)往水流上游迁移距离也呈增加趋势。这也侧面说明，污染物在吸附力度较强的土壤中，弥散作用对重金属As(V)的迁移影响不可忽视。

图3 不同渗透系数重金属砷纵向浓度变化曲线Fig.3 Longitudinal concentration curves of heavy metal arsenic with different hydraulic conductivity

3.2.2 弥散度对重金属砷溶质迁移的影响

弥散度是用来描述含水层中污染物弥散作用的参数，弥散系数除以孔隙平均流速便是弥散度，其取值与土的孔隙结构有关，且具有明显的尺度效应[13-14]。当弥散度增加时，重金属As(V)溶质在土壤中迁移的变化趋势与渗透系数增加一致，重金属As(V)溶质纵向迁移距离增加，污染源中心浓度减小。图4为其他参数不变，αL分别取1.750、2.625、3.500、4.375、5.250 m时，重金属溶质通入第10年的纵向浓度变化。纵向弥散度取1.750 m时，污染源中心浓度为11.14 mg/L，与对照组弥散度取3.5 m相比，污染源中心浓度增加了28.4%，往下游迁移距离减少了4.5 m。纵向弥散度取5.250 m时，污染源中心浓度为7.13 mg/L，向下游迁移距离为23.5 m，与对照组相比，中心浓度减少17.8%，往下游迁移距离增加3.5 m。这是因为在地下水流速一定情况下，随着弥散度增大，弥散系数也相应增大，土壤中重金属As(V)溶质受到弥散作用越强，从而污染源溶质向四周弥散越快，表现为污染源中心浓度越低，迁移距离越远，该数值模拟结果上游以及下游都是如此情况。

图4 不同弥散度重金属砷纵向浓度变化曲线Fig.4 Longitudinal concentration curve of heavy metal arsenic with different dispersion

3.2.3 分配系数对重金属砷溶质迁移的影响

图5 不同分配系数重金属砷纵向浓度变化曲线Fig.5 Longitudinal concentration curves of heavy metal arsenic with different partition coefficients

图5为其他参数不变，分配系数分别取10、15、20、25、30 L/kg，重金属溶质通入第10年的纵向浓度变化。结果表明，重金属As(V)在土壤中的迁移距离随吸附分配系数的增加而减小，污染晕中心浓度与分配系数的大小也呈负相关变化。Kd为30 L/kg时，运移10年污染源中心浓度6.9 mg/L，向下游迁移17 m，相对参照组Kd=20 L/kg，污染源中心浓度降幅27.6%，向下游迁移距离减少了3 m；Kd为10 L/kg时，污染晕中心浓度为11.6 mg/L，向下游迁移28 m，与参照组相比，污染源中心浓度增幅34.4%，迁移距离增加了8 m。这是由于分配系数越大，土壤对污染物的吸附能力越强，地下水中越多重金属As(V)溶质被吸附到土壤中，导致地下水中As(V)浓度减少；土壤对重金属As(V)污染物吸附增强，也表现为砷在土壤中受到的阻滞作用越强，溶质在土壤中的迁移速度越慢，因而迁移距离越短[9]。

3.2.4 孔隙率对重金属砷溶质迁移的影响

孔隙率对溶质迁移的影响在两方面，一是孔隙率直接决定渗流速度，而渗流速度控制对流迁移；二是决定模型单元中储存溶质的孔隙体积大小。在给定渗透系数情况下，即达西流速不变，孔隙率对溶质的运移作用主要表现在污染物浓度分散上面，孔隙率越大，孔隙实际平均流速越小，水动力弥散系数也越小，污染晕中心的浓度越高，污染物迁移距离越短[15]。在分配系数较小或者可以忽略的情况下的确如此，然而对于土壤吸附非常强的重金属As(V)而言，污染物并非以水的渗流速度v做对流，其孔隙流速为v/R，对应的弥散系数D也是如此，即D/R[16]。在0～1改变θ值，对污染物渗流速度以及弥散度几乎没有影响。图6为其他参数不变，孔隙率分别为0.210、0.315、0.420、0.525、0.630时，重金属砷通入第10年的纵向浓度变化。图6中孔隙率θ在0.21～0.63变化，污染晕中心浓度相对于参照组变化幅度为-0.04～0.04 mg/L，迁移距离几乎无差别。由此可见，在土壤对污染物吸附力度较强的土壤中，相对来说孔隙率的变化对污染物迁移的影响较弱。

图6 不同孔隙率重金属砷纵向浓度变化曲线Fig.6 Longitudinal concentration curve of heavy metal arsenic with different porosity

3.2.5 各参数对重金属溶质迁移影响的比较

各参数变化对污染物迁移距离的影响如表2所示，可以看出，在0～50%波动范围，渗透系数、弥散度以及吸附分配系数对重金属As(V)迁移距离的影响程度相当，敏感指数分别为18.0%、16.3%、18.0%，孔隙率对重金属As(V)的迁移影响极弱。这与不考虑吸附作用或者分配系数较小的溶质运移规律存在较大偏差[17-18]，其原因是重金属As(V)溶质分配系数比较大，土壤对其的吸附能力比较强，导致在弥散度与孔隙率对重金属As(V)溶质迁移规律的影响上存在变化。一是重金属溶质迁移过程中，由于土壤的吸附作用被大量阻滞在土壤中，迁移距离受弥散作用影响较大，从而弥散度对重金属的敏感性增强；二是重金属溶质在土壤中的实际对流迁移以及弥散系数都直接受阻滞因子大小的影响，而孔隙率的变化对其影响极小，所以在本文的模型中，孔隙率对溶质的迁移作用影响可忽略。

表2 各参数不同幅度下重金属砷迁移距离拟合偏差Table 2 Fitting deviation of migration distance of heavy metal arsenic under different ranges of each parameter

同时可以注意到，各参数的变化对重金属As(V)在土壤中迁移距离的影响并不是呈线性变化，是因为其迁移距离受多种因素的影响，除了所探讨的参数外，还有水压差，降雨入渗系数等等。渗透系数及弥散度变化幅度在±25%以内，可近似看作线性变化。

4 结论

基于南京某废弃化工厂重金属As(V)在土壤中的迁移，运用有限元方法进行三维数值模拟，再在该模型基础上，改变渗透系数、弥散度、分配系数及孔隙率的大小，对比分析重金属溶质在土壤中的迁移规律，并引入敏感性指数，比较各参数对重金属溶质迁移影响的强度。得出如下结论。

(1)该污染场地在连续通入重金属As(V)污染物前7年，污染晕中心浓度持续上升到7.12 mg/L，污染范围随时间推移不断扩大，第7年纵向迁移距离38.5 m，远大于横向迁移距离23.7 m。停止通入污染物后，在地下水对流作用下，污染范围整体向下游移动，同时由于弥散作用，污染晕中心浓度持续下降，在第40年时污染晕中心往下游迁移了6.5 m，浓度下降到0.33 mg/L。

(2)重金属砷As(V)在土壤中的迁移距离随渗透系数、弥散度的增加而增加，随分配系数增加而减小。并且由于土壤对重金属As(V)的吸附能力强，在土壤阻滞作用下，重金属As(V)在土壤中实际迁移速度较小，分配系数对重金属的迁移影响较大，而孔隙率的变化却影响极小，几乎可忽略。

(3)通过对各参数的敏感性分析，得出在该模型中渗透系数、弥散度、分配系数对重金属As(V)在土壤中的迁移分布影响程度相当，敏感性指数分别为18.0%、16.3%、18.0%，孔隙度为0.006%。