刘 玲 ，陈 坚 ，牛浩博 ，李 璐 ，殷乐宜 ，魏亚强,3

（1.生态环境部环境规划院/长江经济带生态环境联合研究中心，北京 100012；2.中国地质大学(北京)水资源与环境学院，北京 100083；3.上海交通大学环境科学与工程学院，上海 200240）

污染物下渗及其迁移转化过程是地下水污染防治工作和地下水环境研究的焦点，已引起国内外学者的广泛关注[1]。已有学者将各土壤水运动模拟软件与地下水模拟软件MODFLOW相结合建立了土壤-地下水水流运动耦合模型[2-8]。Yakirevich等[9]采用有限差分法建立了拟三维饱和-非饱和水流运动和溶质运移耦合模型，有效模拟场地尺度的水流运动及溶质运移。林琳等[10]将三维Richards方程分成一维非饱和方程和三维饱和方程，并将饱和模型与非饱和模型在潜水面处耦合建立了基于迦辽金有限单元法的拟三维数学模型，为研究大范围饱和-非饱和水分运动提供了有效工具。查元源等[11-12]将一维非饱和水流模型同三维地下水模型结合，模拟了包含大气边界的区域地下水运动问题。并在2014年全面分析了在不同水流条件下各迭代数值模型的稳定性、质量守恒、计算成本以及适用范围。并釆用线性化的方法，将各种类型方程迭代模型转化为非迭代模型，提高了模型的稳定性和计算效率。

由于数据采集难度大并且模拟计算量受限，现有研究多是用非饱和带模型计算饱和带的补给流量，将其代入饱和带方程计算新的地下水位，更新非饱和带的下边界。这种方式不利于模型之间的计算反馈，易出现计算误差[13]。因此，针对大范围场地尺度的三维土壤-地下水水流及溶质运移模型的构建和应用显得尤为迫切。

FEFLOW软件基于变饱和的控制方程，可以实现完全三维的饱和-非饱和数值模拟，但是现有研究多是将FEFLOW应用于饱和带[14-16]。为了精确刻画实际场地土壤-地下水系统中污染物迁移规律，揭示变饱和反应溶质迁移模型的参数敏感性，本文以某铬污染场地为例，建立了完全三维土壤-地下水水流及溶质运移模型，采用统一的变饱和水流及溶质运移控制方程，以水头作为变量，实现了大范围场地的土壤-地下水耦合模拟。

1 研究区概况

1.1 研究区水文地质条件

研究区位于湘乡市丘陵地区（图1），多年平均降雨量为1 484 mm，平均蒸发量为1 142 mm。隶属华南加里东褶皱带，地壳运动包括雪峰、加里东、印支——燕山等。

图1 研究区位置图Fig.1 Location of the study area

研究区地下水主要为孔隙水和岩溶水，以及少量裂隙水（图2）。地下水补给来源有大气降水、地表水，孔隙水径流途径短，大部分地区就地排泄；裂隙水沿岩石裂隙、孔隙径流，向地势低洼的沟谷排泄；岩溶水沿地下管道、裂隙径流，向地势低洼的沟谷排泄。

图2 研究区水文地质平面图（本图来源为全国地质资料馆1∶20万水文地质图G4902幅）Fig.2 Hydrogeological map of the study area (The source of this figure is G4902 1∶200 000 hydrogeological figure of the national geological data center)

研究场区内地层由上至下共分为5层：①杂填土（Q4ml）：由粉质黏土、碎石构成，厚3～4 m，孔隙度较大，透水性佳；②粉质黏土（Q4al）：以黏土为主，含少量砂，厚2～4 m，透水性差，含水量较少，为相对阻水层；③中-粉砂（Q4al）：由砂和黏土构成，孔隙发育，厚1 m左右，含水量较少；④圆砾（Q4al）：由砂砾构成，厚3～4 m，孔隙水发育，以侧向补给为主；本层为模型的主要含水层；⑤泥质粉砂岩（Kdnd）：由砂泥构成，深度约9～10 m，弱透水，为隔水底板（图3）。

图3 场地水文地质剖面图Fig.3 Hydrogeological profiles of the site

1.2 场地污染特征

本研究结合省地方标准与场地的实际情况，采用重金属污染场地修复标准DB43/T1165-2016的工业用地标准作为土壤污染物浓度的参考标准值。按照调查技术规范要求，调查共布设297个土壤采样点、45个地下水采样点。场地污染情况如表1所示。

表1 六价铬检测结果Table 1 Statistics of soil test results

2 研究方法

2.1 数学模型构建

2.1.1 水流数学模型

以水头h作为主要变量，三维变饱和多孔介质中的水流运动控制方程为[17]：

式中：S0——贮水系数/m-1；

h——水力水头/m；

∇h——水力梯度，无量纲；

Sw—— 饱和度，非饱和带为水头h的 函数，0＜Sw≤1，饱和带Sw=1，无量纲；

krw—— 相对渗透率，非饱和带krw是饱和度的函数，由Van Genuchten-Mualem模型计算可得；饱和带为1，无量纲；

K——介质渗透系数/（m·s-1）；

t——时间/s；

V——达西流速/（m·s-1）；

φ——孔隙度，无量纲；

e——单位法向量，无量纲；

Q——流量系统的源汇项/s-1；

χ——浮力系数，表示流体的密度效应，无量纲；

∇V——流速梯度/（s-1）。

2.1.2 溶质运移数学模型

三维变饱和溶质运移控制方程为[18]：

式中：c——污染物源浓度/（mg·L-1）；

c*——渗入地下水中的污染物浓度/（mg·L-1）；

Dij——水动力弥散系数/（m2·s-1）；

D0——源溶液的分子扩散系数/（m2·s-1）；

|V|——达西流速的绝对值/（m·s-1）；

R——阻滞系数，表征土壤吸附能力，无量纲；

ρs——土壤颗粒密度/（kg·m-3）；

kd——分配系数/（L·kg-1）；

λ——反应系数/（10-4·s-1），本文中主要考虑土壤及地下水中铬的氧化还原反应；

αL、αT——纵向和横向弥散度/m；

τ——弯曲度，由Millington-Quirk方程可得，无量纲；

qc——边界浓度通量/（mg·m-2·s-1）；

ni——边界上的单位法向量，无量纲；

xi、xj—— 空间坐标（xi、xj=x,y,z）/m；

Vi、Vj——xi、xj方向的达西流速/（m·s-1）。

2.2 耦合方法及数值求解过程

2.2.1 耦合方法

由于水流和溶质运移控制方程均与流体密度相关密切，因此，可将水流和溶质运移通过表示流体密度效应的浮力项非线性耦合，利用FEFLOW软件实现土壤-地下水耦合数值模拟。

模型在非饱和带基于以下假设[19]：

（1）气液分离，唯一动态相是液相；

（2）流体运动符合达西定律，在模拟中不断根据瞬时条件改变潜水面的位置。

2.2.2 数值求解过程

采用Galerkin有限单元（FEM）网格及上游加权法对水流及溶质运移的控制方程进行求解，采用预处理共轭梯度法（PCG）求解水流控制方程，用预处理正交最小化（ORTHOMIN）法求解溶质运移控制方程[20]。

2.3 数值模型构建

2.3.1 概念模型

研究区面积约为11.16 km2，由于研究区不是完整的水文地质单元，因此上游没有完整的自然边界。根据水文地质调查结果，将涟水河西北侧约2.5 km处等水位线作为模拟区上游边界，下游以涟水河为边界。垂向上，以地表作为模拟的上边界，底部泥质粉砂岩为下边界，总深度约为23 m。根据研究场地渗透系数不同（表2），将模型概化为5层，由于第一层中碎石土主要存在于上部1 m范围内，导致第一层上部和下部孔隙度不同，因此将模型第一层细化为2个亚层，每层厚度由软件插值可得。区内地下水主要接受降雨及地表水的补给，流向为自西北向东南流向涟水河。

表2 研究区水文地质参数取值表Table 2 Values of hydrogeological parameters in the study area

污染物主要是由于上层渣土长期淋溶导致六价铬不断进入含水层所导致。六价铬和三价铬在运移过程中可能会发生相互转化。本次模拟将场地调查获取的浓度数据作为模拟的初始条件，在污染源全部清理后，无人工干预的情况下，模拟土壤和地下水中六价铬的时空变化趋势。

2.3.2 初始和边界条件

初始条件：模型初始水位为48 m，土壤含水量由模型计算可得。

垂直边界：模型上边界概化为降雨入渗、蒸发边界；本文中补给由年平均降雨量与蒸发量差值计算所得，设为343 mm/a；下边界为隔水边界。

侧向边界：上游边界概化为定水头边界，水头最大值为56 m，最小值为54 m；下游（涟水河）概化为定水头边界，水头最大值为45 m，最小值为37 m。

溶质边界：本次模拟中为自由出流边界。

2.3.3 模型参数选取

根据水文地质条件将研究区概化为非均质、各向异性三维非稳定渗流系统。通过土壤试验及相关文献获得土壤相关参数[21]，含水层的水文地质参数主要通过前期的水文地质调查获得（表2）。

溶质运移考虑对流、弥散、吸附、反应衰减等过程。弥散系数与土壤性质有关，通过试验由“三点公式”求得[22-23]。纵向弥散度为弥散系数与平均孔隙水流速度之比，横向弥散度取纵向弥散度的1/5。根据线性等温平衡吸附规律，利用试验确定六价铬的分配系数，阻滞系数为分配系数与介质密度的乘积（式4）。反应衰减主要考虑六价铬的氧化还原反应，反应系数在野外分析的基础上结合了前人的研究资料进行整理分析综合确定[24-29]（表3）。

表3 溶质运移模型参数取值表Table 3 Parameter values of the solute transport model

2.4 模型识别验证

模型参数校验主要采用“试错法”进行调整。初始流场是将研究区参数初始值输入模型，经过稳定流计算得到天然流场，然后根据实际观测水位对天然流场进行参数校正，得到校正后的地下水初始流场（图4），将模拟值与实际值做拟合（图5），得到其相关系数为0.999 6，模拟值与实际值相比，均方差为0.092 1，均方根误差为0.303 5，表明模型初始流场基本符合研究区实际水文地质条件，反映了实际流场特征，故可利用该模型得到的流场作为非稳定流的初始流场并以此为基础进行溶质运移模拟。

图4 研究区模拟流场图与实际流场图Fig.4 The simulated flow field and actual flow field in the study area

图5 模拟值与实际值拟合图Fig.5 Fitting diagram of the simulated value and actual value

3 结果和讨论

3.1 模拟预测结果

根据场地污染现状（图6），上层土壤污染晕中心浓度最高值为3 410 mg/L，下层土壤为3 430 mg/L（第一层顶底板），地下水（第四层底板）污染晕中心浓度最高值为109 mg/L，并分别在主要含水层底板位于场区周围区域、场地下游方向以及涟水河西侧河岸设置4处观测点监测浓度随时间的变化，分别为CG-1、CG-2、CG-3、CG-4。

图6 初始污染羽分布图Fig.6 Distribution figures of the initial contaminate plume

3.1.1 土壤模拟预测结果

土壤初始污染晕自二分厂向四周呈递减趋势。在降水淋滤作用下，土壤中六价铬不断向下运移，并向四周扩散。第28天时模型第三层下部六价铬浓度达到最大值（图7a），污染晕中心浓度约为3 382.94 mg/L。此后，土壤中六价铬浓度逐渐降低。第378天时土壤中污染晕的分布范围达到最大（图7b），场地东南侧迁移距离最大，约为300 m，北侧最短，约为90 m，此时，污染晕中心浓度为357.94 mg/L。而后污染晕范围逐渐减小，在1 900 d时土壤中污染晕已完全消失。

图7 土壤中（第三层）六价铬的运移模拟预测结果Fig.7 Simulation and prediction results of hexavalent chromium transport in soil (the third layer)

3.1.2 饱和地下水模拟预测结果

饱和带初始浓度中心为东、西渣场。地下水流动是六价铬在含水层运移的主要驱动力，结合研究区地下水流场可知污染物主要向污染场地东南侧的涟水河迁移（图8），第10天时运移至场地南侧边界GC-2附近，此时污染晕中心浓度值为53.01 mg/L，第25天时运移至东北侧边界GC-1附近，污染晕中心浓度值为422.45 mg/L，第49天时主要含水层六价铬浓度达到最大值，为1 533.65 mg/L。在此之前，由于土壤中六价铬不断下渗，东、西两渣场以及两场之间区域浓度均较高，而后随着入渗含水层中污染物浓度逐渐降低。在第585天时污染羽迁移至涟水河，污染晕中心浓度为25.79 mg/L，第917天时主要含水层污染晕分布范围最大，此时污染晕中心位于西渣场，浓度为7.91 mg/L，而后污染晕逐渐减小（表4）。第2 000 天时仅在场地边界处存在小范围浓度超标区域。

表4 污染晕面积表Table 4 Contaminant halo area

图8 地下水六价铬运移模拟预测结果Fig.8 Simulation and prediction results of hexavalent chromium migration in groundwater

土壤中六价铬在场地附近以垂向运移为主，水平运移速率较小；随着埋深的增加，水平运移逐渐加强（图9）。由于土壤中六价铬的下渗，地下水中的污染物浓度先是逐渐增大，当污染物全部进入含水层后，地下水中污染物浓度开始逐渐减小。土壤中的六价铬约需10 h下渗至潜水面处（图10a），全部污染物下渗至隔水底板约需2 d（图10b）。

图9 污染物迁移至涟水河时污染羽三维分布图Fig.9 3D distribution of the contaminant plume when contaminants migrate to the Lianshui River

图10 场地六价铬垂向浓度模拟结果Fig.10 Simulation results of chromium vertical concentrations at the site

3.2 参数敏感性分析

土壤是地下水污染的重要媒介，污染物主要是通过大气降水、地表水或灌溉水的入渗淋滤从土壤进入地下污染地下水。因此本文将溶质运移模型简化，考虑初始污染物仅存在于上层土壤的情况对参数敏感性进行分析。

已有研究[30]表明降雨量、阻滞系数以及反应常数会对污染物运移产生影响。本文分别改变降雨量、阻滞系数和反应常数研究其影响。

3.2.1 降雨量对污染物运移的影响

为了讨论降雨入渗变化对污染物运移的影响，本文将模型降雨入渗条件设置为时变的情况。依据研究区近年的降雨量和蒸发量统计资料，将丰水期（每年的7—8月）入渗补给量设为200 mm/a，枯水期（每年的9月至次年6月）设为15 mm/a。

结果显示，每一个丰水期过后约15 d，非饱和带土壤底部的污染物高浓度区域面积均会有所回升，上升约8 d后，含水层达到新的收支平衡，此后非饱和带中污染物浓度再次开始下降（图11）。主要是由于降雨首先要补给包气带，地下水对降雨补给的响应存在滞后性。当补给到达含水层之后，地下水位上升，非饱和带中污染物浓度随之增大，随着污染物的不断下渗以及降雨补给的持续淋洗，非饱和带和饱和带中的污染物浓度均逐渐降低。

图11 雨季后土壤中污染物浓度变化Fig.11 Changes in concentrations of pollutants in soil after rainy season

3.2.2 阻滞系数对溶质运移的影响

本文中阻滞系数是实验室测定的六价铬的分配系数与固体体积百分数的乘积，无量纲，与含水率密切相关。饱和含水率为定值，而非饱和带含水率变化较大，因此本文仅考虑非饱和带土壤阻滞系数对溶质运移的影响。

由图12、图13可以看出，当阻滞系数相差两个数量级时，曲线变化明显，污染物运移速度变缓，同一距离处污染物的浓度降低。阻滞系数的大小表征土壤颗粒对污染物吸附能力的强弱，阻滞系数越大，吸附作用对污染物的阻滞作用就越大，污染物迁移越慢。阻滞系数对污染物运移的影响可能由于研究区土壤中含有黏性土，黏性土颗粒高度分散，电荷不均衡，表面能较大，可以吸附水中的六价铬[31]。另外，在较大时空的模拟中，阻滞系数一般是由实验室测定的分配系数通过溶质运移控制方程中的阻滞系数计算公式求得，然而通过这样计算求得的阻滞系数值往往偏高，可能会导致模拟结果不够精准，因此在有条件的情况下，应尽可能采用试验直接测定阻滞系数的方式进行模拟[32]。

图12 不同阻滞系数下场区至河流浓度随距离的变化Fig.12 Variation of concentration with distance under different retardation coefficients from site to river

图13 不同阻滞系数河岸处观测点浓度随时间的变化Fig.13 Variation of concentration with time at the observation point on the bank with different retardation coefficients

3.2.3 反应常数对污染物运移的影响

文中反应常数是指六价铬被还原成三价铬的速率，单位为10-4s-1。由图14、图15可见，当反应常数增大两个数量级时，曲线变化及其明显，污染物运移速率变缓，同一距离处浓度显著降低。研究发现[33]，正常pH值条件下，在天然水体中，三价铬和六价铬之间可以相互转化，六价铬可被某些还原性物质还原为三价，三价铬也可被氧化成六价。六价铬极易溶于水，而三价铬在水溶液中则以难溶的沉淀形式存在。因此，当六价铬被还原成三价铬时，形成沉淀后析出，地下水中的六价铬浓度降低。

图14 不同反应常数下场区至河流浓度随距离的变化Fig.14 Variation of concentration with distance under different reaction constants from site to river

图15 不同反应常数河岸处观测点浓度随时间的变化Fig.15 Variation of concentration with time at the observation point on the river bank with different reaction constants

4 结论

（1）包气带的存在使地下水对降雨的响应存在滞后性。土壤和地下水中的浓度均呈现先升高后降低的趋势，场地边界处土壤中的污染物浓度比地下水中的污染物浓度低。在模拟期内的夏季降雨量增加期间地下水中污染物浓度呈现明显的下降趋势，之后缓慢回升；而土壤中的污染物浓度则在降雨量增大时呈现上升趋势，随后缓慢回落。降雨会引起潜水面波动，可能导致地下水中的污染物进入土壤，从而影响土壤中的污染分布。

（2）阻滞系数和反应常数的影响均较大，当变化两个数量级时，曲线变化明显；但当反应常数仅增大至10-6时，迁移出场区的污染物浓度约减少2 000 mg/L，较难迁移至涟水河，影响最为强烈。但是，在较大时空的模拟中，往往很难将阻滞系数和反应常数的影响分开，在今后的研究中需要探索更为有效的方式将二者分别讨论。

（3）采用基于FEFLOW的数值模型，能够解决各系统之间交互性差的问题，很好的模拟潜水面在土壤和地下水系统之间的波动，以及溶质在二者之间的运移，提供较为精确的模拟结果，实现了大范围场地的完全三维土壤-地下水耦合模拟，为提出科学有效的地下水污染防控措施提出定量化的理论依据。